Legal Notice, Disclaimer, and Book Ordering Information:

(* Legal Notice, Disclaimer, and Book Ordering Information:

Last Modified: FIO June 8, 2002

Mathematica For Physics: 2nd Edition

Robert L. Zimmerman and Fredrick I. Olness

WebSite:  
  www.physics.smu.edu/~olness
  darkwing.uoregon.edu/~phys600/

MathSource Number: 0206-862

ISBN 0-201-53796-6
For ordering information, call 1-800-282-0693
Addison-Wesley Publishing Company

Communication with the authors:
Fredrick I. Olness:  olness@mail.physics.smu.edu
Robert L. Zimmerman: bob@zim.uoregon.edu

Copyright 2002, Addison-Wesley Publishing Company, Inc. The material in
this file may be distributed freely  so long as the content remains
unchanged, and the copyright and  reference notices are included. The
publisher grants permission for the noncommercial use of these programs
and program segments.  All other uses require the prior written consent
of the publisher.    

All rights are reserved with regard to the  material in the
"Mathematica For Physics" text, and may not be reproduced, stored in a
retrieval system, or transmitted,  in any form or by any means,
electronic, mechanical, photocopying,  recording, or otherwise, without
the prior written permission  of the publisher. 

The programs and applications presented in this book have been included
for  their instructional value. They have been tested with care but are
not  guaranteed for any particular purpose. The publisher does not
offer any  warranties or representations, nor does it accept any
liabilities with  respect to the programs or applications. 

 *)

7. Orbiting Bodies

7.1 Introduction

Off[General::"spell"]; 
Off[General::"spell1"];

7.2 Two Body Problem

Overview of the Two-Body Problem

Problem 1: Equivalent One-body Problem

Remarks and outline

Solution

Clear["Global`*"]

Part a

eq1 = {M0 == m1 + m2, μ == ((m1 m2)/M0)};

mSol= Solve[eq1,{m1,m2}][[1]]

eq2 = {R[t] == ((m1 r1[t] + m2 r2[t])/M0), 
          r[t] == r1[t] - r2[t]};

sol1 = ((Solve[eq2, {r1[t], r2[t]}][[1]] //. 
                Rule -> Equal //. {M0 -> m1 + m2}) // 
          MapAll[Apart, #] &)

rSol = {r1 -> (((m2 r[#1])/(m1 + m2)) + R[#1] &), 
          r2 -> (-((m1 r[#1])/(m1 + m2)) + 
                  R[#1] &)};

T=(1/2) m1 ((r1'[t])^2)+(1/2) m2 ((r2'[t])^2);

La =  T - V[r0[t]] /.  rSol /. mSol // Simplify

Needs["Calculus`VectorAnalysis`"]

rRule = 
      r -> ((CoordinatesToCartesian[{r0[#], θ[#], ϕ[#]}, 
                  Spherical] // Evaluate) &)

RRule = {R -> ({Rx[#], Ry[#], Rz[#]} &)};

r/:VectorQ[r]=True; 
R/:VectorQ[R]=True;

vectorRule = 
{x_[   t_]^2 :> 
            x[   t] . x[   t] /; 
              VectorQ[x], 
x_'[
                t_]^2 :> x'[t] . x'[t] /; 
              VectorQ[x]};

La //.vectorRule

Lb= La //.vectorRule //.rRule //.RRule //Simplify

Needs["Calculus`VariationalMethods`"]

EulerEquations[Lb,{Rx[t],Ry[t],Rz[t]},t]

L = 
      Lb /. {Rx -> (0 &), Ry -> (0 &), 
            Rz -> (0 &)} // Simplify

eqMotion = (EulerEquations[Lb, {r0[t], θ[t], ϕ[t]}, 
                    t] // Solve[#, {r0''[t], θ''[
                          t], ϕ''[t]}] &  // Flatten  // 
              Simplify ) /. 
          Rule -> Equal;
eqMotion // ColumnForm

Part b

θRule = θ -> (π/2 &);

eqMotion =  eqMotion //. θRule;
eqMotion // ColumnForm

dSol = DSolve[{eqMotion[[3]], 
              r0[0] == Sqrt[el/(μ ϕ'[0])]}, {r0[t]}, 
            t] // PowerExpand //Flatten

ϕRule= Solve[dSol /.{Rule->Equal} ,ϕ'[t]][[1]]

rEq=eqMotion[[1]] /.ϕRule //ExpandAll

energy= Integrate[(rEq[[1]]-rEq[[2]]) r0'[t],t]

eq3 = FirstIntegrals[L /.  θRule, {r0[t],   ϕ[t]} , 
                t] /. {FirstIntegral[ϕ] -> -el , 
                FirstIntegral[t] -> En  μ} /. 
            Rule -> Equal // ExpandAll;
eq3 // ColumnForm

Solve[eq3 , {ϕ'[t], En}] // Flatten // ExpandAll // 
      ColumnForm

Part c

uRule = {r0 -> ((1/u)[ϕ[#1]] &)};

ϕRule1 = {ϕRule, D[ϕRule, t]}  //. uRule  // 
        Flatten

eq3 = ( Solve[rEq //. uRule //. ϕRule1, u''[ϕ[
                      t]]][[1]] )  /. 
            Rule -> Equal /. ϕ[t] -> ϕ // 
        Simplify

eq4 = eq3//.{V'[
              x_] -> -f[x]}

Problem 2: Kepler Orbits

Remarks and outline

Clear[diffSeriesOne];

    
diffSeriesOne[eqin_, z_, t_, order_, 
          initList__List] := 
Module[{eq, zSer, eSer, eqs, 
            vars, sol, zSol}, 
eq = 
            eqin[[1]] - eqin[[2]] == 
              0; 
zSer = 
            Series[z[t], {t, 0, order}]; 
eSer = 
            eq /. (z[t] -> zSer) // 
              ExpandAll; 
eqs = 
            Thread[CoefficientList[eSer[[1]], t] == 
                0]; 
vars = (Table[
                  D[z[t], {t, j}], {j, 2, order + 2}] /. 
                t -> 0);      sol = Solve[eqs, 
                vars][[1]]; 
zSol = zSer /. sol /. 
              initList; 
Return[
            zSol]
];

    
Protect[diffSeriesOne];

Solution

Clear["Global`*"]

Part a

eq1=u''[ϕ]==-((m f[(1/u)[ϕ]])/((el^2) u[ϕ]^2)) - 
            u[ϕ];

forceRule={\(V -> \((\(-\(k\/#1\)\) &)\)\),f[(1/u)[ϕ]]->-V'[r],r -> (1/u)[ϕ]};

eq2= eq1 //.forceRule

eq3 = DSolve[{eq2, u'[0] == 
                  0 }, {u[ϕ]}, ϕ][[1]] // ExpandAll // 
        Simplify

eq4 = eq3 /. C[1] -> ((e k m)/(el^2)) // Simplify

eq5 = eq4 //. {el -> (Sqrt[a (1 - (e^2)) k m])}

Part b

Needs["Graphics`Graphics`"]

?PolarPlot

plotOrbit[a_, e_, opts___] := 
        PolarPlot[
          a (1 - e^2)/(1 + e Cos[ϕ]), {ϕ, 0, 2 π}, 
          opts, DisplayFunction -> Identity];

list = 
        Table[plotOrbit[1, 
            eIn , 
PlotStyle ->  {Thickness[0.008], 
                Hue[eIn ]}],  {eIn, 0.1, 0.9, 0.2} ];

Show[ list , DisplayFunction -> $DisplayFunction];

plotOrbit[1,1.3,DisplayFunction->$DisplayFunction];

Part c

eq2

initial={u'[0]->0,u[0] -> (((1 + e) k m)/(el^2))};

seriesSol = diffSeriesOne[eq2, u, ϕ, 6, initial] // Normal // 
        ExpandAll

cterm[
          n_] = (((((-1))^n) e k m (ϕ^(2 n)))/((el^2) (2 n)!));

seriesSol == (((k m)/(el^2)) + 
          Sum[cterm[i], {i, 0, 3}])

answer = ((k m)/(el^2)) + (∑_(n = 0)^∞) cterm[
                n]   // PowerExpand // Simplify

(u[ϕ]  /. eq4 ) == answer

Problem 3: Precessing Ellipse

Remarks and outline

firstOrderPert[eq1_, x_, {x0_, v0_}, eps_:  eps, 
          t_:  t] := 
Module[{xRule, eq2, eqEps, initCon0, 
            initCon1, dsol1, 
            dsol2}, 
xRule = {x -> (x[0][#] + 
                      eps x[1][#] &)}; 
eq2 = 
            eq1[[1]] - eq1[[2]] /. xRule; 
eqEps = 
            CoefficientList[eq2, eps] == 0 // 
              Thread; 
initCon0 = {x[0][0] == 
                x0, (x[0])'[0] == 
                v0}; 
initCon1 = {x[1][0] == 
                0, (x[1])'[0] == 
                0}; 
dsol1 = (DSolve[
                  Join[{eqEps[[1]]}, initCon0], {x[0][t]}, t] // 
                Flatten); 
dsol2 = (DSolve[
                  Join[{eqEps[[2]]}, initCon1] //. dsol1, {x[1][t]}, 
                  t] // Flatten); 
Return[x[t] /. 
                xRule //. 
              Join[dsol1, 
                dsol2]]
];

    
Protect[firstOrderPert];

Solution

Clear["Global`*"]

Part a

eq1=u''[ϕ]==-((m f[(1/u)[ϕ]])/((el^2) u[ϕ]^2)) - 
            u[ϕ];

forceRule={\(V -> \((\(-\(k\/#1\)\) - b\/#1\^2 &)\)\),f[(1/u)[ϕ]]->-V'[r],r -> (1/u)[ϕ]};

eq2 = eq1 //. forceRule // ExpandAll

eq3 =  (DSolve[{ eq2, u'[0] == 0} , 
              u[ϕ], ϕ])  // Simplify  // Flatten

CRule =  
        C[1] ->  ((e k m)/( (el^2) - 2 b m));
uSol =  eq3 /. CRule // Simplify

shiftRule={((Sqrt[(el^2) - 2 b m])/el) -> 1 - Δ};

eq3 /. CRule /. shiftRule // FullSimplify

elRule = Solve[{((k m)/((el^2) - 2 b m)) == 
                (1/(a (1 - (e^2))))}, el][[2]]  // 
        Simplify

eq3 /. CRule /. shiftRule /. elRule // Simplify

Part b

plotorbit[start_, finish_, {a_, e_, Δ_}, opts___ ] := 
      ParametricPlot[ {((a (1 - (e^2)) Cos[ϕ])/(1 + e Cos[(1 - Δ) ϕ])), ((a (1 - (e^2)) Sin[ϕ])/(1 + e Cos[(1 - Δ) ϕ]))} // Evaluate , {ϕ, start, 
          finish}, opts ]

values = { 1, 0.5, 0.1  };

plotorbit[0, 6 π, values, 
        PlotStyle -> {Thickness[0.008], Hue[0.8]}];

Part c

eq2

eq4 = firstOrderPert[eq2, 
                u, {(((1 + e) k m)/((el^2) - 2 b m)), 0}, 
                b, ϕ]   + O[b]^2  // Normal // Simplify  // 
        Expand

u[ϕ]/. uSol

eq5 = ( u[ϕ] /.  uSol) + O[b]^2  // Normal // 
          ExpandAll//PowerExpand

eq4 == eq5 // Simplify

Problem 4: Numerical Solution

Remarks and outline

Solution

Clear["Global`*"]

Part a

eq1=u''[ϕ]==-((m f[(1/u)[ϕ]])/((el^2) u[ϕ]^2)) - 
            u[ϕ];

forceRule={\(V -> \((\(-\(k\/#1\)\) &)\)\),f[(1/u)[ϕ]]->-V'[r],r -> (1/u)[ϕ]};

eq2= eq1 //. forceRule

values = {k -> 1, el -> 1, m -> 1};

initial1={u[0]==1.6, u'[0]==0};

ndsol1 = 
        NDSolve[{eq2 /. values, initial1} // Flatten, 
            u[ϕ], {ϕ, 0, 2 π}] // Flatten;

u1[ϕ_]=u[ϕ]/.ndsol1;
u1[1.1]

Needs["Graphics`Graphics`"]

PolarPlot[
        (1/u1)[ϕ] // Evaluate, {ϕ, 0, 2 π}];

initial2= {u[0]==3, u'[0]==0};

ndsol2 = 
        NDSolve[{eq2 /. values, initial2} // Flatten, 
            u[ϕ], {ϕ, 0, 2 π}] // Flatten;

u2[ϕ_]=u[ϕ] //.ndsol2;
u2[1.1]

PolarPlot[
        (1/u2)[ϕ] // Evaluate, {ϕ, 0, 2 π}];

Part b

numOrbit[V_, r_, {r0_:  1, dr0_:  0, el_:  1, m_:  1}, angle_ , 
        opts___] := 
      Module[{u, ϕ, uforce, eq1, initial, eqs, 
          ndsol}, 
uforce = -D[V, r] /. 
            r -> 
              1/u[ϕ]; 
eq1 = (u''[ϕ] == -u[ϕ] - m uforce/(el^2 u[ϕ]^2) // 
              ExpandAll); 
initial = {u[0] == 
              1/r0, u'[0] ==  -dr0/
                r0^2 }; 
eqs = 
          Join[{eq1}, initial]; 
ndsol = 
          NDSolve[eqs, u[ϕ], {ϕ, 0, angle}] // 
            Flatten;   PolarPlot[
          1/u[ϕ] /. ndsol // Evaluate, {ϕ, 0, angle}, 
          opts, 
DisplayFunction -> 
            Identity]
]

numOrbit[-(1/r), r, { .6, 0, 1, 1}, 2 π, 
        PlotStyle -> {Thickness[ .008], Hue[0.6]}, 
        DisplayFunction -> $DisplayFunction];

lista = numOrbit[#, r, {0.6, 0, 1, 1}, 6 π, 
              PlotStyle -> {Thickness[0.008], 
                  Hue[0.6]} ] & /@  {-(1/(r^1.1))`, -(1/(r^0.9))`};

Show[GraphicsArray[lista]];

\(\(listb = \(numOrbit[#, r, {0.6, 0, 1, 1}, 6\ π, 
              PlotStyle -> {Thickness[0.008], 
                  Hue[0.6]}\ ] &\) /@ {\(-\(0.02`\/r\^3\)\) - 1\/r, 
            0.02`\/r\^3 - 1\/r};\)\)

Show[GraphicsArray[listb]];

7.3 Restricted Three Body Problem

Overview of the Three Body Problem and Initialization of PMotion and Mgraph

PMotion[μ1_, {x0_, y0_, vx0_, vy0_}, tmax_, 
        step_:  1000] := 
Module[{rrule, Urule, eqMotion, 
          r1, r2, μ, 
          U}, 
rrule = {r1 -> (Sqrt[(((μ + x[t]))^2) + y[t]^2]), 
            r2 -> (Sqrt[(((-1 + μ + x[t]))^2) + y[t]^2])}; 
Urule = {U -> ((1 - μ)/r1) + (μ/r2) + 
                (1/2) (x[t]^2 + 
                      y[t]^2)}; 
eqMotion = 
{x''[t] - 2 y'[t] == 
              D[U /. Urule /. rrule, x[t]], y''[t] + 
                2 x'[t] == 
              D[U /. Urule /. rrule, 
                y[t]]}; 
NDSolve[{eqMotion /. μ -> μ1, x[0] == x0, 
                y[0] == y0, x'[0] == vx0, y'[0] == 
                  vy0} // Flatten, {x, y}, {t, 0, tmax}, 
            MaxSteps -> step] // Flatten]

Mgraph[fx_, fy_, t_, tfinal_, 
          Opts___] := 
ParametricPlot[{fx[t], fy[t]} // 
            Evaluate, 
{t, 0, 
            tfinal}, 
Opts, 
AspectRatio -> Automatic, 
DisplayFunction -> 
            Identity];

Protect[PMotion,Mgraph];

Problems on the Equal Mass Primaries (μ=1/2 )

Problem 1: Lagrangian Points for Equal Mass Binaries (μ=1/2 )

Solution

Clear["Global`*"]

Part a

Urule = 
        U -> ((1 - μ)/r1) + (μ/r2) + 
            (1/2) (x[t]^2 + y[t]^2);

rrule = {r1 -> (Sqrt[(((μ + x[t]))^2) + y[t]^2]), 
          r2 -> (Sqrt[(((-1 + μ + x[t]))^2) + y[t]^2])};

(eq1 = D[U  /. Urule /. rrule, x[t]] == 0 ) // 
      TraditionalForm  

    
(eq2 = D[U  /. Urule /. rrule, y[t]] == 0 ) // 
      TraditionalForm

{L[2], L[3], L[4], L[5], 
          L[1]} = {eq1, eq2}  /. μ -> (1/2)  // 
          NSolve[#, {x[t], y[t]}] &;
{L[2], L[3], L[4], L[5], L[1]}  // ColumnForm

CLpt = L[#] &  /@  Range[5];

{x[t], y[t]}  /. {L[1 ], L[2], L[3]}  // 
      TableForm[#, 
          TableHeadings -> { {"L[1]", "L[2]", "L[3]"} , { 
                x[t], y[t]} }] &

{L[4], L[5]} // ColumnForm

r1 == r2 == 1 /. rrule /. 
        L[4] /. μ -> 1/2

\(\(Lgraph = \ 
        Graphics[{{PointSize[0.02], 
              Hue[0.9], \ \ Map[
                Point, {x[t], y[t]} /. 
                  CLpt\ ]}\ , 
\ \ \ \ \ \(Text[\ 
                  L\_#\ \ , {x[t] + 0.1, y[t]} /. CLpt[\([#]\)]\ ] &\) /@ \ 
              Range[5]}\ ];\)\)

\(\(equilateral = \ 
        Graphics[{{PointSize[0.02], \ \ \ Point[{\(-\ 1\)/2, 0}\ \ ], \ 
              Point[{\ 1/2, 0}\ \ ]}, 
{Thickness[0.002], 
              Hue[ .4], \ \ \ Line[\ {{\(-1\)/2, 0}\ , {x[t], y[t]} /. 
                    L[4], {1/2, 0}, {\(-1\)/2, 0}, {x[t], y[t]} /. 
                    L[5]\ \ \ , {\ 1/2, 0}}\ ]\ \ \ \ \ }, \(Text[\ 
                  L\_#\ \ , {x[t] + 0.1, y[t]} /. CLpt[\([#]\)]\ ] &\) /@ 
              Range[5]}\ ];\)\)

Show[Lgraph, equilateral, 
        AspectRatio -> Automatic];

Part b

eq3 = r1 == r2 == 1   /. rrule

eq4=Solve[eq3, {x[t], y[t]} , 
        InverseFunctions -> True]

{eq1, eq2} /. eq4//Simplify

Part c

jacobi =C== -( ((x'[t])^2)+((y'[t])^2))+2  U  /.Urule

jacGC = 
        ContourPlot[ 2  U   /. Urule /. rrule /. μ -> 
                1/2 // Evaluate, 
{x[t], -2, 
            2} , {y[t], -2 , 
            2 },  
ContourShading -> 
            False, 
Contours -> 
            Range[3, 4.5 , 0.1], 
PlotPoints -> 
            66, 
DisplayFunction -> 
            Identity ];

Show[jacGC, Lgraph, equilateral, 
        DisplayFunction -> $DisplayFunction ];

Problem 2: Looping Motion in an Equal Mass Binary System (μ=1/2)

Solution

Clear["Global`*"]

Part a

sol = 
        PMotion[1/2,  {0.3, 0, 0, 1.781 }, 6 π , 1000];

fx[t_] := x[t] /. sol;  
fy[t_] := y[t] /. sol ;

pt1 =   Mgraph[fx, fy, t, 
          6 π  , 
PlotPoints -> 
            100, 
PlotStyle -> 
            Hue[0.9] ,  
AspectRatio -> 
            Automatic, 
Epilog -> {AbsolutePointSize[
                7], Point[{1/2, 
                  0} ]}, 
DisplayFunction -> $DisplayFunction];

Show[pt1, PlotRange -> {{0.692, 0.709}, {-0.01, 0.01 }}, 
        Ticks -> False];

ParametricPlot[{{(1/2) Cos[t], 
              (1/2) Sin[t] }, {fx[t] Cos[t] - fy[t] Sin[t], 
              fx[t] Sin[t] + fy[t] Cos[t]}} // Evaluate, {t, 
          0,    2 π}, 
PlotStyle -> {GrayLevel[
              0.3], Hue[0.6]} , 
AspectRatio -> 
          Automatic];

Part b

sol=PMotion[1/2,  { 1.35 , 0 , 0, - 1.568 }, 
        10 π , 1000]

fx[t_] := x[t] /. sol;  
fy[t_] := y[t] /. sol ;

Mgraph[fx, fy, t, 
        10 π , 
PlotPoints -> 
          100, 
PlotStyle -> 
          Hue[0.8], 
PlotRange -> {{ -0.5 , 
              1.5 }, {-1, 
              1}}, 
Epilog -> {AbsolutePointSize[7], 
            Point[{1/2, 
                0} ] }, 
DisplayFunction -> $DisplayFunction];

pt1 = 
ParametricPlot[
{{(1/2) Cos[t], (1/2) Sin[t] }, 1.35 {  Cos[t],  Sin[t] }, 
                  1.1 {  Cos[t],  
                      Sin[t] }, 
{fx[t] Cos[t] - 
                      fy[t] Sin[t], fx[t] Sin[t] + fy[t] Cos[t]}} // 
                Evaluate, 
{t, 
                0, # π}, 
PlotStyle -> {GrayLevel[
                    0.3], Hue[0.6], 
                  Hue[0.9]}, 
AspectRatio -> 
                Automatic, 
DisplayFunction -> 
                Identity
] &  /@  {2, 10};

Show[GraphicsArray[pt1]];

Part c

sol=PMotion[1/2, {-0.01112 , 0, 0,  +1.538 }, 
        10 π, 1000]

fx[t_] := x[t] /. sol ; 
fy[t_] := y[t] /. sol ;

Mgraph[fx, fy, t, 
        10  π , 
PlotPoints -> 
          100, 
PlotStyle -> 
          Hue[0.6]  , 
Epilog -> {AbsolutePointSize[
              7], Point[{1/2, 
                0} ]}, 
DisplayFunction -> $DisplayFunction];

Problem 3: Symmetric orbits about the y-axis for μ=1/2

Solution

Clear["Global`*"]

sol = 
        PMotion[1/2, { 0, 0.12067 , -1.0600745428 , 0 }, 6 π , 
          1000] ;

fx[t_] := x[t] /. sol ;  
fy[t_] := y[t] /. sol ;

CLpt = 
{{y[t] -> 0. , 
            x[t] -> 0. }, 
{y[t] -> 0. , 
            x[t] -> -1.1984}, 
{y[t] -> 0. , 
            x[t] -> 1.1984}, 
{y[t] -> 0.8660, 
            x[t] -> 0. }, 
{y[t] -> -0.8660, 
            x[t] -> 0. }};

Mgraph[fx, fy, t, 
        6 π, 
PlotPoints -> 
          100, 
PlotStyle -> 
          Hue[ .8],  
Epilog -> {AbsolutePointSize[
              7], Point[{1/2 , 0} ], 
            Point[{ -1/2, 
                0} ], Point[{x[t], y[t]} /. CLpt[[#]] ] &  /@  
              Range[4] } , 
DisplayFunction -> $DisplayFunction];

Problem 4: Mass exchange Between Equal Mass Binaries

Solution

Clear["Global`*"]

Part a

sol = 
        PMotion[1/2, { 0.044  , 0 , 0 , 0  }, 8 , 
          1000] ;

fx[t_] := x[t] /. sol;  
fy[t_] := y[t] /. sol ;

pt1 =  
        Mgraph[fx, fy, t,  
          8 , 
PlotPoints -> 
            100, 
PlotStyle -> 
            Hue[ .8],  
Epilog -> {AbsolutePointSize[
                7], Point[{0, 0}] , 
              Point[{1/2, 
                  0}]}, 
DisplayFunction -> $DisplayFunction];

jacobi =C== -( ((x'[t])^2)+((y'[t])^2))+2  U ;

Urule = 
        U -> ((1 - μ)/r1) + (μ/r2) + 
            (1/2) (x[t]^2 + y[t]^2);

rrule = {r1 -> (Sqrt[(((μ + x[t]))^2) + y[t]^2]), 
          r2 -> (Sqrt[(((-1 + μ + x[t]))^2) + y[t]^2])};

( jacobi /. Urule /. rrule /. μ -> (1/2) /. 
                y[t] -> 0 /. x'[t] -> 0  /. 
            x[t] -> 0.044  /. y'[t] -> 
            0  )[[2]]

pt2 =  
        ContourPlot[ 2  U   /. Urule /. rrule /. μ -> 
                1/2 // Evaluate, 
{x[t], -1, 
            1} , {y[t], - .5 ,  .5 }, 
 
          ContourShading -> 
            False, 
AspectRatio -> 
            Automatic, 
Contours -> {  4.033,  
              3.9 , 3.5 },  
PlotPoints -> 
            66,   
Epilog ->  
            pt1[[1]]];

Part b

sol = 
        PMotion[1/2,  { 0.044  , 0 , 0.6  , 0   }, 35 , 
          4000] ;

fx[t_] := x[t] /. sol  ;
fy[t_] := y[t] /. sol;

pt1  =  
        Mgraph[fx, fy, t,  35 , PlotPoints -> 100, 
          PlotStyle -> Hue[0.8]   ];

pt2 =  
        ContourPlot[ 2  U   /. Urule /. rrule /. μ -> 
                1/2 // Evaluate, 
{x[t], -1.1, 
            1.1} , {y[t], -0.6 , 0.6 }, 
 
          ContourShading -> 
            False, 
AspectRatio -> 
            Automatic, 
Contours -> {4, 3.6105,  
              3.5,  3.4 , 3.3 },  
PlotPoints -> 
            66,   
Epilog ->  
            pt1[[1]]];

Problems on the Sun-Jupiter System (μ=.000954)

Problem 5:  Lagrangian Points for the Sun-Jupiter System

Solution

Clear["Global`*"]

Part a

Urule = 
          U -> ((1 - μ)/r1) + (μ/r2) + 
              (1/2) (x[t]^2 + y[t]^2);
 jacobi =C== -( ((x'[t])^2)+((y'[t])^2))+2  U  /.Urule;

rrule = {r1 -> (Sqrt[(((μ + x[t]))^2) + y[t]^2]), 
          r2 -> (Sqrt[(((-1 + μ + x[t]))^2) + y[t]^2])};

jacGJ = 
        ContourPlot[ 2  U   /. Urule /. rrule /. μ -> 
                0.000954 // Evaluate, 
{x[t], -1.2, 
            1.2} , {y[t], -1.2 , 
            1.2 } ,  
ContourShading -> 
            False , 
Contours -> 
            Range[3, 3.05 , 0.01], 
PlotPoints -> 
            150 ];

Part b

{L[5], 
        L[4]} = {{x[t] ->  (1/2) (1 - 2 μ), 
            y[t] -> -((Sqrt[3])/2)}, {x[t] -> 
              (1/2) (1 - 2 μ), 
            y[t] ->  ((Sqrt[3])/2)}} /. {μ -> 0.000954 }

eq1 = D[U   /. Urule /. rrule, x[t]] == 
            0  /.  μ ->  .000954 /. y[t] ->  
          0

        
eq2 = D[U   /. Urule /. rrule, y[t]] == 
            0  /.  μ ->  .000954 /. y[
            t]->0

eq3 = FindRoot[eq1[[1]] /.  μ -> 0.000954 /. 
              y[t] -> 0, {x[t],  # }] &  /@  {1, 
          0.9, -1 }

{L[3], L[1], 
        L[2]} = 
Outer[ List, 
            eq3,  {y[t] -> 0} ] // 
          Flatten//Partition[#, 2] &

JLpt = L[#] &  /@  Range[5];
JLpt  // ColumnForm

\(\(Lgraph = 
Graphics[
{PointSize[
                0.02], {Point[{μ, 0}], 
                Point[{\ 1 - μ, 0}]}, {PointSize[0.02], Hue[0.9], 
                Map[Point, {x[t], y[t]} /. JLpt\ ]}, \(Text[\ 
                    L\_#\ \ , {x[t]\ , y[t] + 0.1} /. \ 
                      JLpt[\([#]\)]\ ] &\)\  /@ \ \ Range[
                  5]\ \ 
}\ 
] /. μ -> 0.000954\ ;\)\)

Show[jacGJ,  Lgraph ];

Problem 6: Numerical Solution for the Trojan Asteroids

Solution

Clear["Global`*"]

Part a

initial =  { R Cos[θ],  R Sin[θ], 0, 0} /. 
          R -> (1/2) (Sqrt[3 + (((1 - 2 μ))^2)]) /. μ ->  
          0.000954

sol = 
        PMotion[ 0.000953875 , 
          initial /. θ -> 35 Degree   // Flatten, 200 , 
          1000];

fx[t_] := x[t] /. sol ;
fy[t_] := y[t] /. sol ;

pt1 = 
        Mgraph[fx, fy, t, 200 , PlotPoints -> 100, 
          PlotStyle -> Hue[0.8]];

Urule = 
        U -> ((1 - μ)/r1) + (μ/r2) + 
            (1/2) (x[t]^2 + y[t]^2);

rrule = {r1 -> (Sqrt[(((μ + x[t]))^2) + y[t]^2]), 
          r2 -> (Sqrt[(((-1 + μ + x[t]))^2) + y[t]^2])};

jacGJ = 
        ContourPlot[ 2  U   /. Urule /. rrule /. μ -> 
                0.000954 // Evaluate, 
{x[t], -1.2, 
            1.2} , {y[t], -1.2 , 
            1.2 } ,  
ContourShading -> 
            False , 
Contours -> 
            Range[3, 3.05 , 0.01], 
PlotPoints -> 
            150 , 
DisplayFunction -> 
            Identity];

JLpt = 
{{x[t] -> 0.9324, 
            y[t] -> 0}, 
{x[t] -> -1.0004, 
            y[t] -> 0}, 
{x[t] -> 1.0688 , 
            y[t] -> 0}, 
{x[t] -> 0.4990, 
            y[t] -> ((Sqrt[3])/2)}, 
{x[t] -> 0.4990, 
            y[t] -> -((Sqrt[3])/2)}};

\(\(Lgraph = 
Graphics[
{PointSize[
                0.02], {Point[{μ, 0}], 
                Point[{\ 1 - μ, 0}]}, 
{PointSize[
                  0.02], Hue[0.9], 
                Map[Point, {x[t], y[t]} /. JLpt\ ]}, \(Text[\ 
                    L\_#\ \ , {x[t]\ , y[t] + 0.1} /. \ 
                      JLpt[\([#]\)]\ ] &\) /@ \ 
                Range[5]\ \ 
}\ 
]\ \  /. μ ->  .000954\ ;\)\)

Show[jacGJ , pt1, 
        Lgraph, 
PlotRange -> {{-0.3, 1.2}, {0, 
              1.2}} , 
DisplayFunction -> $DisplayFunction];

Show[
        pt1, 
PlotRange  -> {{0.5, 0.85}, {0.55, 
              0.9}}, 
Ticks -> 
          False, 
DisplayFunction -> $DisplayFunction];

Part b

sol = 
        PMotion[ 0.000953875 , 
          initial /. θ ->  20 Degree   // Flatten, 200 , 
          1000] ;

fx[t_] := x[t] /. sol ; 
fy[t_] := y[t] /. sol ;

Mgraph[fx, fy, t, 
        160 , 
PlotPoints -> 
          100, 
PlotStyle -> 
          Hue[0.8] , 
Epilog -> 
          Lgraph[[1]], 
DisplayFunction -> $DisplayFunction ];

Problem 7: Perturbative Solution for the Trojan Asteroids

Solution

Clear["Global`*"]

Part a

rrule = {r1 -> (Sqrt[(((μ + x[t]))^2) + y[t]^2]), 
          r2 -> (Sqrt[(((-1 + μ + x[t]))^2) + y[t]^2])};

Urule = 
        U -> ((1 - μ)/r1) + (μ/r2) + 
            (1/2) (x[t]^2 + y[t]^2);

eqMotion = 
{x''[t] - 2 y'[t] == 
            D[U /. Urule /. rrule , 
              x[t]], 
y''[t] + 2 x'[t] == 
            D[U /. Urule /. rrule , y[t]]};

    
eqMotion //  ColumnForm // TraditionalForm

pert = 
{y  -> (y0 + ηy[#] &), 
x  -> (x0 + ηx[#] &)} ;

zeropt = {x0 -> (1/2) (1 - 2 μ), 
          y0 -> ((Sqrt[3])/2)};

eq1 = eqMotion  /. pert;

eqPer = ((MapAt[ 
                        Series[# , {ηx[t] , 0, 1}, {ηy[t] , 0, 
                              1} ] & , eq1[[#]] /. zeropt , 2] // 
                      Normal // 
                    ExpandAll)   /.  ηx[t] ηy[t] -> 
                  0  
) &   /@   {1, 
            2};

    
eqPer // ColumnForm

Part b

perSol = { ηx[t], ηy[t]} /. 
            DSolve[{eqPer, ηx[0] == xi,  ηy[0] == 
                    yi, 
ηx'[0] == 
                    0,  ηy'[0] == 0} // 
                Flatten,  {ηx[t], ηy[t]}, t] // 
          Flatten ;

Short[perSol, 4]

{Ω1, Ω2} = {(Sqrt[-(1/2) - (1/2) (Sqrt[1 - 27 μ + 27 (μ^2)])]), (Sqrt[-(1/2) + (1/2) (Sqrt[1 - 27 μ + 27 (μ^2)])])};

{Ω1, Ω2} /. μ ->  
        0.000954

Part c

Abs[{Ω1, Ω2} /. μ ->  
          0.000954]

11.86 { ((1 )/0.997) , ((1 )/0.0805)  } yrs

Part d

zeropt /. μ->0.000954

trojanpert = {0.499  + perSol[[1 ]], ((Sqrt[3])/2) + 
                      perSol[[2]]} /. {xi -> 0.0001, 
                    yi -> 0.0001} /. μ -> 0.000954  // 
              ExpandAll // ExpToTrig  // Chop;

    
trojanpert // Short[#, 6] &

pt1 = 
        ParametricPlot[ 
          trojanpert // Evaluate,  {t, 0, 80} ];

nSol = 
        PMotion[ 
          0.000954 ,  { 0.499  + 0.0001,   ((Sqrt[3])/2) + 0.0001, 0   ,  
              0  } // Flatten, 80 , 1000];

fx[t_] := x[t] /. nSol;  
fy[t_] := y[t] /. nSol ;

Mgraph[fx, fy, t, 80 , 
PlotPoints -> 100, 
        PlotStyle -> {Hue[0.8], Thickness[0.01], 
            Dashing[{0.03}]} , 
Epilog -> 
          pt1[[1]], 
DisplayFunction -> $DisplayFunction  ];

Problems on the Earth-Moon System (μ= .01215)

Problem 8:  Lagrangian Points for the Earth-Moon System

Solution

Clear["Global`*"]

Part a

Urule = 
        U -> ((1 - μ)/r1) + (μ/r2) + 
            (1/2) (x[t]^2 + y[t]^2);

jacobi =C== -( ((x'[t])^2)+((y'[t])^2))+2  U  /.Urule

rrule = {r1 -> (Sqrt[(((μ + x[t]))^2) + y[t]^2]), 
          r2 -> (Sqrt[(((-1 + μ + x[t]))^2) + y[t]^2])};

jacGE = 
        ContourPlot[ 2  U   /. Urule /. rrule /. μ -> 
                0.01215  // Evaluate, 
{x[t], -1.2, 
            1.2} , {y[t], -1.2 , 
            1.2 },  
ContourShading -> 
            False, 
Contours -> 
            Range[3, 3.55 , 0.05], 
PlotPoints -> 
            150 ];

Part b

eq1 = D[U  /. Urule /. rrule, 
          x[t]]==0 
 eq2 = D[U  /. Urule /. rrule, 
          y[t]]==0

eq3 = FindRoot[eq1[[1]] /. μ -> 0.01215  /. 
              y[t] -> 0, {x[t],  # }] &   /@ {1, 
          0.9, -1 }

{L[1], L[2], 
        L[3]} = 
Outer[ List, 
            eq3,  {y[t] -> 0} ] // 
          Flatten//Partition[#, 2] &

eq4 = 
{{x[t] -> (1/2) (1 - 2 μ), 
            y[t] -> -((Sqrt[3])/2)}, 
{x[t] -> 
              (1/2) (1 - 2 μ), y[t] -> ((Sqrt[3])/2)}};

{L[4], L[5]} = 
      eq4 /. μ->0.01215

LptE = L[#] &  /@  Range[5] ;

Lgraph = 
Graphics[
{PointSize[ .02], {Point[{μ, 0}], 
                Point[{ 1 - μ, 
                    0}]}, 
{PointSize[ .02], Hue[ .9], 
                Map[Point, {x[t], y[t]} /. 
                    LptE]} 
} 
]   /. μ -> 0.01215;

Show[jacGE,  Lgraph ];

Problem 9: Motion about L[4] in the Earth-Moon System

Solution

Clear["Global`*"]

eqPer=
{-2 (ηy)'[t]+(ηx)''[t]==((3 ηx[t])/4) + 
                (3/4) (Sqrt[3]) (1 - 2 μ) ηy[t],2 (ηx)'[t]+(ηy)''[t]==(3/4) (Sqrt[3]) (1 - 2 μ) ηx[
                    t] + ((9 ηy[t])/4)};

perSol = {ηx[t], ηy[
                        t]} /. 
DSolve[{eqPer, ηx[
                              0] == 0.0001, ηy[0] == 
                            0.0001,  ηx'[0] == 
                            0, ηy'[0] == 0} // 
                        Flatten, {ηx[t], ηy[t]}, 
                      t
] /. μ -> 0.01215  // 
                Flatten // ExpandAll // ExpToTrig  // Simplify // 
        Chop

zeropt = {x0 -> (1/2) (1 - 2 μ), 
          y0 -> ((Sqrt[3])/2)} /. μ ->  .01215

position  = {x0, y0} + perSol /. zeropt

pt1 = 
        ParametricPlot[ position // Evaluate, {t, 0, 80} ];

nSol=PMotion[ 
        0.01215 ,  {0.48785 + 0.0001, ((Sqrt[3])/2) + 0.0001, 0 , 0 } // 
          Flatten, 80 , 1000]

fx[t_] := x[t] /. nSol ; 
fy[t_] := y[t] /. nSol ;

Mgraph[fx, fy, t, 
        80 , 
PlotPoints -> 100, 
        PlotStyle -> {Thickness[0.008], Dashing[{0.05}], Hue[0.8] }, 
        DisplayFunction -> $DisplayFunction, 
Epilog -> pt1[[1]]  ];

Problem 10:  Orbit  around  the  Earth and Moon

Solution

Clear["Global`*"]

sol = 
        PMotion[0.01215 , {1.2 , 0 , 0 , -1.5    }, 2  π , 
          3000] ;

fx[t_] := x[t] /. sol ; 
fy[t_] := y[t] /. sol ;

pt1 = 
        Mgraph[fx, fy, t, 
          2 π, 
PlotPoints -> 
            100, 
PlotStyle -> 
            Hue[0.8], 
DisplayFunction -> $DisplayFunction ];

ParametricPlot[
{{  Cos[t],   Sin[
                t] }, 
{fx[t] Cos[t] - 
                fy[t] Sin[t],  
fx[t] Sin[t] + 
                fy[t] Cos[t]}} // Evaluate, 
{t, 
          0,    2 π}, 
PlotStyle -> {GrayLevel[
              0.3], Hue[0.6], 
            Hue[0.9]}, 
AspectRatio -> 
          Automatic];

7.4 Exercises

Exercise 1: Central force problems

Exercise 2: Central force procedure

Exercise 3: Central forces and elliptical solutions

Exercise 4 : Attractive Inverse fifth power force

Exercise 5: Eccentric Anomaly

Exercise 6: Eccentric Anomaly

Exercise 7: Kepler problem with drag

Exercise 8: Lagrange points

Exercise 9:  Orbit around the Sun and Jupiter

Exercise 10: Exact Solution for the Three-Body Problem