Npar = 10000 x = rnorm(Npar,0, 5) xp = rnorm(Npar,0,2) sqrt(var( x)) sqrt(var(xp)) hist(x) hist(xp) plot(x,xp,pch=".") CScalc = function(x,xp){ epsx=sqrt(var(x)*var(xp)-cov(x,xp)^2) betx= var(x )/epsx alfx=-cov(x,xp)/epsx gamx= var(xp)/epsx c(alfx,betx,gamx,epsx) } CScalc(x,xp) cov(x,xp) # drift a distance L L = 3 x = x + xp*L hist(x) hist(xp) plot(x,xp,pch=".") CScalc(x,xp) # focus through a quadrupole, focal length F F = 2 xp = xp-x/F hist(x) hist(xp) plot(x,xp,pch=".") CScalc(x,xp) L = c( 5, 2, 7, 4, 3, 3, 6, 8) F = c(-8, 4,-9, 4,-7, 9,-5, 7) i = 0 V = CScalc(x,xp) # values at s=0 while(i < length(L)){ i = i+1 xp = xp-x/F[i] x = x+xp*L[i] V = rbind(V,CScalc(x,xp)) } V alpha = V[,1] beta = V[,2] gamma = V[,3] eps = V[,4] s = c(0,cumsum(L)) plot(s, beta ,typ="l",ylim=c(0,1.2*max( beta ))) plot(s,sqrt(beta),typ="l",ylim=c(0,1.2*max(sqrt(beta))))