trajectory.mws

>    restart;

Define the central potential we're looking at:
 U := r -> -4/r+1/r^2;

U := proc (r) options operator, arrow; -4/r+1/(r^2) end proc

Values of the masses of the two bodies:
 m1 := 3; m2 := 6;

m1 := 3

m2 := 6

Find total mass and relative mass:

M := M1 + m2; mu := m1*m2/(m1+m2);

M := M1+6

mu := 2

Initial conditions: choose z axis such that r(t=0) is parallel to it -> theta0=0; Let alpha be the angle between v0 and r0.

>    r0 := 2; v0 := 1; alpha := Pi/3;

r0 := 2

v0 := 1

alpha := 1/3*Pi

compute angular momentum and the relative energy:
 l := evalf(mu*r0*v0*sin(alpha));
Erel := evalf(mu*v0^2/2 + U(r0));

l := 3.464101616

Erel := -.7500000000

define the effective potential and plot it:
 Ueff := r -> l^2/(2*mu*r^2)+U(r);
plot({Erel, Ueff(r)},r=0..20, y=-2*Erel .. 2*Erel);

Ueff := proc (r) options operator, arrow; 1/2*l^2/mu/r^2+U(r) end proc

[Maple Plot]

Find the turning points:

>    rmin := fsolve(Erel=Ueff(r), r, 0..r0);
rmax := fsolve(Erel=Ueff(r), r, r0..infinity);

rmin := 1.333333334

rmax := 3.999999999

We now want to let r go from r0->rmax->rmin->rmax-> ... and see what the trajectory is.  The function f(x) I'm defining below will do precisely this:

>    phase := arcsin((2*r0-rmin-rmax)/(rmax-rmin));
f := x -> (rmin+rmax)/2 +(rmax-rmin)/2*sin(2*Pi*x+phase);

phase := -.5235987758

f := proc (x) options operator, arrow; 1/2*rmin+1/2*rmax+1/2*(rmax-rmin)*sin(2*Pi*x+phase) end proc

>    plot({rmin, rmax, f(x)},x=0..10, 'r' =rmin/2.. 1.5*rmax);

[Maple Plot]

How many points should we plot? -> value N
What should be the step Dx between consecutive steps?
The solution will give you N*Dx full runs from rmin->rmax->rmin (the first starts at r0);

>    N := 250; Dx:=0.02;

N := 250

Dx := .2e-1

>    theta0 := 0; r[0] := r0; zp[0] := r0; xp[0] := 0; theta[0] := 0:

theta0 := 0

r[0] := 2

zp[0] := 2

xp[0] := 0

>    for i from 1 to N do
   x := (i-1)*Dx:
   r[i] := evalf(f(x)):
   Dtheta := evalf(int(l/r^2/sqrt(2*mu*(Erel-Ueff(r))),r=r[i-1]..r[i])):
   theta[i] := theta[i-1] + abs(Dtheta);
   zp[i] := r[i]*cos(theta[i]); xp[i] := r[i]*sin(theta[i]):
end do:

>    sol := [seq([xp[i],zp[i]], i=1..N)]:

>    plot(sol,labels=['x','z']);

[Maple Plot]