Modeling the Solar System

Using Kepler's First and Second Laws

An honor's project and Science Seminar presentation

Work by Ali, Darrel, David, Dusty, Jonathan, and Tran

Our Beautiful Mathematica Algorithm

Planatary Data

Data: http://nssdc.gsfc.nasa.gov/planetary/

PlanetaryData={{"Planet","Mass (10^24kg)","Diameter (km)","Density (kg/m^3)","Gravity (m/s^2)","Escape Velocity (km/s)","Rotation Period (hours)","Length of Day (hours)","Distance from Sun (10^6 km)","Perihelion (10^6 km)","Aphelion (10^6 km)","Orbital Period (days)","Max Orbital Velocity (km/s)","Min Orbital Velocity (km/s)","Orbital Inclination (degrees)","Orbital Eccentricity","Axial Tilt (degrees)","θ(t) Velocity Model","Ascending Node", "Perihelion", "Planet θ on 8th August 1997", "Time Adjustment"},{"MERCURY",0.33,4879,5427,3.7,4.3,1407.6,4222.6,57.9,46,69.8,88,58.98,38.86,7,0.205,0.01,"calculate",0.8435856953296872`,1.3518273188396879`,5.48772863676997`,"calculate"},
{"VENUS",4.87,12104,5243,8.9,10.4,-5832.5,2802,108.2,107.5,108.9,224.7,35.26,34.79,3.4,0.007,177.4,"calculate",1.3384738047326794`,2.2996632757202486`,4.135390983656999`,"calculate"},{"EARTH",5.97,12756,5515,9.8,11.2,23.9,24,149.6,147.098074,152.097701,365.2425,30.29,29.29,0,0.017,23.4,"calculate",6.094689747964199`,1.795100806273452`,5.731722873671973`,"calculate"},
{"MOON",0.073,3475,3340,1.6,2.4,655.7,708.7,0.384,0.363,0.406,27.3,1.076,0.964,5.1,0.055,6.7,"calculate",0,0,0,"calculate"},{"MARS",0.642,6792,3933,3.7,5,24.6,24.7,227.9,206.6,249.2,687,26.5,21.97,1.9,0.094,25.2,"calculate",0.8650968783605173`,5.865845667101206`,4.580229856896876`,"calculate"},
{"JUPITER",1899,142984,1326,23.1,59.5,9.9,9.9,778.6,740.5,816.6,4331,13.72,12.44,1.3,0.049,3.1,"calculate",1.7535549887589787`,0.2739782953195659`,5.6297310681731805`,"calculate"},{"SATURN",568,120536,687,9,35.5,10.7,10.7,1433.5,1352.6,1514.5,10747,10.18,9.09,2.5,0.057,26.7,"calculate",1.9833188581377723`,1.550951933199721`,0.3657789487830384`,"calculate"},{"URANUS",86.8,51118,1270,8.7,21.3,-17.2,17.2,2872.5,2741.3,3003.6,30589,7.11,6.49,0.8,0.046,97.8,"calculate",1.2932086905822064`,3.066206647208402`,5.291657999535076`,"calculate"},{"NEPTUNE",102,49528,1638,11,23.5,16.1,16.1,4495.1,4444.5,4545.7,59800,5.5,5.37,1.8,0.011,28.3,"calculate",2.3002130544346264`,0.1257684258987114`,5.233615853529528`,"calculate"},{"PLUTO",0.0125,2390*10,1750,0.6,1.1,-153.3,153.3,5870,4435,7304.3,90588,6.1,3.71,17.2,0.244,122.5,"calculate",1.9265520242166567`,3.923543791714553`,4.114885982939944`,"calculate"}};

PlanetaryData=PlanetaryData[[All,{1,11,10,15,8,3,17,12,18,14,13,19,20,21,22}]]; (*This line redefines the list to only include the needed data*)
PlanetaryData[[2;;11,{2,3}]]=(10^6)PlanetaryData[[2;;11,{2,3}]];
PlanetaryData[[2;;11,{4,7}]]=π/180*PlanetaryData[[2;;11,{4,11}]]//N;
PlanetaryData//TableForm

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=5; (*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>π}}];
PlanetaryData[[planet,15]]=PlanetaryData[[planet,8]]/TotalArea*Area[PlanetaryData[[planet,14]]];
(**********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: 
PlanetaryData[[Row,9]]/.t->value*)
(****************************************************)
For[i=2,i≤11,i++,PlanetaryData[[i,9]]=AnomalyModule[PlanetaryData[[i,2]],PlanetaryData[[i,3]],Floor[PlanetaryData[[i,8]]],i,0]]
PlanetaryData[[All,{9}]]//TableForm

(****************************************************)
(*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*)
(****************************************************)
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=2,i≤11,i++,
aphelion=PlanetaryData[[i,2]];
perihelion=PlanetaryData[[i,3]];
e=(aphelion-perihelion)/(aphelion+perihelion);(*the eccentricity of the orbit*)
r[θ_]=((A[aphelion,perihelion](1-e^2))/(1+e*Cos[θ]));
AppendTo[VelocityTable,{PlanetaryData[[i,1]],PlanetaryData[[i,10]],(1/60)NIntegrate[(Sqrt[(((r[θ]))^2)+(((r'[θ]))^2)]),{θ,PlanetaryData[[i,9]]/.t->((PlanetaryData[[i,8]])/2)-(1/(48*60)),PlanetaryData[[i,9]]/.t->((PlanetaryData[[i,8]])/2)+(1/(48*60))}],PlanetaryData[[i,11]],(1/(60^2))NIntegrate[(Sqrt[(((r[θ]))^2)+(((r'[θ]))^2)]),{θ,PlanetaryData[[i,9]]/.t->0-(1/48),PlanetaryData[[i,9]]/.t->0+(1/48)}],PlanetaryData[[i,15]]}];
];
Print[VelocityTable//TableForm];
]

The Earth's elliptical orbit in 3D (with tilt and everything else so far)

(****************************************************)
(*Build and graph our model of the Solary System. This
makes use of the anomoly models θ(t) found in the sub-
module*)
(****************************************************)
RangeSettings={{(10^8)*{{-3,3},{-3,3},{-0.5,0.5}},{2,3,4,6},1250, (10^7)*2},
{(10^8)*{{-18,18},{-18,18},{-2.5,2.5}},{7,8},1000, (10^7)*15},
{(10^8)*{{-50,70},{-70,70},{-20,20}},{9,10,11},9000,(10^7)*60},
{(10^8)*{{-50,70},{-70,70},{-20,20}},{7,8,9,10,11},3000, (10^7)*45}};

A[aph_,per_]=((aph+per)/2); (*Major axis of the ellipse*)
B[aph_,per_]=(1/2)(Sqrt[(((aph+per))^2)-(((aph-per))^2)]);(*Minor axis of the ellipse*)
Shift[aph_,per_]=((aph-per)/2); (*Shifts the orbits to be centered at the Sun (or planet in the case of a moon)*)
T[angle_]={{1,0,0},{0,Cos[angle],-Sin[ angle]},{0,Sin[angle],Cos[ angle]}}//N;

Tx[θ_]={{1,0,0},{0,Cos[θ],-Sin[θ]},{0,Sin[θ],Cos[θ]}}//N;
Ty[θ_]={{Cos[θ],0,Sin[ θ]},{0,1,0},{-Sin[ θ],0,Cos[θ]}}//N;
Tz[θ_]={{Cos[ θ],-Sin[ θ],0},{Sin[θ],Cos[θ],0},{0,0,1}}//N;(*3D Rotation matrix about the x-axis to tilt the planets*)

phelionOrbit[t_,aph_,per_]={A[aph,per]*Cos[t]-Shift[aph,per],B [aph,per]Sin[t],0};(*translation*)

Planet[θ_,ϕ_,r_]={r*Sin[ϕ]Cos[θ],r*Sin[ϕ]Sin[θ],r*Cos[ϕ]};
Orbit[θ_]=Table[(Tx[PlanetaryData[[i,4]]].Tz[PlanetaryData[[i,13]]]).phelionOrbit[θ, PlanetaryData[[i, 2]], PlanetaryData[[i, 3]]],{i,2,11}];

(****************************************************)
(*Manipulate to display the planets, their orbits,
rotation, moons, and some little green men if you look
closely enough*)
(*time is in earth days*)
(****************************************************)
Manipulate[
RangeOption=If[view=="Inner",1,If[view=="Giant",2,If[view=="Outer",3,4]]];
Module[{Plotrange,PlanetSelection,PlanetCount,Anomaly,Rotation,PositionXYZ,Planets},
Plotrange=RangeSettings[[RangeOption,1]];
PlanetSelection=RangeSettings[[RangeOption,2]];
PlanetCount=Length[PlanetSelection];
Anomaly=Table[PlanetaryData[[i,9]]/.t->(time+PlanetaryData[[i,15]]+876),{i,PlanetSelection}];
Rotation=Table[((24*2π*time)/(PlanetaryData[[i,5]])),{i,PlanetSelection}];
PositionXYZ=Table[Orbit[Anomaly[[i]]][[PlanetSelection[[i]]-1]],{i,1,PlanetCount}];
Planets[θ_,ϕ_]=Table[PositionXYZ[[i]]+Tx[PlanetaryData[[PlanetSelection[[i]],7]]].Planet[Rotation[[i]]+θ,ϕ,PlanetaryData[[PlanetSelection[[i]],6]]]*RangeSettings[[RangeOption,3]],{i,1,PlanetCount}];
  
(*graph the elliptical orbit*)
Show[ ParametricPlot3D[Orbit[θ], {θ, 0, 2*π}, PlotRange ->Plotrange, PlotStyle -> Directive[Dashing[0.02],White],Ticks->False,Background->Black,Boxed->False],

  (*plot of the planets*)
  ParametricPlot3D[Planets[ θ, ϕ],{θ,0,1.9π}, {ϕ, 0, π},
(*options for display*)
PlotRange -> Plotrange,Mesh->False,PlotPoints->8
  ,Ticks->True],
Graphics3D[{Glow[Yellow],Sphere[{0,0,0},RangeSettings[[RangeOption,4]]]}]
]]


(*Manipulate Controllers*)
, Style["Model of the Solar System",Bold,12],(*Title*)
{{time,4080,"days"}, If[view=="Inner",4000,0], If[view=="Inner",5000,If[view=="Giant",30000,250000]],1,ImageSize->Small,Animator,DefaultDuration->Dynamic[speed],AnimationRunning->False},{{speed,50,""},{25->"Fast",50->"Medium",100->"Slow"}},(*Time slider*)

Delimiter,
{{view,"Inner","Planets"},{"Inner","Giant","Outer","All"}}(*Buttons to select the window*), 

(*Notes*)
Style["",Plain,8],
Style["Planet size scale",Bold,8],
Style["Mercury: 1250",Plain,8],
Style["Venus: 1250",Plain,8],
Style["Earth: 1250",Plain,8],
Style["Mars: 1250",Plain,8],
Style["Jupiter: 1250",Plain,8],
Style["Saturn: 1250",Plain,8],
Style["Uranus: 9,000",Plain,8],
Style["Neptune: 9,000",Plain,8],
Style["Pluto: 90,000",Plain,8],
ControlPlacement->Left]

Manipulate[x,{x,1,60,Animator,DefaultDuration->Dynamic[speed],AnimationRunning->False},{{speed,60,""},{2->"Fast",60->"Slow"},RadioButtonBar}]

Appearance->Labeled,