# Unexpected behavior when using Cross

I'm trying to calculate the equation of the plane that is formed from the position of the Earth, Moon and Sun using precise coordinates from NASA's SPICE toolkit. I found this extremely clear refresher on plane geometry that lays out how to use the cross product of any two vectors between the three points.

Here is the example they use, which works just as expected:

P = { 1, -2, 0 };
Q = { 3, 1, 4 };
R = { 0, -1, 2 };

PQ = P - Q;
PR = R - P;
QR = R - Q;

normal = Cross[PQ, PR];
Print["Normal: ", normal];

g[{ x_, y_, z_ }] := normal[[1]]x + normal[[2]]y + normal[[3]]z
Print["Intercepts: ", { g[PQ], g[PR], g[QR] }];

v[{ p1_, p2_, p3_ }] := ( -normal[[1]]*p1 - normal[[2]]*p2 - normal[[3]]*p3 )
Print["Sums: ", { v[P], v[R], v[Q] }];


The output is just what we'd expect:

Normal: {-2,8,-5}
Intercepts: {0,0,0}
Sums: {18,18,18}


Okay, so far so good. But when I try to plug in the values for the Moon, Earth and Sun, I cannot get zeroes in the intercepts or consistent values from the right-hand value in the plane equation. (These figures are expressed as Cartesian coordinates in the J2000 reference frame, which is centered on Earth. They hold up under simulations correctly.)

This is the exact same code as above, but with realistic values.

ClearAll["Global*"]
P = { 96299004.0474634, -44374290.33702379, -102360951.2146616 }; (* Position of Sun *)
Q = { -195454.3807010595, 113678.2864233889, 279243.668852435 }; (* Position of Moon *)
R = { 0, 0, 0 }; (* Position of Earth *)

PQ = P - Q;
PR = R - P;
QR = R - Q;

normal = Cross[PQ, PR];
Print["Normal: ", normal];

g[{ x_, y_, z_ }] := normal[[1]]x + normal[[2]]y + normal[[3]]z
Print["Intercepts: ", { g[PQ], g[PR], g[QR] }];

v[{ p1_, p2_, p3_ }] := ( -normal[[1]]*p1 - normal[[2]]*p2 - normal[[3]]*p3 )
Print["Sums: ", { v[P], v[R], v[Q] }];


And the output is haywire!

Normal:{7.55022*10^11,6.88399*10^12,-2.27396*10^12}
Intercepts:{1.98574*10^7,-1.97919*10^7,65024.}
Sums:{-1.97919*10^7,0.,65024.}


For the life of me, I can't think of why this doesn't work, other than the possibility that we're losing precision given the extremely large values. I tried wrapping just about everything in N[expr,30] and got identical results.

Curiously, the value of 65024 is equal to 2^16 - 2^9.

## 3 Answers

It is simply a machine precision rounding issue:

P = Rationalize[{96299004.0474634, -44374290.33702379,    -102360951.2146616}, 0];
(*Position of Sun*)
Q = Rationalize[{-195454.3807010595, 113678.2864233889,
279243.668852435}, 0];
(*Position of Moon*)
R = {0, 0, 0};(*Position of Earth*)
PQ = P - Q;
PR = R - P;
QR = R - Q;
normal = Cross[PQ, PR];
g[{x_, y_, z_}] := normal[[1]] x + normal[[2]] y + normal[[3]] z
Print["Intercepts: ", {g[PQ], g[PR], g[QR]}];


Intercepts: {0,0,0}

incidentally, g is just a dot product so you can do this:

g[x_List] := normal.x


Wrapping everything in N[...,30] doesn't increase precision:

Precision[N[1.0, 30]]
(* MachinePrecision *)


Use SetPrecision or Rationalize.

george2079 has given you the answer you need. This is an addendum to his answer, offered with the intent of introducing you to methods for doing vector computations in a more efficient and concise way.

First I recommend you make a careful reading of the following documentation articles:

An experienced Mathematica user would perform your calculations using the tools and programming style advice discussed in the the linked articles. The resulting code would look something like this:

sun = Rationalize[{96299004.0474634, -44374290.33702379, -102360951.2146616}, 0];
moon = Rationalize[{-195454.3807010595, 113678.2864233889, 279243.668852435}, 0];
earth = {0, 0, 0};

sunMoon = sun - moon;
sunEarth = sun - earth;
moonEarth = moon - earth;

normal = Cross[sunMoon, sunEarth]

{-(1184253413193392152257837111260/1568501642905352793),
-(613314201971573899984964150516/89092826182611295),
561634409141053254820817385379/246985574148174535}

#.normal & /@ {sunMoon, sunEarth, moonEarth}


{0, 0, 0}

#.normal & /@ -{sun, earth, moon}


{0, 0, 0}`

I hope this will help you to advance your Mathematica skills. If you take the time to do so, you will enjoy it even more.

• Thank you so much! Extremely helpful. – Chris Wilson Jan 22 '18 at 18:29
• @ChrisWilson. Extremely helpful, but not worth an up-vote? – m_goldberg Jan 22 '18 at 23:17
• Completely my oversight. To much wrangling with Mathematica after 4AM :) Corrected. – Chris Wilson Jan 22 '18 at 23:57