# Simple R script to visualize longitudinal dynamics and phase stability
# in an accelerator

# number of particles to simulate
Np = 400
# generate initial energy deviations [eV]
dEmin = -1e8; dEmax = 1e8
dE0 = runif(Np,dEmin,dEmax)
# generate initial phases [rad]
dphimin = 0; dphimax = 6*pi
dphi0 = runif(Np,dphimin,dphimax)

# speed of light [m/s]
clight = 299792458
# rf frequency [Hz]
frf = 53e6
# rf voltage kick [eV]
qV = 1e6
# synchronous phase [rad]
phis = 135*pi/180
# phase slip factor
eta = 0.02
# synchronous energy [eV]
Es = 8.938272e9
# kinematic quantities
mp = 938.272e6 # proton mass [eV]
pc = sqrt(Es^2 - mp^2)
beta = pc/Es
# machine circumference [m]
L = 484
# revolution time [s]
tau = L / beta / clight

# synchrotron frequency
fs = sqrt(-eta*2*pi*frf*qV*cos(phis)/tau/beta^2/Es) /2/pi
Ns = 1/fs/tau

# print out main parameters
txtmsg = sprintf('RF frequency = %s MHz\nRF voltage = %s MV\nSynchronous phase = %s deg\nSlip factor = %s\nSynchronous energy = %s GeV\nRevolution time = %s us\nSynchrotron frequency = %s kHz\nSynchrotron period = %s turns',
                 formatC(frf/1e6, di=3),
                 formatC(qV/1e6, di=3),
                 formatC(phis/pi*180, di=3),
                 formatC(eta, di=3),
                 formatC(Es/1e9, di=3),
                 formatC(tau*1e6, di=3),
                 formatC(fs/1e3, di=3),
                 formatC(Ns))
message(txtmsg)

# number of turns
N = 181

# colors
library(RColorBrewer)
matlab.colors =colorRampPalette(
  c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
    "#FF7F00", "red", "#7F0000"))
my.seq.pal = colorRampPalette(brewer.pal(9, "YlGnBu"))
# color represents time
co = rev(my.seq.pal(N))

# initial conditions
dE = dE0
dphi = dphi0

# prepare canvas
wi = 8; he = 6; pt = 10; fo = "Times"; fn = "tmp"
dev.new(width=wi, height=he)
par(oma=rep(0, 4), mar=c(4, 4, 0.5, 0.5), mgp=c(2, 0.5, 0), tcl=-0.2,
    las=1, bg="black",
    fg="white", col="white", col.axis="white", col.lab="white",
    col.main="white", col.sub="white")

plot(range(dphi0*180/pi), range(dE0)/1e6, type="n",
  xlab="phase [deg]", ylab="energy deviation [MeV]", axes=FALSE)
axis(1, at=seq(0, by=90, length.out=20))
axis(2)

# calculate turn-by-turn map
# of longitudinal dynamics
for (i in 1:N) {
  points(dphi*180/pi, dE/1e6, pch=".", col=co[i])
  dE = dE + qV*(sin(dphi)-sin(phis))
  dphi = dphi + (2*pi*frf*tau*eta/beta^2/Es)*dE
  Sys.sleep(0.05)
}

legend('topleft', legend=strsplit(txtmsg, '\n')[[1]], bg='black', cex=0.6)

box()

# save graphics to PDF file
dev.copy2pdf(file=paste(fn, ".pdf", sep=""), pointsize=pt, family=fo)