# Program to plot longitudinal phase space motion
#  modelled on MI parameters

Nturns = 200

# Some Parameters 
Ws = 8000 # MeV
phis = 30*pi/180  # synchronous phase angle
eV = 1 # MeV
gamma = (938+Ws)/938
beta = sqrt(1-1/gamma^2)
eta = -1/gamma^2
h = 588
k = 2*pi*h*eta/beta^2*(gamma-1)/gamma/Ws

# initialize the phase space plot
phi   = 0
dW    = 0
plot(phi, dW, xlim=c(-180,180), ylim=c(-1,1)*50, typ="n")

trk = 1
while (trk < 30) {
# initialize particle positions in phase space
  u0 <- locator(1)
  phi   <- u0$x/180*pi
  dW    <- u0$y
# track the particle...
  i = 1
  while (i < Nturns+1) {
    phi = phi + k*dW
    dW  = dW + eV*(sin(phi)-sin(phis))
    points(phi*360/2/pi, dW, pch=21,col="red")
    i = i + 1
  }
trk = trk + 1
}