Synchrotron Stability -- Fermilab Booster

Difference equation for phase:

Clear[harmonic,slippage,beta,Es,deltaE]
deltaphi[harmonic_,slippage_,beta_,Es_,deltaE_] =
	harmonic*slippage/beta^2/Es * deltaE

Difference equation for energy:

deltadeltaE[Volt_,phase_,phases_] =
	Volt*(Sin[phase] - Sin[phases])

Booster parameters

harmonic = 84.;
voltage = .86;
synchronousE = 8000.;
phases = Pi;
massp = 938.27231;
beta = Sqrt[synchronousE^2-massp^2]/synchronousE
gamma = 1./Sqrt[1-beta^2]
slippage = 1/5.45^2-1/gamma^2

Subroutine to create table of phases and energies:

phasespace[phi_,ee_] := Module[{u,v,umax,vmax},
	umax=10;
	vmax=50;
	v= Table[{phi,ee},{vmax}];
	u = Table[{phi,ee},{umax}];
	Do[
	Do[u[[n]] = If[n == 1, 
		{v[[m-1,1]]+deltaphi[harmonic,slippage,
		beta,synchronousE,v[[m-1,2]]],
		v[[m-1,2]]+deltadeltaE[voltage,v[[m-1,1]],phases]},
		{u[[n-1,1]]+deltaphi[harmonic,slippage,
		beta,synchronousE,u[[n-1,2]]],
		u[[n-1,2]]+deltadeltaE[voltage,u[[n-1,1]],phases]}],
		{n,1,umax}];
		v[[m]]=u[[umax]];,{m,2,vmax}];
	v]

Tables of different trajectories for different energies:

plot25=phasespace[N[Pi],25.];

plot50=phasespace[N[Pi],50.];

plot75=phasespace[N[Pi],75.];

plot100=phasespace[N[Pi],100.];

plot125=phasespace[N[Pi],125.];

plot150=phasespace[N[Pi],150.];

Make plots from tables:

p25 = ListPlot[plot25,
	PlotJoined->True,
	PlotRange->{{0,3 Pi},{-150,150}},
	DisplayFunction->Identity];

p50 = ListPlot[plot50,
	PlotJoined->True,
	PlotRange->{{0,3 Pi},{-150,150}},
	DisplayFunction->Identity];

p75 = ListPlot[plot75,
	PlotJoined->True,
	PlotRange->{{0,3 Pi},{-150,150}},
	DisplayFunction->Identity];

p100 = ListPlot[plot100,
	PlotJoined->True,
	PlotRange->{{0,3 Pi},{-150,150}},
	DisplayFunction->Identity];

p125 = ListPlot[plot125,
	PlotJoined->True,
	PlotRange->{{0,3 Pi},{-150,150}},
	DisplayFunction->Identity];

p150 = ListPlot[plot150,
	PlotJoined->True,
	PlotRange->{{0,3 Pi},{-150,150}},
	DisplayFunction->Identity];

Show plots:

Show[p25,p50,p75,p100,p125,p150,
	DisplayFunction->$DisplayFunction,
	AxesLabel->{"Phase (rad)","Energy (MeV)"}]