108 double xba, yba, zba, xca, yca, zca;
109 double balength, calength;
110 double xcrossbc, ycrossbc, zcrossbc;
112 double xcirca, ycirca, zcirca;
122 balength = xba * xba + yba * yba + zba * zba;
123 calength = xca * xca + yca * yca + zca * zca;
130 xcrossbc = orient2d(b[1], b[2],
c[1],
c[2],
a[1],
a[2]);
131 ycrossbc = orient2d(b[2], b[0],
c[2],
c[0],
a[2],
a[0]);
132 zcrossbc = orient2d(b[0], b[1],
c[0],
c[1],
a[0],
a[1]);
135 xcrossbc = yba * zca - yca * zba;
136 ycrossbc = zba * xca - zca * xba;
137 zcrossbc = xba * yca - xca * yba;
141 denominator = 0.5 / (xcrossbc * xcrossbc + ycrossbc * ycrossbc +
142 zcrossbc * zcrossbc);
145 xcirca = ((balength * yca - calength * yba) * zcrossbc -
146 (balength * zca - calength * zba) * ycrossbc) * denominator;
147 ycirca = ((balength * zca - calength * zba) * xcrossbc -
148 (balength * xca - calength * xba) * zcrossbc) * denominator;
149 zcirca = ((balength * xca - calength * xba) * ycrossbc -
150 (balength * yca - calength * yba) * xcrossbc) * denominator;
151 circumcenter[0] = xcirca;
152 circumcenter[1] = ycirca;
153 circumcenter[2] = zcirca;
155 if (xi != (
double *) NULL) {
165 if (((xcrossbc >= ycrossbc) ^ (-xcrossbc > ycrossbc)) &&
166 ((xcrossbc >= zcrossbc) ^ (-xcrossbc > zcrossbc))) {
167 *xi = (ycirca * zca - zcirca * yca) / xcrossbc;
168 *
eta = (zcirca * yba - ycirca * zba) / xcrossbc;
169 }
else if ((ycrossbc >= zcrossbc) ^ (-ycrossbc > zcrossbc)) {
170 *xi = (zcirca * xca - xcirca * zca) / ycrossbc;
171 *
eta = (xcirca * zba - zcirca * xba) / ycrossbc;
173 *xi = (xcirca * yca - ycirca * xca) / zcrossbc;
174 *
eta = (ycirca * xba - xcirca * yba) / zcrossbc;