FODO Cell

Transport matrix for cell

Hlaf length focussing quad:

mhalff = {{1, 0}, {-(1/(2 f)), 1}}

Full length focussing quad:

mf = mhalff . mhalff

Full length defocussing quad:

md = mf /. f -> -f

Drift space:

mdrift = {{1, L}, {0, 1}}

Total cell:

mcell = Simplify[mhalff . mdrift . md . mdrift . mhalff]

Various functions at center of focussing quad.  Because of symmetry, alpha is zero.  Therefore, the Cos[phi] is:

csphi = mcell[[1, 1]]

The beta function is:

betaplus = 
      Simplify[mcell[[1, 2
            ]]/(Sqrt[1 - csphi csphi])]

Alpha is:

alphaplus = 0.

And gamma is:

gammaplus = (1/betaplus)

Transport of the Beta function along the beam line (as derived from the beam matrix):

betanew[mm_, beta_, alpha_, gamma_] = 
      mm[[1, 1]]^2 beta - 
        2 mm[[1, 1]] 
          mm[[1, 2]] alpha + 
        mm[[1, 2]]^2 gamma

Check to see that the beta function is the same after drifting one cell:

Simplify[betanew[mcell, betaplus, alphaplus, gammaplus]]

alphanew[mm_, beta_, alpha_, gamma_] = 
      -mm[[1, 1]] 
          mm[[2, 1]] beta + 
        (mm[[2, 1]] 
                mm[[1, 2]] + 
              mm[[1, 1]] 
                mm[[2, 2]]) alpha - 
        mm[[1, 2]] 
          mm[[2, 2]] gamma

Simplify[alphanew[mcell, betaplus, alphaplus, gammaplus]]

gammanew[mm_, beta_, alpha_, gamma_] = 
      mm[[2, 1]]^2 beta - 
        2 mm[[2, 1]] 
          mm[[2, 2]] alpha + 
        mm[[2, 2]]^2 gamma

Simplify[gammanew[mcell, betaplus, alphaplus, gammaplus]]

Single particle plots

emmitance = ((deltax^2)/betaplus)

cond = {deltax -> 1. , L -> 1. }

x[phi_, beta_, e_] = (Sqrt[e beta]) Cos[phi]

xp[phi_, beta_, alpha_, e_] = 
      -(Sqrt[(e/beta)]) (alpha Cos[phi] + Sin[phi])

plotplus = 
      ParametricPlot[{x[phi, betaplus, emmitance], 
            xp[phi, betaplus, alphaplus, emmitance]} /. 
          Flatten[{cond, f -> 1. }], {phi, 0, 2 π}]

plotmf = 
      ParametricPlot[{
            x[phi, betanew[mhalff, betaplus, alphaplus, gammaplus], 
              emmitance], 
            xp[phi, betanew[mhalff, betaplus, alphaplus, gammaplus], 
              alphanew[mhalff, betaplus, alphaplus, gammaplus], 
              emmitance]} /. Flatten[{cond, f -> 1. }], {phi, 0, 
          2 π}]

plothdrift1 = 
      ParametricPlot[{
            x[phi, betanew[(mdrift /. L -> (L/2)) . mhalff, betaplus, 
                alphaplus, gammaplus], emmitance], 
            xp[phi, betanew[(mdrift /. L -> (L/2)) . mhalff, 
                betaplus, alphaplus, gammaplus], 
              alphanew[(mdrift /. L -> (L/2)) . mhalff, betaplus, 
                alphaplus, gammaplus], emmitance]} /. 
          Flatten[{cond, f -> 1. }], {phi, 0, 2 π}]

plotdrift1 = 
      ParametricPlot[{
            x[phi, betanew[mdrift . mhalff, betaplus, alphaplus, gammaplus], 
              emmitance], 
            xp[phi, betanew[mdrift . mhalff, betaplus, alphaplus, gammaplus], 
              alphanew[mdrift . mhalff, betaplus, alphaplus, gammaplus], 
              emmitance]} /. Flatten[{cond, f -> 1. }], {phi, 0, 
          2 π}]

plotminus = 
      ParametricPlot[{
            x[phi, betanew[(md /. f -> 2 f) . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus], emmitance], 
            xp[phi, betanew[(md /. f -> 2 f) . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus], 
              alphanew[(md /. f -> 2 f) . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus], emmitance]} /. 
          Flatten[{cond, f -> 1. }], {phi, 0, 2 π}]

plotmd = 
      ParametricPlot[{
            x[phi, betanew[md . mdrift . mhalff, betaplus, alphaplus, 
                gammaplus], emmitance], 
            xp[phi, betanew[md . mdrift . mhalff, betaplus, alphaplus, 
                gammaplus], 
              alphanew[md . mdrift . mhalff, betaplus, alphaplus, gammaplus], 
              emmitance]} /. Flatten[{cond, f -> 1. }], {phi, 0, 
          2 π}]

plothdrift2 = 
      ParametricPlot[{
            x[phi, betanew[
                (mdrift /. L -> (L/2)) . md . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus], emmitance], 
            xp[phi, betanew[
                (mdrift /. L -> (L/2)) . md . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus], 
              alphanew[(mdrift /. L -> (L/2)) . md . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus], emmitance]} /. 
          Flatten[{cond, f -> 1. }], {phi, 0, 2 π}]

plotdrift2 = 
      ParametricPlot[{
            x[phi, betanew[mdrift . md . mdrift . mhalff, betaplus, 
                alphaplus, gammaplus], emmitance], 
            xp[phi, betanew[mdrift . md . mdrift . mhalff, betaplus, 
                alphaplus, gammaplus], 
              alphanew[mdrift . md . mdrift . mhalff, betaplus, alphaplus, 
                gammaplus], emmitance]} /. Flatten[{cond, f -> 1. }], {
          phi, 0, 2 π}]

Show[plotplus, plotmf, plothdrift1, plotdrift1, plotminus, plotmd, 
      plothdrift2, plotdrift2]

Peak in the emmitance

For stable operation L/2f<1.  For a given aperature in the focussing quads (deltax), the emmitance is maximized at some point:

Plot[emmitance /. cond, {f,  .5, 10}]

Find this point:

fmax = 
      Simplify[Solve[
          (∂_f)(emmitance /. {L -> 1, deltax -> 1}) == 
            0, f]]

N[fmax]

Plots of beta functions vs. L

fmax:

plot1 = 
      Plot[If[l ≤ 1, 
            betanew[(mdrift /. L -> l) . mhalff, betaplus, 
                alphaplus, gammaplus] /. cond, 
            betanew[(mdrift /. L -> l - 1) . md . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus] /. cond, Null] /. 
          fmax[[2]], {l, 0, 2}, 
        DisplayFunction -> Identity]

plot1a = 
      Plot[If[l ≤ 1, 
            betanew[(mdrift /. L -> l) . mhalff, betaplus, 
                alphaplus, gammaplus] /. cond, 
            betanew[(mdrift /. L -> l - 1) . md . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus] /. cond, Null] /. 
          f -> 
            -(f /. fmax[[2]]), {
          l, 0, 2}, DisplayFunction -> Identity]

Show[plot1, plot1a, DisplayFunction -> $DisplayFunction]

long focal length example:

plot2 = 
      Plot[If[l ≤ 1, 
            betanew[(mdrift /. L -> l) . mhalff, betaplus, 
                alphaplus, gammaplus] /. cond, 
            betanew[(mdrift /. L -> l - 1) . md . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus] /. cond, Null] /. 
          f -> 5, {l, 0, 2}, DisplayFunction -> Identity]

plot2a = 
      Plot[If[l ≤ 1, 
            betanew[(mdrift /. L -> l) . mhalff, betaplus, 
                alphaplus, gammaplus] /. cond, 
            betanew[(mdrift /. L -> l - 1) . md . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus] /. cond, Null] /. 
          f -> -5, {l, 0, 2}, DisplayFunction -> Identity]

Show[plot2, plot2a, DisplayFunction -> $DisplayFunction]

And a very short focal length:

plot3 = 
      Plot[If[l ≤ 1, 
            betanew[(mdrift /. L -> l) . mhalff, betaplus, 
                alphaplus, gammaplus] /. cond, 
            betanew[(mdrift /. L -> l - 1) . md . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus] /. cond, Null] /. 
          f ->  .55, {l, 0, 2}, DisplayFunction -> Identity]

plot3a = 
      Plot[If[l ≤ 1, 
            betanew[(mdrift /. L -> l) . mhalff, betaplus, 
                alphaplus, gammaplus] /. cond, 
            betanew[(mdrift /. L -> l - 1) . md . mdrift . mhalff, 
                betaplus, alphaplus, gammaplus] /. cond, Null] /. 
          f -> - .55, {l, 0, 2}, DisplayFunction -> Identity]

Show[plot3, plot3a, DisplayFunction -> $DisplayFunction]

Show[plot1, plot2, plot3, DisplayFunction -> $DisplayFunction]