(*  This is a simulation of an Earth to Mars trajectory  *)
(*  It illustrates how one might use a Runge-Kutta method to integrate the equations of motion *)

(*  These are the coefficients of a fourth-order RK method  *)
(*  You are not to use these coefficients for your method   *) 
a[1]=1/2;a[2]=1/2;a[3]=1;b[1,0]=1/2;b[2,0]=1/4;b[2,1]=1/4;b[3,0]=0;b[3,1]=-1;b[3,2]=2;c[0]=1/6;c[1]=0;c[2]=2/3;c[3]=1/6;

(* This is the right-hand-side of the system of differential equations *)

f[t_,x_]:={x[[2]],-x[[1]]/(x[[1]]^2+x[[3]]^2)^1.5,x[[4]],-x[[3]]/(x[[1]]^2+x[[3]]^2)^1.5,x[[6]],-x[[5]]/(x[[5]]^2+x[[7]]^2)^1.5,x[[8]],-x[[7]]/(x[[5]]^2+x[[7]]^2)^1.5,x[[10]],-x[[9]]/(x[[9]]^2+x[[11]]^2)^1.5,x[[12]],-x[[11]]/(x[[9]]^2+x[[11]]^2)^1.5};

(*  These are parameters that could be changed if you like  *)

fps=1;                                   (*  frames/second for the animation  *)
sizeplot=1.7;                  (*  size of the plot  *)
imagesz=0.75 72*8;   (*  size of the image  *)
sma=1.52369;                   (*  semi-major axis of Mars  *)
smasv=0.5*(sma+1.0);   (*  semi-major axis of space vehicle  *)pmars =2.0*Pi*sma^-1.5;psv=2.0*Pi*(smasv)^-1.5;nsteps=53;                    (*  number of steps for the simulation  *)ssv=Sqrt[2.0-1.0/smasv]-0.005;
smars=Sqrt[1.0/sma];
x0={0.0,1.0,1.0,0.0,
sma/Sqrt[2],0.805/Sqrt[2],sma/Sqrt[2],-0.805/Sqrt[2],
0.0,ssv,1.0,0.0};    (*  initial conditions for all three bodies  *)ptlist={{x0[[1]],x0[[3]]}};  (*  these lists keep track of their positions  *)ptlist2={{x0[[5]],x0[[7]]}};ptlist3={{x0[[9]],x0[[11]]}};t0=0.0;     
tf=2*Pi;      (* t0 is the initial time and  tf is the final time  *)
h=(tf-t0)/nsteps;  (* the stepsize is h  *)
For[i=0,i<nsteps,i++,     (*  for each step of size h, we use rk4  *)
f0=f[t0,x0];                              (*  to advance the positions and velocities  *)
x[1]=x0+h*b[1,0]*f0;
f1=f[t0+a[1]*h,x[1]];
x[2]=x0+h*(b[2,0]*f0+b[2,1]*f1);
f2=f[t0+a[2]*h,x[2]];
x[3]=x0+h*(b[3,0]*f0+b[3,1]*f1+b[3,2]*f2);
f3=f[t0+a[3]*h,x[3]];
x[4]=x0+h*(c[0]*f0+c[1]*f1+c[2]*f2+c[3]*f3);
x0=x[4]; t0=t0+h;
ptlist=Append[ptlist,{x0[[1]],x0[[3]]}];
ptlist2=Append[ptlist2,{x0[[5]],x0[[7]]}];
ptlist3=Append[ptlist3,{x0[[9]],x0[[11]]}];
lps[i]=ListPlot[{ptlist,ptlist2,ptlist3},PlotRange->{{-sizeplot,sizeplot},{-sizeplot,sizeplot}},ImageSize->imagesz,PlotStyle->{{Darker[Cyan],PointSize[0.008]},{Red,PointSize[0.008]},{Lighter[Purple],PointSize[0.007]}},AspectRatio->Automatic,Background->RGBColor[0.0,0.2,0],PlotLabel->Style[Framed["Earth-to-Mars Trajectory\nTime = "<>ToString[i]<>" weeks"],12,White,Background->Darker[Gray]]];
]

lp=ListPlot[ptlist,Joined->True,AspectRatio->Automatic,PlotStyle->Blue,PlotRange->{{-sizeplot,sizeplot},{-sizeplot,sizeplot}},Background->RGBColor[0,0.5,0,.5]];lp2=ListPlot[ptlist2,Joined->True,AspectRatio->Automatic,PlotStyle->Red,PlotRange->{{-sizeplot,sizeplot},{-sizeplot,sizeplot}},Background->RGBColor[0,0.5,0,.5]];lp3=ListPlot[ptlist3,Joined->True,AspectRatio->Automatic,PlotStyle->Purple,PlotRange->{{-sizeplot,sizeplot},{-sizeplot,sizeplot}},Background->RGBColor[0,0.1,0]];


Show[lp,lp2,lp3]    (*  if you want to trace the paths that the bodies take  *)

ListAnimate[Table[lps[i],{i,1,nsteps-1}],fps]   (*  the animation of the spacecraft, Earth, and Mars over one year   *)