Kepler’s First Law

Andy Long
Spring, 2025

Based on material from section 13.4 of Stewart

Definitions:

I chose the following values so as to produce nice plots which illustrate the motion of the planet: they don’t reflect the actual values for our planet!

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 scales the eccentricity -- that is, how much the ellipse varies from a circle (in which case c=0).

h=1;
c=10;
eccentricity=.6;(* Equals c/(GM) in reality *)(* So eccentricity is dependent 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[θ]
p[θ_]:={x[θ],y[θ]}; (* position vector *)
pPerp[θ_]:={-y[θ],x[θ]}; (* perpendicular to position *)
v[θ_]:=p'[θ];(* velocity vector *)
a[θ_]:=v'[θ];(* acceleration vector *)

Plotting the orbit, and initial motion:

p1=ParametricPlot[p[θ],{θ,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[{p[θ]},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 p[th0],{t,0,1}];
p5=ParametricPlot[
t((r'[th0]/r[th0])p[th0]+pPerp[th0])/r[th0]^2
+
p[th0],
{t,0,1},PlotStyle->Orange];
Show[p3,p1,p4,p5]

Fixing the speed

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

We start from
					 (UnderBar[h])=(r^2)((OverHat[u]) x (d/(dt))((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)| (d/(dt))((OverHat[u])) | .

The magnitude of the vector change (d/(dt))((OverHat[u])) will be ((dθ)/dt), because (OverHat[u]) is a unit vector -- only its angle will change with time, not its magnitude (or length). 

Therefore   |(d/(dt))((OverHat[u])) | = ((dθ)/dt).
 
Finally we have that 
					 ((dθ)/dt)=(h/(r^2))=(h/((r[θ])^2))
					 
(which I couldn’t solve analytically). It’s separable, since  
 
 					((r[θ])^2)dθ =  hdt 
 
 So we integrate both sides with respect to time, to get
 
 				∫((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 didn’t give anything nice and illuminating. 
 
Integrate[(eccentricity d/(1+eccentricity Cos[th]))^2, {th,0,θ}]

So let’s 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) : we go back to the original equation,

					 ((dθ)/dt)=(h/(r^2))=(h/((r[θ])^2))
					 
and integrate both sides with respect to time, to find θ(t):

end=.0442 (* this value works for the parameter values I chose, to produce one period of the motion. *)
s=NDSolve[{θ'[t]==h/r[θ[t]]^2,θ[0]==0},θ,{t,0,end}]

What we obtain is a really good numerical solution for θ(t). Note that it takes about end = .0442 units of time for θ(t) to go through 2 π :


We see from this graph that the time rate of change of θ(t) is greatest (has greatest slope) around zero, and the least time rate of change around .022 (when the value of θ(t) is π -- i.e., when we’re farthest from the sun).

Plot[Evaluate[θ[t]/.s],{t,0,end},PlotRange->All]
{θ[0]/.s,θ[end]/.s,2.0Pi}

Plotting the orbit with the correct speed

Now let’s 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[{p[point]},PlotMarkers->{}]
]
],{t,0,end,.00001}
]

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

Manipulate[
Module[
{rvector,
vvector,
avector, 
secondVvector, 
angle,
th0=θ[t]/.s,
dthdt,
normv,
normr,
p4, p5, p6,p7},
th0=th0[[1]];
rvector=p[th0];
normr=Norm[rvector];
dthdt=h/normr^2;
vvector=v[th0] dthdt; (* v = dr/dt = dr/dθ dθ/dt *)
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=
a[th0]dthdt^2
+
vvector (-2  Sin[th0]/d)/normr;

angle=ArcCos[Dot[vvector,rvector]/normv/normr];

(* Now we'll plot these velocity and acceleration vectors from the position of the point in motion: *)
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]p[th0]/r[th0]^3+pPerp[th0]/r[th0]^2;
p7=ParametricPlot[t secondVvector+p[th0],{t,0,.005},PlotStyle->{Orange,Thick}];
Show[p5,p1,p2,p4,p6,
ListPlot[p[th0],PlotMarkers->{}]
]
],
{t,0,end,.000001}
]