Physics 2300 Mathematica session #8

Here are some visualization tools to help with the homework.  The Math08 homework task is to visit the sections regarding Morin 7.2 and 7.4 below, and use them to make plots of three orbits.

Morin 5.86 (Drag Force on a sphere)

A round object (sphere or cylinder) of radius R moves at velocity V through a medium with a
certain number density ([n]=#/L^3) of tiny particles each of mass m.
Equivalently, the round object is at rest while a “wind” of particles streams by.
Let’s make a 2D movie for the case of the cylinder.

Begin["Morin522`"]
num=100;
Do[(
x0[j]=Random[Real,{0,10}]; (*For each particle, *)
y0[j]=Random[Real,{-1,1}];(* pick a random starting place*)
θ0[j] = ArcSin[y0[j]];   (* Determine where on the disk this one hits *)
t0[j] = x0[j]-Cos[θ0[j]]; (*And when it will hit*)
pt[j,t_] = If[t<t0[j],{x0[j]-t,y0[j]},{Cos[θ0[j]],Sin[θ0[j]]}+(t-t0[j]){Cos[2θ0[j]],Sin[2θ0[j]]}];),
{j,1,num}] (*and put together the two pieces of the motion *)
Animate[Show[Graphics[{Blue,Disk[]}],ListPlot[Table[pt[j,t],{j,1,num}],PlotStyle->Red],PlotRange->{{-1,5},{-1.5,1.5}}],{t,0,5}]
End[];

Next let’s move in to three dimensions.

Begin["Morin586`"]
sphere=ParametricPlot3D[{Cos[ϕ],Sin[ϕ]Cos[θ],Sin[θ]Sin[ϕ]},{ϕ,0,π},{θ,0,2π},PlotStyle->Opacity[0.8]];
ϵ=.1;
Manipulate[Show[sphere,ParametricPlot3D[{Cos[α],Sin[α]Cos[β],Sin[α]Sin[β]},{β,0,2π},{α,θ,θ+ϵ},PlotStyle->Red]],{{θ,π/6},0,π/2}]

This visualization tool shows the rings which one integrates over e.g. in Morin problem 5.86.
As you tackle that problem, it is a good idea to  express the area dA of such a little ring, and verify that if you integrate dA as θ ranges from 0 to π, you get the expected 4πR^2 surface area for a sphere.  Also note that the flux of particles is proportional to the vertical projection dA Cos[θ], so that e.g. the largest rings collect almost no particles, since they are sideways to the wind.
Verify that if you integrate dA Cos[θ] as θ ranges over the front half of the sphere (0 to only π/2), the total
projected area is π(R^2), i.e. the size of the shadow or “cross-section” for the sphere.
Then fold in the momentum transfer Δ(Subscript[p, x])=mv[1+Cos[2θ]) and integrate to find the force on the sphere.

And BTW here is 3D version of the movie.

num=100;
Do[(
x0[j]=Random[Real,{0,10}]; (* For each particle, *)
y0[j]=Random[Real,{-1,1}];(* pick a random starting place*)
z0[j]=Random[Real,{-1,1}];(* with x in {0,10}, and y and z in {-1,1}*)
b0[j] = Sqrt[y0[j]^2+z0[j]^2]; 
ϕ0[j]= ArcTan[y0[j],z0[j]];
θ0[j] = ArcSin[b0[j]];   (* Determine where on the disk we hit *)
t0[j] = x0[j]-Cos[θ0[j]]; (*And when it will hit*)
t0[j] = If[Element[t0[j],Reals],t0[j],1000];(*Some won't hit ever*)
pt[j,t_] = If[t<t0[j],{x0[j]-t,y0[j],z0[j]},{Cos[θ0[j]],Sin[θ0[j]]*Cos[ϕ0[j]],Sin[θ0[j]]*Sin[ϕ0[j]]}+(t-t0[j]){Cos[2θ0[j]],Sin[2θ0[j]]*Cos[ϕ0[j]],Sin[2θ0[j]]*Sin[ϕ0[j]]}];),
{j,1,num}] (*and put together the two pieces of the motion *)
Animate[Show[Graphics3D[{Blue,Sphere[]}],ListPointPlot3D[Table[pt[j,t],{j,1,num}],PlotStyle->Red],PlotRange->{{-1,5},{-1.5,1.5},{-1.5,1.5}}],{t,0,5}]

End[]

Morin 5.60 (Solar Escape)

The question is posed what velocity we need to have with respect to Earth so that a
satellite can escape not only Earth but the whole solar system, i.e. have positive total energy.
If we call the Earth’s orbital speed vO (=2πR/year), then we should launch in the direction
the Earth, and if that launch speed relative to Earth is vL, then the speed relative to the Sun
as we launch is (vO+vL).
Now here is the wrong logic: in the Sun’s frame we want the satellite to have total energy zero,
so the kinetic energy must balance the negative potential energies from the Earth and the Sun,
and we find vL by solving ...

Begin["Morin560`"]
Solve[(1/2)m (vO+vL)^2 - G Me m /Re - G Ms m/Rs==0,vL]
s=%[[2]];

Here is a package which defines all the constants we need.

Needs["PhysicalConstants`"]
vL/.s//. {G->GravitationalConstant,Me->EarthMass, Ms->SolarMass,Re->EarthRadius,Rs->AstronomicalUnit,vO->2π AstronomicalUnit/Year}//Convert[#,Meter/Second]&

Pay no attention to the warning message.  Here is the new and verbose way to express the same thing:

vL/.s/.{G->Quantity["GravitationalConstant"],Me->Quantity["EarthMass"], Ms->Quantity["SolarMass"],Re->
Quantity["EarthRadius"],Rs->Quantity["AstronomicalUnit"],vO->2π Quantity["AstronomicalUnit"]/Quantity["Year"]}
UnitConvert[%]

BTW what is wrong with the above logic?  And of course, how do we find the correct answer?

End[]

Effective Potential

Suppose a particle moves with no external force.
Here is the movie with coordinates drawn relative to some origin.

Begin["Veff`"];
state[x_,y_] := Show[Graphics[{
{Disk[{0,0},.1]},
{Red,Disk[{x,y},.1]},
{Dashed,Blue,Thin,Line[{{0,0},{x,0}}]},
{Dashed,Blue,Thin,Line[{{x,0},{x,y}}]},
{Thin,Green,Line[{{0,0},{x,y}}]},
{Green,Thin,Circle[{0,0},.2,{0,ArcTan[x,y]}]},
Text["r",{x/2,y/2+.2}],
Text["x",{x/2,-.2}],
Text["y",{x+.2,.5}],
Text["θ",{.3,.1}]
}],PlotRange->{{-4,4},{-1,2}}]
Animate[state[x,1],{x,-4,4}]

In the 1D effective theory, the effective potential has only the centrifugal term.
Here is the plot with a line drawn for an arbitrary particular energy.

Ueff[r_] =  L^2/ ( 2 m  r^2)
Plot[{Ueff[r]/.{L->1,m->1},1},{r,0,3},PlotStyle->{{Red},{Blue,Dashed}}]

In terms of L and E0 we can find the turning point, the minimum radius:

Solve[Ueff[r]==E0,r]
turningpoint= r/.%[[2]]

YOUR JOB: Express L and E0 in terms of the initial velocity v0 and impact parameter b,
and verify that the turning point is in fact exactly what you expect.

To find r[t] we need to DSolve this equation

rdiffeq =m r''==-Ueff'[r] /. {r''->r''[t],r->r[t]}

Can you think of a function where taking two derivatives is the same thing as raising it to -3 power?
Ask Mathematica for the solution, maybe choosing t=0 to the the time when r’[0]==0, and go “of course”.
Then with r[t] safely in hand, ask to DSolve for θ[t].  Recall the diff eqn:

θdiffeq = θ'[t] == L/(m r[t]^2)

End[]

Celestial Orbits

Start by writing down the effective potential, which gives the diffeq for r:

Begin["Celestial`"];
Ueff[r_] =  -α/r + l^2/ ( 2 m  r^2)
reqn = m r''[t] == - D[Ueff[r[t]],r[t]]

Plot for a choice of parameters

p = {l->1,m->1,α->1}
Plot[Ueff[r]/.p,{r,.2,3}]

The mininum of this potential is located at:

sr0 = Solve[Ueff'[r]==0,r][[1]]
er0 = Ueff[r]/.sr0

BTW if you modify the potential above, you may want to do this numerically:

sUmin = FindMinimum[Ueff[r]/.p,{r,.8}]
r0 = r /. sUmin[[2]]
e0 = sUmin[[1]]

Choose an energy for our motion somewhat bigger than the minimum,
and finally make our diagram:

e = e0+ .1 (* For ellipse*)
Plot[{Ueff[r]/.p,e},{r,.1 r0,3 r0}]

We might as well locate the turning points:

turn1 = (r/.FindRoot[(Ueff[r]/.p)==e,{r,.5r0}])
turn2 = (r/.FindRoot[(Ueff[r]/.p)==e,{r,2r0}])

Now to find the actual motion we need to set up initial conditions. 
Let us agree to launch at r[0]=r0 with whatever speed we need to make the energy right:

eeqn = (1/2) m(rdot)^2 + Ueff[r0] == e
rdot0 = (rdot/.(Solve[eeqn/.p,rdot])[[2]])

Now we are ready to solve, save a copy of the solution as “myr[t]”, and plot.

nsoln = NDSolve[{reqn/.p,r'[0]==rdot0,r[0]==r0},r[t],{t,0,100}]
myr[t_] = (r[t]/.nsoln[[1]])
Plot[myr[t],{t,0,30}]

Now that we have r[t], find θ[t].

θeqn = θ'[t] == l/(m (myr[t])^2);
θsoln = NDSolve[{θeqn/.p,θ[0]==0},θ[t],{t,0,100}]
myθ[t_] = θ[t]/.θsoln[[1]]
Plot[myθ[t],{t,0,30}]

And now we can draw the curve.  Note that if you change parameters above,
you might want to play with the aspect ratio and plot range here.

p1=ParametricPlot[{myr[t]*Cos[myθ[t]],myr[t]*Sin[myθ[t]]},{t,0,100},AspectRatio->1,PlotRange->{{-1.5,1.5},{-1,2}},PlotStyle->Red];
p2 = ListPlot[{{0,0},{0,1}},PlotStyle->{Black,PointSize[.02]}];
Show[p1,p2]

This tool lets you draw an ellipse with foci at (0,0) and (0,a).

Manipulate[Show[p1,p2,ContourPlot[(Sqrt[(x-0)^2 + (y-0)^2]+Sqrt[(x-0)^2+(y-a)^2])==1+Sqrt[1+a^2],{x,-1.5,1.5},{y,-1,2}]],{a,0,2}]

Now who can resist the temptation to make a movie?!

orbit[tt_] := Show[{ParametricPlot[{myr[t]*Cos[myθ[t]],myr[t]*Sin[myθ[t]]},{t,0,tt},DisplayFunction->Identity,PlotRange->{{-1.5,1.5},{-1,2}},AspectRatio->1],
ListPlot[{{myr[tt]*Cos[myθ[tt]],myr[tt]*Sin[myθ[tt]]}},PlotStyle->{PointSize[.03],RGBColor[1,0,0]},DisplayFunction->Identity]},DisplayFunction->$DisplayFunction]

Animate[orbit[t],{t,0,50}]

Your job: play!  Change the energy and get a hyperbola or a parabola.
 Also change the force law by adding some e.g. small 1/r^4 force.

End[];

Morin 7.2

Consider a force with associated potential

V[r_] = (-C/3)/r^3

We launch a particle from a point essentially at ∞ with some v0 and with some “impact parameter” b,
meaning our starting position is (x,b) for some very large x, and with velocity (-v0,0).
Among other things we then know our starting angular momentum and energy:

E0 = (1/2) m v0^2
L = m v0 b

Choosing parameters we can then plot the effective potential for our choice of b

Veff[r_] = L^2/(2 m r^2) + V[r]
p = {C->1,m->1,v0->1};
Manipulate[Plot[{Veff[r]/.p/.b->ImpactParam,E0/.p},{r,0,4},PlotRange->{-.4,.8}],{ImpactParam,1.1,1.3}]

Your job in the homework: for a given v0, find the condition on b which determines
whether the particle gets over the hump (and falls in to the singularity) or scatters off.
That is, determine the position and height of the peak by taking a derivative:

peakeqn =Veff'[r]==0
Solve[peakeqn,r]
peakE = Veff[r]/.%[[1]]

Then equate this to the initial energy (1/2) m v0^2 to determine the critical b.
Note that the critical value for our numerical parameters turns out to be (3)^(1/6).

To make trajectories, we need to solve the equation of motion for r.
We’ll pick some b=1.3 to start, and call r=10 close enough to ∞ for our purposes.

brule = b->1.3;
reqn = r''[t]==(-Veff'[r]/.p/.r->r[t]/.brule)
sr = NDSolve[{reqn,r[0]==10,r'[0]==-1},r[t],{t,0,15}]
myr[t_] = r[t]/.sr[[1]]
Plot[myr[t],{t,0,15}]

To make 2D plots we need θ[t], which we find from

θeqn = θ'[t] == ( L / (m myr[t]^2) /.p/.brule)

We want to start a distance b off of the x-axis, so the starting angle should some Arctan:

sθ = NDSolve[{θeqn,θ[0]==ArcTan[b/10]/.brule},θ[t],{t,0,15}];
myθ[t_]= θ[t]/.sθ[[1]]

Mathematica Homework: With r[t] and θ[t] in hand, make a ParametricPlot of the motion.  Change parameters and draw a trajectory which is “captured” and one which escapes to ∞.  These can be “Show”n on the same plot for handing in with the homework.

Morin 7.4

Let's find some orbits with a central force of the form β r^k.
Start with the effective potential.  If you defined L above, we clear it (and others) now with:

L=.;r=.;k=.
Ueff[r_] =  L^2/ ( 2 m  r^2) + β r^k
reqn = m r''[t] == - D[Ueff[r[t]],r[t]]

Can we DSolve?!?!  You may want to learn about the "Abort Evaluation" menu item under "Evaluation" above.

DSolve[reqn,r[t],t]

So let's choose specific numbers, plot and find the minimum, which we call r0.

p = {m->1,L->1,k->7,β->1}
Plot[Ueff[r]/.p,{r,.4,1.3}]
sUmin = FindMinimum[Ueff[r]/.p,{r,.8}]
r0 = r /. sUmin[[2]]
e0 = sUmin[[1]]

Choose an energy for our motion somewhat bigger than the minimum, and find the turning points:

e = e0+ .2
Plot[{Ueff[r],e}/.p,{r,.4,1.2}]
r1 = (r/.FindRoot[(Ueff[r]/.p)==e,{r,.5r0}])
r2 = (r/.FindRoot[(Ueff[r]/.p)==e,{r,2r0}])

Now to find the actual motion we need to set up initial conditions. 
Let us agree to launch at r[0]=r0 with the speed appropriate for the given energy.

eeqn = ((1/2) m(rdot)^2 + Ueff[r0] == e)/.p
rdot0 = (rdot/.(Solve[eeqn,rdot])[[2]])

Now we are ready to solve, save a copy and plot:

nsoln = NDSolve[{reqn/.p,r'[0]==rdot0,r[0]==r0},r[t],{t,0,10}]
myr[t_] = (r[t]/.nsoln[[1]])
Plot[myr[t],{t,0,10}]

Now that we have r[t], find θ[t].

θeqn = θ'[t] == (L/(m (myr[t])^2))/.p
θsoln = NDSolve[{θeqn,θ[0]==0},θ[t],{t,0,10}]
myθ[t_] = θ[t]/.θsoln[[1]]
Plot[myθ[t],{t,0,10}]

And now we can draw the curve.  Note that you might want to play with the aspect ratio.

p1=ParametricPlot[{myr[t]*Cos[myθ[t]],myr[t]*Sin[myθ[t]]},{t,0,10},AspectRatio->1,PlotRange->{{-1.5,1.5},{-1.5,1.5}}]

And of course we need a movie:

orbit[tt_] := Show[{ParametricPlot[{myr[t]*Cos[myθ[t]],myr[t]*Sin[myθ[t]]},{t,0,tt},DisplayFunction->Identity,PlotRange->{{-1.5,1.5},{-1,2}},AspectRatio->1],
ListPlot[{{myr[tt]*Cos[myθ[tt]],myr[tt]*Sin[myθ[tt]]}},PlotStyle->{PointSize[.03],RGBColor[1,0,0]},DisplayFunction->Identity]},DisplayFunction->$DisplayFunction]

Animate[orbit[t],{t,0,10}]

Mathematica Homework: 
The book claims (and you are to show in your own words in the homework) that for small amplitude oscillations (i.e. e close to e0) there are certain special values of k (e.g. k=7), for which the frequency of radial oscillation is a rational multiple of the orbital freqency, which means the orbit eventually closes.  Pick your own special value of “k”, a smallish energy, and ParametricPlot an approximately closed orbit for the homework.