# Plotting 3D Intersections

I'm trying to plot the intersections of four spheres. I've got the spheres plotted using:

ContourPlot3D[{(-9.2877 - x)^2 + (9.3049 - y)^2 + (5436354.04 - z)^2 == 21496269.296^2,
(20.40241 - x)^2 + (204.918 - y)^2 + (23272267.679 - z)^2 == 20095995.0541^2,
(-39.29329 - x)^2 + (282.248 - y)^2 + (20240909.994 - z)^2 == 22488938.185^2,
(8.341136 - x)^2 + (48.1826 - y)^2 + (23018246.984 - z)^2 == 19410828.319^2},
{x, -20000000, 20000000}, {y, -20000000, 30000000}, {z, -20000000, 40000000},
AxesLabel -> {"x", "y", "z"}]


Is there a simple way to find the intersections and plot those or highlight them?

• Possible duplicate: 82473 – shrx Jun 24 '15 at 13:30
• Four spheres in 3D will not in general have a common point of intersection, and in particular, yours don't. Do you want the set of intersections of all six pairs of spheres (a set of circles)? The set of intersections of all four triplets of spheres (a set of points)? Something else entirely? – Michael Seifert Jun 24 '15 at 13:30
• I thought they might not. I'm looking at graphing pseudoranges from GPS satellites. Ideally if the pseudoranges were the geometric ranges there would be a point of intersection. However, there are errors that leads to either more overlap or less overlap of the spheres. Getting the set of intersections of the pairs of spheres would be useful. – David M Jun 24 '15 at 13:35

Considering your GPS-inspired goal, here's an approach which visualizes a least-squares solution to distances from sphere shells, and the corrected standard deviation as size of the red sphere:

Module[{spheres, sol},
spheres = {
Sphere[{-9.2877, 9.3049, 5436354.04}, 21496269.296],
Sphere[{20.40241, 204.918, 23272267.679}, 20095995.0541],
Sphere[{-39.29329, 282.248, 20240909.994}, 22488938.185],
Sphere[{8.341136, 48.1826, 23018246.984}, 19410828.319]};
sol = Sphere[{x, y, z} /. #2, Sqrt[#1/(Length[spheres] - 1)]] & @@
Minimize[
Total[RegionDistance[#, {x, y, z}]^2 & /@ spheres], {x, y, z}, Reals];
Graphics3D[{
Opacity[0.5], spheres,
Opacity[1], Red, sol}]]


This solution requires v10 due to its use of geometric regions functionality.

• Wow that's great, thanks! I'm getting a couple of issues running 10.1: Minimize::ivar : 99999.985 is not a valid variable and ReplaceAll::reps : {100000.,100000.,100000.} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing – David M Jun 24 '15 at 14:50
• @DavidM Have you defined variables x, y or z earlier on your session? Restarting your kernel might help. – kirma Jun 24 '15 at 14:57
• Restarting the kernel fixed it. Thank you – David M Jun 24 '15 at 16:33

Do you mean the intersection of the regions enclosed by the spheres?

RegionPlot3D[And @@ {
(-9.2877 - x)^2 + (9.3049 - y)^2 + (5436354.04 - z)^2 <= 21496269.296^2,
(20.40241 - x)^2 + (204.918 - y)^2 + (23272267.679 - z)^2 <= 20095995.0541^2,
(-39.29329 - x)^2 + (282.248 - y)^2 + (20240909.994 - z)^2 <= 22488938.185^2,
(8.341136 - x)^2 + (48.1826 - y)^2 + (23018246.984 - z)^2 <= 19410828.319^2},
{x, -20000000, 20000000}, {y, -20000000, 30000000}, {z, -20000000, 40000000},
AxesLabel -> {"x", "y", "z"}]


• I think that will do quite well. For some reason though it's not showing up in Mathematica for me. I'm just getting an empty bounding box – David M Jun 24 '15 at 14:20

If you're actually looking for the pairwise intersections of the spheres, here's some code that does the trick. First, let's define the set of equations:

eqns = {(-9.2877 - x)^2 + (9.3049 - y)^2 + (5436354.04 - z)^2 == 21496269.296^2,
(20.40241 - x)^2 + (204.918 - y)^2 + (23272267.679 - z)^2 ==  20095995.0541^2,
(-39.29329 - x)^2 + (282.248 - y)^2 + (20240909.994 - z)^2 == 22488938.185^2,
(8.341136 - x)^2 + (48.1826 - y)^2 + (23018246.984 - z)^2 == 19410828.319^2};


We can use Reduce to find out whether these equations have a solution over the reals:

Reduce[eqns, Reals]

(* False *)


So no common solution to these over the real numbers. Similarly, we can try to see if the triplets of spheres have any common intersections:

Reduce[#, Reals]& /@ Subsets[eqns, {3}]

(* {False, False, False, False} *)


So no luck there, either. Let's go to the pairwise intersections of the spheres:

intersections = Reduce[#, Reals]& /@ Subsets[eqns, {2}];


This yields a set of inequalities & equalities, of which two are identically False (i.e., the spheres don't intersect at all) and the rest of which are non-trivial. So there are pair-wise intersections of these four spheres.

Actually plotting them is trickier. Mathematica is pretty crap at dealing with implicit curves in 3D, so we have to massage the data into a form where it can be fed into ParametricPlot3D. Basically, we discard the first and last points of each element of intersections (which are always single points); extract the valid regions of z, along with x and y as functions of z; and then feed these data into ParametricPlot3D. The code is pretty ugly, and I'm sure it could be improved:

intersections = Reduce[#, Reals, Backsubstitution -> True] & /@ Subsets[eqns, {2}];
paramcurves[inters0_] :=
Module[{inters = inters0, zranges, xyfns, plotdata},
inters = Drop[Drop[inters, 1], -1];
zranges = ({z, First[Minimize[{z, First[#]}, z]], First[Maximize[{z, First[#]}, z]]} & /@ inters) /. Or -> List;
xyfns = (({x, y} /. {ToRules[Drop[#, 1]]}) & /@ inters) /. Or -> List;
plotdata = Table[{xyfns[[i]], zranges[[i]]}, {i, 1, Length[zranges]}];
ParametricPlot3D[Append[Evaluate[#[[1, 1]]], z], Evaluate[#[[2]]],  PlotStyle -> {Thick, Red}] & /@ plotdata
]
intersections = DeleteCases[intersections, False];
curves = paramcurves /@ intersections;
Show[contplot, Flatten[curves], PlotRange -> All]


• "pretty crap at dealing with implicit curves in 3D" - but have you seen this? – J. M. is away Jun 24 '15 at 15:00
• I had not in fact seen that, and it's pretty useful; the top-rated answer there could easily be adapted to the present case. – Michael Seifert Jun 24 '15 at 15:20