Kepler’s First Law

Andy Long
Spring, 2019

This is from Chapter 13 of Stewart (end of section 13.4).

I chose the following at (somewhat) random: you want a little “eccentricity” to your ellipse (not too much like a circle). h is related to the angular momentum. c is the length of the constant vector we obtain in the derivation, and determines the eccentricity -- that is, how much the ellipse varies from a circle (c=0).

h=1;
c=10;
eccentricity=.6; (* Equals c/(GM) in reality -- so dependent as well on universal gravitation and the mass of the sun. *)
d=h^2/c;
r[θ_]:=eccentricity d/(1+eccentricity Cos[θ])
x[θ_]:=r[θ]Cos[θ]
y[θ_]:=r[θ]Sin[θ]
p1=ParametricPlot[{x[θ],y[θ]},{θ,0,2π}];
p2 = ListPlot[{{0,0}},PlotMarkers->{""}];
Show[p1,p2]

If you plot r(θ) as θ changes uniformly, you’ll see that the “planet” is moving at the wrong speeds: it should be faster when it’s near the sun, and slower when it’s far away. That’s because the angular change in the plot is constant, which is not the case for planetary motion:

Manipulate[
Show[
p1,
p2,
ListPlot[{{x[θ],y[θ]}},PlotMarkers->{""}]
],{θ,0,2Pi}]

The velocity vector runs along the curve traced out by r:

th0=2Pi/3;
dydx[t_]:=y'[t]/x'[t]
tangentLine[t_]:=y[th0]+dydx[th0](x[t]-x[th0])
p3=ParametricPlot[{x[t],tangentLine[t]},{t,th0-Pi/12,th0+Pi/12},PlotLabel->StringJoin["θ =",ToString[th0,InputForm]],
PlotRange->{{-.2,.1},{-.1,.1}}];
p4=ParametricPlot[t{x[th0],y[th0]},{t,0,1}];
p5=ParametricPlot[
t(
r'[th0]/r[th0]{x[th0],y[th0]}+{-y[th0],x[th0]}
)/r[th0]^2
+
{x[th0],y[th0]},
{t,0,1},PlotStyle->Orange];
Show[p3,p1,p4,p5]

To get the speed right, we would need to solve a “differential equation” to find θ as a function of time t -- that is, θ(t). (Make sure that you take DEs one day!):

From (UnderBar[h])=(r^2)((OverHat[u]) x (OverHat[u])’):

Since the length of (OverHat[u]) can’t change, the only change possible is in direction -- perpendicular to the direction of (OverHat[u]) (i.e., increasing or decreasing θ, with no change in the length “r” -- which is one in this case). Hence, computing the norm of (UnderBar[h]), we get

h=(r^2)| (OverHat[u])’ | .

The change  (OverHat[u])’ will be ((r dθ)/dt), but the length r is 1 (this is a unit vector). Therefore
 | (OverHat[u])’ | = ((dθ)/dt)
Finally,
 ((dθ)/dt)=(h/(r^2))=(h/((r[θ])^2))
 (which I couldn’t solve analytically). It’s separable, so  ∫((r[θ])^2)dθ =  h∫dt = h t + C.
 
 By the way, this proves Kepler’s equal areas  formula, too: if we turn these into definite integrals, then it says for a given time span T,
 
 (∫_0^T)((r[θ(t)])^2)dθ(t) =  h(∫_0^T)dt = h T 
 
 which says the area traced out is ((h T)/2) -- equal areas are traced out in equal times. 
 
 I tried using Mathematica to compute this integral, but it just spun:
 
Integrate[(eccentricity d/(1+eccentricity Cos[th]))^2, {th,0,θ}]

So we’ll do it numerically (approximately): the end time was chosen to give us values of θ from 0 to 2π (a complete period of the planet’s orbit) :

end=.0442 (* this works for the parameter values I chose! *)
s=NDSolve[{θ'[t]==h/r[θ[t]]^2,θ[0]==0},θ,{t,0,end}]
Plot[Evaluate[θ[t]/.s],{t,0,end},PlotRange->All]
{θ[0]/.s,θ[end]/.s,2.0Pi}

What we obtain is a really good numerical solution for θ(t). 

Now we want to plot the orbit as a function of time, including the sun’s location (indicated by the dodecahedron), and the planet’s position (indicated by the icosahedron). These markers are chosen in honor of Kepler, i.e. 
https://upload.wikimedia.org/wikipedia/commons/1/19/Kepler-solar-system-1.png

Manipulate[
Module[{point=θ[t]/.s},
point=point[[1]];
Show[
p1,
p2,
ListPlot[{{r[point]Cos[point],r[point]Sin[point]}},PlotMarkers->{}]
]
],{t,0,0.044}
]

Now I’ll add in a position vector, velocity, and acceleration, and we’ll have it all!

Manipulate[
Module[{rvector,vvector,secondVvector, avector, 
angle,th0=θ[t]/.s,dthdt,
norma, normv,normr,normsecond,
p4, p5, p6,p7},
th0=th0[[1]];
rvector={x[th0],y[th0]};
normr=Norm[rvector];
dthdt=h/normr^2;
(* v = dr/dt = dr/dθ dθ/dt *)
vvector={x'[th0],y'[th0]} dthdt;
normv=Norm[vvector];
(* 
a = (d^2)r/(dt^2) = (d^2)r/((dθ)^2) (((dθ/dt))^2) + dr/dθ (d^2)θ/(dt^2) 
D[gg[ff[t]],{t,2}]
So
(d^2)θ/(dt^2) = (1/r^2)' = -2(1/r^3) dr/dt
        = -2(1/r^3) r^2Sin[θ]/d * dθ/dt
       = -2/r Sin[θ]/d * dθ/dt
And
dr/dθ (d^2)θ/(dt^2) = dr/dθ -2/r Sin[θ]/d * dθ/dt
				= vvector -2Sin[θ]/d /r 
*)
avector=
{x''[th0],y''[th0]}dthdt^2
+
vvector (-2  Sin[th0]/d)/normr;
norma=Norm[avector];
angle=ArcCos[Dot[vvector,rvector]/normv/normr];
p5=ParametricPlot[t vvector+rvector,{t,0,.005},PlotStyle->{Red,Thick},
PlotRange->{{-.25,.1},{-.1,.1}},PlotLabel->angle];
p6=ParametricPlot[t avector+rvector,{t,0,.00003},PlotStyle->{Black,Thick}];
p4=ParametricPlot[t rvector,{t,0,1}];secondVvector=r'[th0]{x[th0],y[th0]}/r[th0]^3+{-y[th0],x[th0]}/r[th0]^2;
normsecond=Norm[secondVvector];
p7=ParametricPlot[t secondVvector+{x[th0],y[th0]},{t,0,.005},PlotStyle->{Orange,Thick}];
Show[p5,p1,p2,p4,p6,
ListPlot[{{r[th0]Cos[th0],r[th0]Sin[th0]}},PlotMarkers->{}]
]
],
{t,0,end}
]