35 double circumcenter[3];
40 double xba, yba, zba, xca, yca, zca, xda, yda, zda;
41 double balength, calength, dalength;
42 double xcrosscd, ycrosscd, zcrosscd;
43 double xcrossdb, ycrossdb, zcrossdb;
44 double xcrossbc, ycrossbc, zcrossbc;
46 double xcirca, ycirca, zcirca;
59 balength = xba * xba + yba * yba + zba * zba;
60 calength = xca * xca + yca * yca + zca * zca;
61 dalength = xda * xda + yda * yda + zda * zda;
63 xcrosscd = yca * zda - yda * zca;
64 ycrosscd = zca * xda - zda * xca;
65 zcrosscd = xca * yda - xda * yca;
66 xcrossdb = yda * zba - yba * zda;
67 ycrossdb = zda * xba - zba * xda;
68 zcrossdb = xda * yba - xba * yda;
69 xcrossbc = yba * zca - yca * zba;
70 ycrossbc = zba * xca - zca * xba;
71 zcrossbc = xba * yca - xca * yba;
78 denominator = 0.5 / orient3d(b,
c,
d,
a);
81 denominator = 0.5 / (xba * xcrosscd + yba * ycrosscd + zba * zcrosscd);
85 xcirca = (balength * xcrosscd + calength * xcrossdb + dalength * xcrossbc) *
87 ycirca = (balength * ycrosscd + calength * ycrossdb + dalength * ycrossbc) *
89 zcirca = (balength * zcrosscd + calength * zcrossdb + dalength * zcrossbc) *
91 circumcenter[0] = xcirca;
92 circumcenter[1] = ycirca;
93 circumcenter[2] = zcirca;
95 if (xi != (
double *) NULL) {
101 *xi = (xcirca * xcrosscd + ycirca * ycrosscd + zcirca * zcrosscd) *
103 *
eta = (xcirca * xcrossdb + ycirca * ycrossdb + zcirca * zcrossdb) *
105 *
zeta = (xcirca * xcrossbc + ycirca * ycrossbc + zcirca * zcrossbc) *