# Program to plot particle phase space motion
#  in presence of a sextupole field

nu  <- 0.47001014 # zero-amplitude tune
Nturns = 200
amax = 1

# Rotation matrix elements for linear propagation:
a <-  cos(2*pi*nu)
b <-  sin(2*pi*nu)
c <- -b
d <-  a

# initialize the phase space plot
u    = 0
du   = 0
plot(u, du, xlim=c(-amax,amax), ylim=c(-amax,amax), typ="n")
text(0.75*amax,amax,paste("nu = ",round(nu,3)))

trk = 1
while (trk < 15) {
# initialize particle positions in phase space
  u0 <- locator(1)
  u  <- u0$x
  du <- u0$y
# track the particle...
  i = 1
  while (i < Nturns+1) {
    du1  <-  du  + u*u/2
	 u1   <- a*u  + b*du1
	 du   <- c*u  + d*du1 + u1*u1/2
    u    <- u1
    points(u, du, pch=".")
    i = i + 1
  }
trk = trk + 1
}