Modeling the Solar System

Using Kepler's First and Second Laws

An honor's project, Science Seminar presentation, and AMATYC/ORMATYC Math Conference Presentation

Work by Ali, Darrel, David, Dusty, Jonathan, and Tran
a.k.a Stardust and the Kepler 5

Our Beautiful Mathematica Algorithm

Planatary Data

Build the Data or

Null

Use the Prebuilt Data

PData={{"Mercury", 4879.4`5., 4.60012719234314290031325072`7.8999879594440605*^7, 
   6.9817079426097492`7.931569737795803*^7, 0.1222580451741751943`6., 
   0.0001745329251994329`1., 87.9090454999999999962`7., 1407.5087999999999818909`6., 
   3.9851394280771610283`7.674334351530899, 1.3518700794063618915`7., Null, Null, 58.98, 
   38.86, 5.48772863676997}, {"Venus", 12103.7`5., 
   1.0747600156617407`7.774163177215992*^8, 1.08941849447193341`7.776538793655386*^8, 
   0.0592488666503767056`6., 3.0955159613371429376`5., 224.5469999000000000002`8., 
   -5832.24`5., 1.3383305132010906804`7.+π, 2.2956835759598511235`8., Null, Null, 35.26, 
   34.79, 4.135390983656999}, {"Earth", 12734.94`5.297422546827782, 
   1.470980735410081476042978466`8.931464628195396*^8, 
   1.520977007525234044156968593`8.93354146590687*^8, 8.726646259971`1.*^-7, 
   0.4092797095926702858`4., 365.0063509999999999961`8., 23.93447232`8., 
   2.9450574097080189643`8.175653323165909, 1.7967674211761812588`8., Null, Null, 30.29, 
   29.29, 5.731722873671973}, {"Mars", 6772.`4., 
   2.066445448522859580926597118`8.692210361367064*^8, 
   2.492287296039747543272096664`8.73158132830309*^8, 0.0322992376703322625`6., 
   0.4396484385773716204`4., 686.5093739999999999994`8., 24.622962`9., 
   0.8653087613317094869`7., 5.8650190791674618627`8., Null, Null, 26.5, 21.97, 
   4.580229856896876}, {"Jupiter", 138346.`5., 7.40742598178938399`8.821156522347685*^8, 
   8.16081455277686681`8.834877674399578*^8, 0.022781782726281984`6., 
   0.0546288055874225158`3., 4329.8544749999999999953`8., 9.92496`5., 
   1.7550359006292960409`8., 0.2575032598453653944`7., Null, Null, 13.72, 12.44, 
   5.6297310681731805}, {"Saturn", 114632.`5., 
   1.349467375380085103`7.2346909247874365*^9, 
   1.503983449624582877`7.280917915189853*^9, 0.0433620071340983205`6., 
   0.4665265090580842959`4., 10748.3367699999999999989`8., 10.5`4., 
   1.9847018570370526761`8., 1.613241687005847518`7., Null, Null, 10.18, 9.09, 
   0.3657789487830384}, {"Uranus", 50532.`5., 2.735555034715551676`9.224734532786503*^9, 
   3.006389404878433276`9.258508596104893*^9, 0.0134365917794035456`5., 
   1.7064084096748560273`4., 30666.1487900000000000094`8., -15.6`2., 
   4.4371484629500816902`7.534647833994085, 2.9838888911697904043`8., Null, Null, 7.11, 
   6.49, 5.291657999535076}, {"Neptune", 49105.`5., 
   4.459631495776823041`8.057435137957565*^9, 4.536874325210059688`8.064806205642425*^9, 
   0.0308778415275080805`6., 0.4942772441647941362`4., 60148.8318000000000000045`8., 
   18.431999999999987065`3., 2.2989771867912896229`8., 0.7848981265667519342`7., Null, 
   Null, 5.5, 5.37, 5.233615853529528}, {"Pluto", 2346.`4., 
   4.436824612745022184`8.466842102737218*^9, 7.375927931417028733`8.679155718017864*^9, 
   0.2991799770537379902`7., 2.1385519324686520305`5., 90491.0372499999999999946`8., 
   -153.2807999999999689559`5., 5.0667513814645829651`8.420263052008503, 
   3.9107027062759296667`8., Null, Null, 6.1, 3.71, 4.114885982939944}};
PlanetTextures={,,,,,,,,};

Modeling the Anomaly as a function of time

(****************************************************)
(*The AnomalyModule is a subprogram that calculates
the positional angle theta (anomaly) as a function of 
time using Kepler's second law.  It does this using 
numerical integration to create (t,θ) data points and 
then uses finds a regression model to fit the data.*)  
(*
Note: Changing the option to "1" for an individual
planet will display the error in the model rather 
than the model itself.*)
(****************************************************)
AnomalyModule[aphelion_,perihelion_,period_,planet_,option_:0]:=
Module[{n,AREA,A0,Area,A,B,RPS,TotalArea,e,r,θTimeTable,BestθModel,θModel},

(**********Set up*********)
Steps=1000;
n=10; (*Number of terms in the Fourier Series*)
A[aph_,per_]=((aph+per)/2); (*Major axis of the ellipse*)
B[aph_,per_]=(Sqrt[aph*per]);(*Minor axis of the ellipse*)
RPS=((N[2π])/Steps); (*Radians traveled in an average planetary step*)
TotalArea=A[aphelion,perihelion]*B[aphelion,perihelion]*π(*the area of the ellipse bounded by the orbit*);
e=(aphelion-perihelion)/(aphelion+perihelion);(*the eccentricity of the orbit*)
r[θ_]=A[aphelion,perihelion](1-e^2)/(1+e*Cos[θ]);(*the polar equation of the ellipse with one foci at the sun*)
AREA[t_]=Chop[Simplify[∫(1/2)(((r[t]))^2)\[DifferentialD]t]];
Area[t_]=Piecewise[{{AREA[t],t<π},{AREA[t]+2AREA[3.14159],t>π}}];
PData[[planet,12]]=PData[[planet,7]]/TotalArea*Area[PData[[planet,15]]];
(**********Construct a model for θ(t)*********)
θTimeTable=Table[{(Steps/TotalArea) *Area[θMax],θMax},{θMax,0,2π,RPS}];(*Kepler's second law*)

(*Series expansion for θ(t)*)
BestθModel=FindFit[θTimeTable, RPS t+(∑_(i=1)^n)((Subscript[a, i])Sin[i RPS t]),Table[(Subscript[a, i]),{i,1,n}],t,MaxIterations->1000];(*Series expansion for θ(t)*)

θModel[t_]= RPS t+(∑_(i=1)^n)((Subscript[a, i])Sin[i RPS t])/.BestθModel;
(*Equation for θ(t)*)

(**********Display Error of option ≠ 0*********)
If[option==0,
θModel[(Steps/period) t],
GraphicsRow[{
(*Time vs θ*)
Show[ListPlot[θTimeTable],Plot[θModel[A],{A,0,period}],AxesLabel->{"Time","θ"}],
(*Time vs θ Error*)
Show[
ListPlot[Table[{θTimeTable[[i,1]],θModel[θTimeTable[[i,1]]]-θTimeTable[[i,2]]},{i,1,Length[θTimeTable],1}]],
AxesLabel->{"Time","θ Error"}],
(*Time vs Log of θ Error*)
Show[ListPlot[Table[{θTimeTable[[i,1]],Log[10,Abs[θModel[θTimeTable[[i,1]]]-θTimeTable[[i,2]]]]},{i,1,Length[θTimeTable],1}]],
Plot[Log[10,Abs[θModel[f]-2*π/period*f]],{f,0,period}],GridLines -> Automatic,AxesLabel->{"Time","Log of θ Error"}]
}]]
]

(****************************************************)
(*Calculates the anomaly models θ(t) for each of
the planets/moons given in the data set.*)
(*These models can be evaluated using the code: 
PData[[Row,11]]/.t->value*)
(****************************************************)
For[i=1,i≤9,i++,PData[[i,11]]=AnomalyModule[PData[[i,4]],PData[[i,3]],Floor[PData[[i,7]]],i,0]]

(****************************************************)
(*In order to check the anomaly models, we use them
to calculate the velocity of the planet at its
perihelion and aphelion and compare these results to
the experimental results published by NASA*)
(****************************************************)
Off[NIntegrate::"precw"];
VelocityCheck=Module[{A,aphelion,perihelion,e,r,VelocityTable},
VelocityTable={{"Planet","NASA min velocity", "Model min velocity", "NASA max velocity", "Model max velocity"}};
A[aph_,per_]=((aph+per)/2);
For[i=1,i≤9,i++,
aphelion=PData[[i,4]];
perihelion=PData[[i,3]];
e=(aphelion-perihelion)/(aphelion+perihelion);(*the eccentricity of the orbit*)
r[θ_]=((A[aphelion,perihelion](1-e^2))/(1+e*Cos[θ]));
AppendTo[VelocityTable,{PData[[i,1]],PData[[i,14]],(1/60)NIntegrate[(Sqrt[(((r[θ]))^2)+(((r'[θ]))^2)]),{θ,PData[[i,11]]/.t->((PData[[i,7]])/2)-(1/(48*60)),PData[[i,11]]/.t->((PData[[i,7]])/2)+(1/(48*60))}],PData[[i,13]],(1/(60^2))NIntegrate[(Sqrt[(((r[θ]))^2)+(((r'[θ]))^2)]),{θ,PData[[i,11]]/.t->0-(1/48),PData[[i,11]]/.t->0+(1/48)}],PData[[i,12]]}];
];
Print[VelocityTable//TableForm];
];

Setting our Setings

RangeSettings={{(10^8) {{-3,3},{-3,3},{-0.5,0.5}},{1,2,3,4},1250,(10^7) 2},{(10^8) {{-18,18},{-18,18},{-2.5,2.5}},{4,5,6},1000,(10^7) 15},{(10^8) {{-50,70},{-70,70},{-20,40}},{7,8,9},9000,(10^7) 60},{(10^8) {{-50,70},{-70,70},{-20,40}},{5,6,7,8,9},3000,(10^7) 45}};

Building the Equations

A[aph_,per_]=((aph+per)/2);
B[aph_,per_]=(1/2) (Sqrt[(((aph+per))^2)-(((aph-per))^2)]);
Shift[aph_,per_]=((aph-per)/2);
Tx[θ_]=N[{{1,0,0},{0,Cos[θ],-Sin[θ]},{0,Sin[θ],Cos[θ]}}];
Ty[θ_]=N[{{Cos[θ],0,Sin[θ]},{0,1,0},{-Sin[θ],0,Cos[θ]}}];
Tz[θ_]=N[{{Cos[θ],-Sin[θ],0},{Sin[θ],Cos[θ],0},{0,0,1}}];
phelionOrbit[t_,aph_,per_]={A[aph,per] Cos[t]-Shift[aph,per],B[aph,per] Sin[t]};
Planet[θ_,ϕ_,r_]={r Sin[ϕ] Cos[θ],r Sin[ϕ] Sin[θ],r Cos[ϕ]};
e={1,0,0};
OrbitBases=Table[Orthogonalize[{Tz[PData[[i,10]]].Ty[PData[[i,5]]].e,Tz[PData[[i,9]]].e}],{i,1,9}];
Orbit[θ_]=Table[phelionOrbit[θ,PData[[i,4]],PData[[i,3]]].OrbitBases[[i]],{i,1,9}];

Plotting the Planets

PositionXYZ={0,0,0};
Manipulate[
If[False,{time,reset}={4080,False}];RangeOption=If[view=="Inner",1,If[view=="Giant",2,If[view=="Outer",3,4]]];Module[{Plotrange,PlanetSelection,PlanetCount,Anomaly,Rotation,Planets},Plotrange=RangeSettings[[RangeOption,1]];PlanetSelection=RangeSettings[[RangeOption,2]];PlanetCount=Length[PlanetSelection];Anomaly=Table[PData[[i,11]]/. t->PData[[i,12]]+time+876,{i,PlanetSelection}];Rotation=Table[((24 2 π time)/(PData[[i,8]])),{i,PlanetSelection}];PositionXYZ=Table[Orbit[Anomaly[[i]]][[PlanetSelection[[i]]]],{i,1,PlanetCount}];Planets[θ_,ϕ_]=Table[RangeSettings[[RangeOption,3]] Tx[PData[[PlanetSelection[[i]],6]]].Planet[θ+Rotation[[i]],ϕ,PData[[PlanetSelection[[i]],2]]]+PositionXYZ[[i]],{i,1,PlanetCount}];Show[ParametricPlot3D[Orbit[θ],{θ,0,2 π},PlotStyle->Directive[Dashing[0.02],White]],Table[ParametricPlot3D[Planets[θ,ϕ][[i]],{θ,0,2 π},{ϕ,π,0},PlotStyle->Directive[Texture[PlanetTextures[[PlanetSelection[[i]]]]]],TextureCoordinateFunction->({#4,-#5}&),Mesh->False,PlotPoints->10],{i,1,PlanetCount}],Graphics3D[Flatten[{Glow[Yellow],Sphere[{0,0,0},RangeSettings[[RangeOption,4]]],Table[Text[Style[PData[[RangeSettings[[RangeOption,2,i]],1]],Red],PositionXYZ[[i]],{0,2}],{i,PlanetCount}]}]],
PlotRange->Plotrange,Ticks->False,Background->Black,Boxed->False,Lighting->{{"Point",White,{0,0,0}}}
]],
Style["Model of the Solar System",Bold,12],
{{time,4080,"days"}, If[view=="Inner",4000,0], If[view=="Inner",5000,If[view=="Giant",30000,250000]],0.1,ImageSize->Small,Animator,DefaultDuration->Dynamic[speed],AnimationRunning->False},{{speed,50,""},{500->"Slow",50->"Medium",25->"Fast"}},
Delimiter,
{{view,"Inner","Planets"},{"Inner","Giant","Outer","All"}},
Style["",Plain,8],
Style["",Plain,8],
Style["Planet Key",Bold,8],
ControlPlacement->Left]