Back to Kepler

Kepler's 3rd law

Data below is a sample set of positional data of Explorer 35

in its eleiptical orbit around the Moon

(Reference: Activities in Astronomy by D.Hoff, L.Kelsey and J.Neff)

The following program is not a compact code. I purposely coded loosly for the purpose of learnig Mathematica in steps. Please email me if you are interested to know more about how to use Mathematica in education or in research.

Click for more Mathematica programming Problems/solutions with Mathematica

H.Tahsiri

Off[General::spell];
Off[General::spell1];
Clear["Global`*"];

(*  ydata as a function of t, given in the above reference,is enterd as {t,y} *)

ydata={{0,1.04},{.25,.63},{.50,.20},{.75,-.22},
{1,-.65},{1.25,-1.03},{1.50,-1.37},{1.75,-1.58},
{2.00,-1.59},{2.25,-1.32},{2.50,-.79},{2.75,-.11},
{3.00,.58},{3.25,1.22},{3.50,1.82},{3.75,2.35},
{4.00,2.81},{4.25,3.22},{4.50,3.59},{4.75,3.90},
{5.00,4.16},{5.25,4.40},{5.50,4.58},{5.75,4.74},
{6.00,4.86},{6.25,4.95},{6.50,5.01},{6.75,5.03},
{7.00,5.04},{7.25,5.00},{7.50,4.95},{7.75,4.87},
{8.00,4.77},{8.25,4.65},{8.50,4.50},{8.75,4.33},
{9.00,4.14},{9.25,3.93},{9.50,3.69},{9.75,3.42},
{10.00,3.15},{10.25,2.85},{10.50,2.52},{10.75,2.20},
{11.00,1.83},{11.25,1.46},{11.50,1.06},{11.75,.65}}

     {{0, 1.04}, {0.25, 0.63}, {0.5, 0.2}, {0.75, -0.22}, {1, -0.65}, {1.25, -1.03}, {1.5, -1.37}, 
      
       {1.75, -1.58}, {2., -1.59}, {2.25, -1.32}, {2.5, -0.79}, {2.75, -0.11}, {3., 0.58}, 
      
       {3.25, 1.22}, {3.5, 1.82}, {3.75, 2.35}, {4., 2.81}, {4.25, 3.22}, {4.5, 3.59}, {4.75, 3.9}, 
      
       {5., 4.16}, {5.25, 4.4}, {5.5, 4.58}, {5.75, 4.74}, {6., 4.86}, {6.25, 4.95}, {6.5, 5.01}, 
      
       {6.75, 5.03}, {7., 5.04}, {7.25, 5.}, {7.5, 4.95}, {7.75, 4.87}, {8., 4.77}, {8.25, 4.65}, 
      
       {8.5, 4.5}, {8.75, 4.33}, {9., 4.14}, {9.25, 3.93}, {9.5, 3.69}, {9.75, 3.42}, {10., 3.15}, 
      
       {10.25, 2.85}, {10.5, 2.52}, {10.75, 2.2}, {11., 1.83}, {11.25, 1.46}, {11.5, 1.06}, 
      
       {11.75, 0.65}}

(* xdata as a function of t, given in the Manual, is enterd as {t,x} *)

xdata={{0,-3.62},{.25,-3.46},{.50,-3.25},{.75,-2.97},
{1,-2.60},{1.25,-2.14},{1.50,-1.55},{1.75,-.85},
{2.00,-.03},{2.25,.78},{2.50,1.45},{2.75,1.87},
{3.00,2.09},{3.25,2.16},{3.50,2.11},{3.75,1.99},
{4.00,1.82},{4.25,1.61},{4.50,1.37},{4.75,1.11},
{5.00,.85},{5.25,.58},{5.50,.28},{5.75,0},
{6.00,-.27},{6.25,-.56},{6.50,-.84},{6.75,-1.12},
{7.00,-1.38},{7.25,-1.64},{7.50,-1.89},{7.75,-2.14},
{8.00,-2.37},{8.25,-2.59},{8.50,-2.80},{8.75,-2.99},
{9.00,-3.17},{9.25,-3.33},{9.50,-3.49},{9.75,-3.59},
{10.00,-3.69},{10.25,-3.77},{10.50,-3.81},{10.75,-3.83},
{11.00,-3.81},{11.25,-3.75},{11.50,-3.65},{11.75,-3.51}}

     {{0, -3.62}, {0.25, -3.46}, {0.5, -3.25}, {0.75, -2.97}, {1, -2.6}, {1.25, -2.14}, 
      
       {1.5, -1.55}, {1.75, -0.85}, {2., -0.03}, {2.25, 0.78}, {2.5, 1.45}, {2.75, 1.87}, 
      
       {3., 2.09}, {3.25, 2.16}, {3.5, 2.11}, {3.75, 1.99}, {4., 1.82}, {4.25, 1.61}, {4.5, 1.37}, 
      
       {4.75, 1.11}, {5., 0.85}, {5.25, 0.58}, {5.5, 0.28}, {5.75, 0}, {6., -0.27}, {6.25, -0.56}, 
      
       {6.5, -0.84}, {6.75, -1.12}, {7., -1.38}, {7.25, -1.64}, {7.5, -1.89}, {7.75, -2.14}, 
      
       {8., -2.37}, {8.25, -2.59}, {8.5, -2.8}, {8.75, -2.99}, {9., -3.17}, {9.25, -3.33}, 
      
       {9.5, -3.49}, {9.75, -3.59}, {10., -3.69}, {10.25, -3.77}, {10.5, -3.81}, {10.75, -3.83}, 
      
       {11., -3.81}, {11.25, -3.75}, {11.5, -3.65}, {11.75, -3.51}}

y[t_]=Interpolation[ydata][t]

     InterpolatingFunction[{{0, 11.75}}, <>][t]

x[t_]=Interpolation[xdata][t]

     InterpolatingFunction[{{0, 11.75}}, <>][t]

p1=Plot[Evaluate[y[t]],{t,0,11.56},
		Epilog->Map[{PointSize[0.02],Point[#]}&,ydata]];

Clear[vy]

vy[t_]=y'[t]

     InterpolatingFunction[{{0, 11.75}}, <>][t]

vytable=Table[{t,vy[t]},{t,0,11.56,.2}];

p11=Plot[Evaluate[vy[t]],{t,0,11.56}];

p2=Plot[Evaluate[x[t]],{t,0,11.56},
		Epilog->Map[{PointSize[0.02],Point[#]}&,xdata]];

vx[t_]=x'[t]

     InterpolatingFunction[{{0, 11.75}}, <>][t]

p22=Plot[Evaluate[vx[t]],{t,0,11.56},PlotRange->All];

vxtable=Table[{t,vx[t]},{t,0,11.56,.2}];

Clear[v]

v[t_]=Sqrt[vx[t]^2+vy[t]^2]

                                                    2                                             2
     Sqrt[InterpolatingFunction[{{0, 11.75}}, <>][t]  + InterpolatingFunction[{{0, 11.75}}, <>][t] ]

vtable=Table[{t,v[t]},{t,0,11.56,.2}];

p3=ParametricPlot[Evaluate[{x[t],y[t]}],{t,0,11.56},AspectRatio->Automatic];

s1=Map[Take[#,{2,2}]&,xdata];

s2=Map[Take[#,{2,2}]&,ydata];

s3=Flatten[Thread[{s1,s2}]];

s4=Partition[s3,2];

(* r=Sqrt[x^2+y^2] *)

r1=Sqrt[s4[[1,1]]^2+s4[[1,2]]^2]
     3.76643

r2=Sqrt[s4[[2,1]]^2+s4[[2,2]]^2]
     3.51689

Clear[radi]

radi=Table[{Sqrt[s4[[i,1]]^2+s4[[i,2]]^2]},{i,1,Length[s4]}];

(* As shown bellow, rmin=1.53323 and rmax=5.32633*)

Sort[radi];

(* semimajor Axis=a *)

a=((rmin+rmax)/2)/.{rmin->1.53323,rmax->5.32633}

     3.42978

(*ecen=e *)

Clear[e]

e=(rmax-rmin)/(2a)/.{rmax->5.3263,rmin->1.53323,a->3.42978}

0.552961

(* semiminor=b *)

b=a Sqrt[(1-e^2)/.{a->3.42978,e->.552961}]

2.85772

Clear[a,e,b,r]

r=Sqrt[(x-a e)^2+y^2]

(* This is Vis Viva orbit equations *)

v1=Sqrt[G M ((2/r)-(1/a))]

given={G->8.642 10^(-13),M->7.35 10^22,a->3.42978*1738,e->.552961}

v2=v1/.given//N

(* x and y values at t=1 hr. See lab table *)

v3=v2/1000/.{x->-2.6*1738,y->-.65*1738}

2.3308

(* the abov theoretical value of speed is in an
	excellent agreement with the experiment *)

s5=ListPlot[s4,PlotJoined->True,PlotRange->All,AspectRatio->Automatic];

s6=ListPlot[s4,PlotJoined->True,AspectRatio->Automatic,	
				Epilog->Map[{PointSize[0.02],Point[#]}&,s4]];

cont=	 Epilog->{{Text["x(t) = ",{1.7,5.5}]},{Text[x[i],{3.2,5.5}]},{Text["y(t) = ",{1.7,5}]},{Text[y[i],{3.2,5}]},
			{Text["t = ",{2,4.5}]},{Text[i,{3,4.5}]},
			{Text["v(t) = ",{3.1,4}]},{Text[v[i],{4.6,4}]},
			Text["Moon and Satellite",{-3,6}],
			Text["e = .553",{3.9,3.5}],
			Text["a = 3.43",{3.9,3}],
			Text["b = 2.86",{3.9,2.5}],
			Text["rmax = 5.33",{3.9,2}],
			Text["rmin = 1.53",{3.9,1.5}]}

s6=Do[Show[p3,
			Graphics[{GrayLevel[.5],PointSize[.19],Point[{0,0}],RGBColor[1,0,0],PointSize[0.04],Point[{x[i],y[i]}] } ], cont,PlotRange->{{-5,5},{-2,6}},AxesLabel->{x,  y},Background->RGBColor[0,0,0],Frame->False],
		{i,0,0}];

Back to Kepler