Planetary Motion in the Sitnikov System

Equations of motion

System of differential equations 
for motion of a planet along the axis 
of a binary star with small eccentricity ecc and initial phase phi0.
Choose units for which 2 G M = 1, average radius of orbits of stars is 1

equationz=z'[t]==zdot[t]
equationzdot=zdot'[t]==-z[t]/(  (1 + ecc*Cos[t+phi0])^2+z[t]^2  )^(3/2)

Time parameters:
binary star year
number of years to evolve the system,
number of positions to save during each year

year=2 Pi
Nyears = 10
Nsave=10

Binary star with circular orbit

Values of parameters
Initial conditions for planet in the plane of the binary star (zstart = 0) with specified velocity zdotstart

ecc=0.;
phi0=0.;
zstart= 0.
zdotstart =1.

Solve the system of differential equations numerically one year at a time for a total of  Nyear years,
saving the position z(t) of the planet Nsave times during each year.
Keep track of  each time the planet crosses the orbital plane of the binary star.

positionlist={{0.,0.}};
crossings={0.};
t0=0.;
Do[
initialz=z[t0]==zstart;
initialzdot=zdot[t0]==zdotstart;
solution=NDSolve[{equationz,equationzdot,initialz,initialzdot},{z[t],zdot[t]},{t,t0,t0+year}];
position[t_]=z[t]/.solution[[1]];
velocity[t_]=zdot[t]/.solution[[1]];
told=t0;
zold=ztart;
Do[tnew=told+(year/Nsave);
znew=position[tnew];
positionlist=Append[positionlist,{(tnew/year),znew}];If[zold znew<0,crossings=Append[crossings,((znew told-zold tnew)/((znew-zold)*year))]];
told=tnew;
zold=znew,
{m,1,Nsave}];
zstart=position[t0+year];
zdotstart=velocity[t0+year];
t0=t0+year,
{n,1,Nyears}]

Plot the position of the planet as a function of time

ListPlot[positionlist,PlotRange->Automatic,AxesLabel->{"years","z"},Joined->True]

List the times when the planet crosses the plane of the binary star,
the time interval between crossings, 
and the maximum time interval.

crossings
crossingintervals=Delete[crossings,1]-Delete[crossings,-1]
maxinterval=Max[crossingintervals]

Find the escape velocity (zdotstart) for ecc = 0 to 3 digits accuracy.
Plot the position of the planet as a function of time for a velocity above and a velocity below the escape velocity.
Show your results to the instructor.

Binary star with elliptical orbit with small eccentricity

Change the eccentricity to ecc = 0.1.
Set the initial velocity equal to the escape velocity for ecc = 0.
Plot z(t) as a function of t for 6 values of the initial phase phi0 ranging from 0 to 2 Pi.

Vary the initial phase phi0 and the initial velocity zdotstart, 
looking for solutions for which the time interval between the planet crossing the orbit of the binary star is as large as possible 
and for solutions for which the final crossing of the planet before escaping is as large as possible.  
Report your record values to  the instructor.
Note the extreme sensitivity to initial conditions.