# How to add moon to orbit Earth that Orbit Sun

I want to add moon(Blue) that orbiting earth(Red). while Earth is orbiting Sun. this is my code so far that consist of earth, sun and moon. how to write a code so that Moon is orbiting earth instead of sun? Tmoon is formula of period of the moon and earth. T is Period for sun and Earth. e is eccentricity. a is Semimajor.

(*"Measuring the Earth-Sun system";*)

h = 0; k = 0; G = 6.63*10^(-11);
e = 0.0167086;
a = 1.496*10^(11); b = a*Sqrt[1 - e^2];
MS = 1.989*10^30(*mass Sun*) ;
ME = 5.9723*10^(24)(*Mass Earth*);
MM = 0.07346*10^(24)(*Mass Moon*);
T = 2 Pi*Sqrt[a^3/(G*(MS + ME))](*Formula Period and semimajor for circulating around central object*);
Tmoon = 2 Pi*Sqrt[aMoon^3/(G*(MM + ME))];
omega0 = 2 Pi/T;
omega1 = 2 Pi/Tmoon;
aMoon = (3.844*10^(8))*10^(2);
bMoon = aMoon*Sqrt[1 - eMoon^2];
eMoon = 0.0549;

(*** define linear eccentricity f *)
f1 = Sqrt[a^2 - b^2]; f2 = -f1;
f1Moon = Sqrt[aMoon^2 - bMoon^2]; f2Moon = -f1Moon;
(*f1=f1+h;f2=h+f2;*)
F1 = {f1 + h, k}; F2 = {f2 + h, k};
F1Moon = {f1Moon + h, k}; F2Moon = {f2Moon + h, k};
(*** end of define focal points*)

(* Define the end points of major and minor axes *)
esMl = {h - a, k}; esMr = {a + h, k};
esmu = {h, b + k}; esmd = {h, k - b};
lsM = Graphics[Line[{esMl, esMr}]];
lsm = Graphics[Line[{esmu, esmd}]];

(**** for the E-Sun ****)
x[t_] := a*Cos[omega0*t] + h;
xMoon[t_] := aMoon*Cos[omega1*t] + h;
y[t_] := b*Sin[omega0*t] + k;
yMoon[t_] := bMoon*Sin[omega1*t] + k;
prp1 = ParametricPlot[{x[t], y[t]}, {t, 0, T},
PlotRange -> {{h - a, h + a}, {k - b, k + b}}];
prp1Moon =
ParametricPlot[{xMoon[t], yMoon[t]}, {t, 0, T},
PlotRange -> {{h - aMoon, h + aMoon}, {k - bMoon, k + bMoon}}];
pp2 = Graphics[{PointSize -> 0.03, Point[{h, k}]}];
pp2Moon = Graphics[{PointSize -> 0.03, Point[{h, k}]}];
ppf = Graphics[{PointSize -> 0.02, Point[{F1, F2}]}];
ppfMoon = Graphics[{PointSize -> 0.02, Point[{F1Moon, F2Moon}]}];
(**** for the E-Sun ****)

(*"Simulate the motion of a moving point overlapping on the ellipse \
using Manipulate";*)
Manipulate[
gp = Graphics[{Red, PointSize -> 0.02, Point[{x[t], y[t]}]}];
gpMoon =
Graphics[{Blue, PointSize -> 0.02, Point[{xMoon[t], yMoon[t]}]}];
Show[{prp1, prp1Moon, ppf, ppfMoon, pp2, pp2Moon, lsM, gp, gpMoon}]
(*Show[{prp1,ppf,pp2,lsM,gp}]*)
,
{t, 0.0, 10 T, 0.005 T}
]

All help is much appreciated

• If you don't want to use PlanetaryMoonData[] or PlanetData[], do you at least have a reference for computing the heliocentric coordinates of the Moon? Commented Mar 24, 2016 at 15:17
• thanks for the info. already add the necessary reference Commented Mar 24, 2016 at 15:33

h = 0; k = 0; G = 6.63*10^(-11);
e = 0.0167086;(*eccentricity Earth*)
a = 1.496*10^(11)(*semimajor Earth*);
MS = 1.989*10^30 ;(*Mass Sun*)
ME = 5.9723*10^(24);(*mass Earth*)
T = 2 Pi*Sqrt[a^3/(G*(MS + ME))];
omega0 = 2 Pi/T;

MM = 0.07346*10^(24);(*Mass Moon*)
aMoon = 3.844*10^(8)*(10^(2));(*semimajor Moon*)
(*bMoon=aMoon*Sqrt[1-eMoon^2];*)
eMoon = 0.0549;(*eccentricity Moon*)
Tmoon = 2 Pi*Sqrt[aMoon^3/(G*(MM + ME))];
omega3 = 2 Pi/Tmoon;

f1 = e*a; f2 = -f1;
b = Sqrt[a^2 - f1^2];

f1 = f1 + h; f2 = h + f2;
F1 = {f1, k}; F2 = {f2, k};

x[t_] := a*Sin[-omega0*t] + h;

y[t_] := b*Cos[-omega0*t] + k;

prp1 = ParametricPlot[{x[t], y[t]}, {t, 0, T}];

pp2 = Graphics[{PlotRange -> {{h - a - aMoon,
h + a + aMoon}, {k - b - aMoon, k + b + bMoon}}, Brown,
PointSize -> 0.03, Point[{h, k}]}];

ppf = Graphics[{PointSize -> 0.02, Point[{F1, F2}]}];

f5 = eMoon*aMoon; f6 = -f5;
bMoon = Sqrt[aMoon^2 - f5^2]; f6 = -f5;
f5 = f5 + x[t]; f6 = x[t] + f6;
F5 = {f5, y[t]}; F6 = {f6, y[t]};

x2[t_] := aMoon*Cos[omega3*t] + x[t];
y2[t_] := bMoon*Sin[omega3*t] + y[t];

Manipulate[
gp = Graphics[{{Red, PointSize -> 0.02, Point[{x[t], y[t]}]}, {Blue,
PointSize -> 0.02, Point[{x2[t], y2[t]}]}}, Frame -> True,
PlotRange -> {{h - a - aMoon, h + a + aMoon}, {k - b - aMoon,
k + b + bMoon}}];

Show[prp1, ppf, pp2, gp]

,
{t, 0.0, T, 0.01 T}
]