Chapter 4 
The Solar System

Sample code 4.1.1. Orbit of Earth around the Sun.

Deltat=0.002 (*in unit of yr *);

(* Planet-spesific data for Earth*)
a=1; (* in A.U. *)
vx0=0.0;(* initial x-component of velocity, in unit of A.U./yr *)
vy0=2Pi/Sqrt[a];(* initial y-component of velocity, in unit of A.U./yr *)
tfinal =a^(3/2);(* in unit if yr *)
(* End of Planet-spesific data *)

NStep =IntegerPart[ (tfinal/Deltat)+1]; (* Number of total time step *)
x0=a;(* initial x-coordinate, in unit of A.U. *)
y0=0.0;(* initial x-coordinate, in unit of A.U. *)
x[0]=x0;t[0]=0;y[0]=y0;
vx[0]=vx0;vy[0]=vy0;

For [i=0,i≤NStep,i++,
r[i]=Sqrt[x[i]^2+y[i]^2];
vx[i+1]=vx[i]-4*(Pi^2)*x[i]*Deltat/r[i]^3;
vy[i+1]=vy[i]-4*(Pi^2)*y[i]*Deltat/r[i]^3;
x[i+1]=x[i]+Deltat*vx[i+1];
y[i+1]=y[i]+Deltat*vy[i+1];
t[i+1]=t[i]+Deltat;
];(* end for *)
txy=Table[{x[i],y[i]},{i,0,NStep}];
ListPlot[txy,PlotRange->{-1,1},Joined->True,AspectRatio->1,AxesLabel->{"x(A.U)","y(A.U)"}]

Sample code 4.1.2. Overlap orbit of Earth around the Sun with that of Mars.

Deltat=0.002 ;(*in unit of yr *)

(* Planet-spesific data for Earth *)
a=1; (* in A.U. *)
vx0=0.0;(* initial x-component of velocity, in unit of A.U./yr *)
vy0=2Pi/Sqrt[a];(* initial y-component of velocity, in unit of A.U./yr *)
tfinal = 0.98a^(3/2);
(* End of Planet-spesific data *)

NStep =IntegerPart[ (tfinal/Deltat)+1]; (* Number of total time step *)
x0=a;(* initial x-coordinate, in unit of A.U. *)
y0=0.0;(* initial x-coordinate, in unit of A.U. *)
x[0]=x0;t[0]=0;y[0]=y0;
vx[0]=vx0;vy[0]=vy0;

For [i=0,i≤NStep,i++,
r[i]=Sqrt[x[i]^2+y[i]^2];
vx[i+1]=vx[i]-4*(Pi^2)*x[i]*Deltat/r[i]^3;
vy[i+1]=vy[i]-4*(Pi^2)*y[i]*Deltat/r[i]^3;
x[i+1]=x[i]+Deltat*vx[i+1];
y[i+1]=y[i]+Deltat*vy[i+1];
t[i+1]=t[i]+Deltat;
];(* end for *)
txyE=Table[{x[i],y[i]},{i,0,NStep}];



(* Planet-spesific data for Mars *)
a=1.52; (* in A.U. *)
vx0=0.0;(* initial x-component of velocity, in unit of A.U./yr *)
vy0=2Pi/Sqrt[a];(* initial y-component of velocity, in unit of A.U./yr *)
tfinal = 0.98a^(3/2);(* End of Planet-spesific data *)


NStep =IntegerPart[ (tfinal/Deltat)+1]; (* Number of total time step *)
x0=a;(* initial x-coordinate, in unit of A.U. *)
y0=0.0;(* initial x-coordinate, in unit of A.U. *)
x[0]=x0;t[0]=0;y[0]=y0;
vx[0]=vx0;vy[0]=vy0;
For [i=0,i≤NStep,i++,
r[i]=Sqrt[x[i]^2+y[i]^2];
vx[i+1]=vx[i]-4*(Pi^2)*x[i]*Deltat/r[i]^3;
vy[i+1]=vy[i]-4*(Pi^2)*y[i]*Deltat/r[i]^3;
x[i+1]=x[i]+Deltat*vx[i+1];
y[i+1]=y[i]+Deltat*vy[i+1];
t[i+1]=t[i]+Deltat;
];(* end for *)
txyMrs=Table[{x[i],y[i]},{i,0,NStep}];

ListPlot[{txyE,txyMrs},PlotRange->{{-2,2},{-2,2}},Joined->True,AspectRatio->1,AxesLabel->{"x(A.U)","y(A.U)"},PlotStyle->{Red,Blue}]

Sample code 4.1.3. Perturbing the initial velocity of Mars to render a non-zero eccentricity to the otherwise circular orbit, e.g., setting vy0 = 0.5 of 2Pi/Sqrt[a]. The period gets modified as well.

Deltat=0.002 ;(*in unit of yr *)

(* Planet-spesific data for Mars *)
a=1.52; (* in A.U. *)
vx0=0.2*2Pi/Sqrt[a];(* initial x-component of velocity, in unit of A.U./yr *)
vy0=0.5*2Pi/Sqrt[a];(* initial y-component of velocity, in unit of A.U./yr *)
tfinal = 0.43a^(3/2);(* End of Planet-spesific data *)

NStep =IntegerPart[ (tfinal/Deltat)+1]; (* Number of total time step *)
x0=a;(* initial x-coordinate, in unit of A.U. *)
y0=0.0;(* initial x-coordinate, in unit of A.U. *)
x[0]=x0;t[0]=0;y[0]=y0;
vx[0]=vx0;vy[0]=vy0;
For [i=0,i≤NStep,i++,
r[i]=Sqrt[x[i]^2+y[i]^2];
vx[i+1]=vx[i]-4*(Pi^2)*x[i]*Deltat/r[i]^3;
vy[i+1]=vy[i]-4*(Pi^2)*y[i]*Deltat/r[i]^3;
x[i+1]=x[i]+Deltat*vx[i+1];
y[i+1]=y[i]+Deltat*vy[i+1];
t[i+1]=t[i]+Deltat;
];(* end for *)
txyMrs=Table[{x[i],y[i]},{i,0,NStep}];

ListPlot[{txyMrs},PlotRange->All,Joined->True,AspectRatio->1,AxesLabel->{"x(A.U)","y(A.U)"},PlotStyle->{Red,Blue}]

Hint: How would you use Mathematica to “measure” the sizes of the semimajor, semiminor, hence the position of the perihelion (furtherst point from the Sun) and apehelion (closest distance to Sun)? How do you determine the planet’s period?

Sample code: 4.1.4. Excercise to abstract info from List.

Print[
Table[txyMrs[[i]],{i,1,20}]
];(*end print*)

Print[txyMrs[[1]]," ",txyMrs[[1,1]]," ",txyMrs[[1,2]]];
Print[txyMrs[[5]]," ",txyMrs[[5,1]]," ",txyMrs[[5,2]]];

Sample code: 4.1.5. Code to scan the distances of all pairs of points (i, j), and select that max of those.

(*Define the functions*)
xdum[idum_]:=txyMrs[[idum,1]];ydum[idum_]:=txyMrs[[idum,2]];
r[i_,j_]:=Sqrt[(xdum[i]-xdum[j])^2+(ydum[i]-ydum[j])^2];
(* The := sign means the function will not be defined until they are called *)

rmax=0;
For[j=1,j≤NStep,j++,
For[i=1,i≤NStep,i++,
If[r[i,j]>rmax,
imax=i;jmax=j;rmax=r[i,j]
];(*end if*)
];(*end for loop i*)
];(* end for loop j*)
Print["imax , jmax , rmax =",imax," ",jmax," ",rmax];
Print["{x[imax],y[imax]},{x[jmax],y[jmax]} =", {x[imax],y[imax]}," ",{x[jmax],y[jmax]}];

Sample code: 4.1.6: Display the part of locus between the two end points  {x[imax],y[imax]}, {x[jmax],y[jmax]}, using the command Graphics. These points are the perihelion and ephelion.

tpoints=Table[{x[idum],y[idum]},{idum,jmax,imax}];
Graphics[Point[tpoints]]