Re: Sphere formula
- To: mathgroup at smc.vnet.net
- Subject: [mg109472] Re: Sphere formula
- From: "David Park" <djmpark at comcast.net>
- Date: Wed, 28 Apr 2010 02:01:14 -0400 (EDT)
I don't know if you are demanding an analytic formula, or if you will settle
for numerical routine. An analytical formula is likely to be quite long and
not exceedingly transparent. If you will accept a routine, here is a
solution using the Presentations package.
We should be able to work in the plane of the points, calculate the center,
and then transform that center back to its 3D position. Generate 3 random
points for the sphere to pass through.
{sp1, sp2, sp3} = RandomReal[{0, 10}, {3, 3}]
To transform to a 2D case:
1) translate sp1 to the origin.
2) find the normal vector to sp2-sp1 and sp3-sp1 and rotate to the z axis.
3) drop the z component and find the circumcenter.
The triangle center in the xy-plane can be found using the
triangleCircumcenter routine from Presentations. Here we can get a simple
formula and you could use that in your own routine.
Needs["Presentations`Master`"]
triangleCircumcenter[{{0, 0}, {x2, y2}, {x3, y3}}] // First
{x2/2 + y2/2 + (y2 (-x2 x3 + x3^2 - x3 y2 + x2 y3 - y2 y3 + y3^2))/(
2 (x3 y2 - x2 y3)), -(x2/2) + y2/2 - (
x2 (-x2 x3 + x3^2 - x3 y2 + x2 y3 - y2 y3 + y3^2))/(
2 (x3 y2 - x2 y3))}
4) add the zero z component and rotate back to the normal.
5) translate the origin back to sp1. That should give us the 3D center.
The following routine calculates the center. The only Presentations routine
it uses is triangleCircumcenter and you could extract the x2,y2,x3,y3
coordinates (x1,y1 are zero) and use the expression above instead.
ClearAll[sphereCenter]
sphereCenter::usage =
"sphereCenter[{p1, p2, p3}] calculates the center of a sphere that \
passes through the three points p1, p2, p3 and lies in the plane of \
the three points.";
SyntaxInformation[sphereCenter] = {"ArgumentsPattern" -> {{_, _, _}}};
sphereCenter[{p1_, p2_, p3_}] :=
Module[{normal, q1, q2, q3, center},
(* Translate p1 to the origin = q1 *)
{q1, q2, q3} = TranslationTransform[-p1]@{p1, p2, p3};
(* Calculate a normal vector *)
normal = q2\[Cross]q3;
(* Rotate the normal vector to the z vector *)
{q1, q2, q3} =
RotationTransform[{normal, {0, 0, 1}}]@{q1, q2, q3} // Chop;
(* Find the center in 2D by dropping z components *)
center = First@triangleCircumcenter[Drop[#, -1] & /@ {q1, q2, q3}];
(* Pad and inverse rotate *)
center =
RotationTransform[{{0, 0, 1}, normal}]@PadRight[center, 3] // Chop;
(* Translate the origin back to p1*)
center = TranslationTransform[p1]@center]
The following generates three random points and then uses sphereCenter to
calculate the sphere center.
{sp1, sp2, sp3} = RandomReal[{0, 10}, {3, 3}]
spcenter = sphereCenter[{sp1, sp2, sp3}]
{{9.5243, 7.49836, 6.09033}, {6.7254, 8.46593, 9.04894}, {4.23268,
3.10471, 2.33983}}
{5.41194, 5.78047, 5.7232}
We can then plot the resulting sphere (using a Riemann sphere with its
annotations to give orientation), an outlined disk passing through the three
points, and mark the three points and the center of the sphere.
Module[{radius = Norm[sp1 - spcenter]},
Draw3DItems[
{(* To give orientation,
represent the sphere as a Riemann sphere scaled and translated to \
position. *)
ColoredRiemannSphere[Opacity[.3, Orange]] //
ScaleOp[radius, {0, 0, 0}] // TranslateOp[spcenter],
(* Draw a disk and circle through the 3 points. *)
{Opacity[.7, ColorData["Legacy"]["DodgerBlue"]],
Disk3D[spcenter, (sp2 - sp1)\[Cross](sp3 - sp1), radius,
Mesh -> None],
Black, Circle3D[spcenter, (sp2 - sp1)\[Cross](sp3 - sp1), radius]},
(* Mark the 3 points and the center *)
AbsolutePointSize[6], Point[{sp1, sp2, sp3, spcenter}]},
ViewRiemann,
ImageSize -> 300]
]
Peter Lindsay from the St. Andrews University Mathematics department keeps
an archive of Presentations solutions to MathGroup questions. These contain
downloadable PDFs and Mathematica notebooks. This should appear there is a
few days.
http://blackbook.mcs.st-and.ac.uk/~Peter/djmpark/html/
David Park
djmpark at comcast.net
http://home.comcast.net/~djmpark/
From: S. B. Gray [mailto:stevebg at ROADRUNNER.COM]
1. The center of a sphere through 4 points has a very nice determinant
form. (http://mathworld.wolfram.com/Sphere.html) What I want is a nice
formula for the center of a sphere through 3 points, where the center is
in the plane of the three points. I have a formula but it's a horrible
mess of hundreds of lines, even after FullSimplify.
2. (Unlikely) Is there a way to get Mathematica to put a long formula into a
matrix/determinant form if there is a nice one?
Any tips will be appreciated.
Steve Gray