Construction of a Simple Collider

The parameters used here are different than in those for your take home exam.
The program is very rudimentary and was only written to help you with your final. If you find errors, please report them to me, this will help future students.

Parameters

Numbers that distinguish different accelerator elements:

Off[General::"spell1"]; 
    Off[General::"spell"];
eleD = 1;
eleQ = 2;
eleB = 3;
eleS = 4;
On[General::"spell1"]; On[General::"spell"];

Formulas

The transport matrices of some elements:

mD[len_, dim_] := 
      Module[{i, m = Table[0, {2  dim}, {2  dim}]}, 
        Do[m[[2  i - 1, 2  i - 1]] = 1; m[[2  i - 1, 2  i]] = len; 
          m[[2  i, 2  i]] = 1, {i, dim}]; m]
mQ[kl_, dim_] := 
      Module[{i, k = {-kl, kl, 0}, m = Table[0, {2  dim}, {2  dim}]}, 
        Do[m[[2  i - 1, 2  i - 1]] = 1; 
          m[[2  i, 2  i - 1]] = k[[i]]; 
          m[[2  i, 2  i]] = 1, {i, dim}]; m]
mB[ϕ_, dim_] := 
      Module[{i, m = Table[0, {2  dim}, {2  dim}]}, 
        Do[m[[i, i]] = 1, {i, 2  dim}]; 
        If[dim == 3, m[[2, 6]] = ϕ; 
          m[[5, 1]] = -ϕ]; m]
mS[k2l_, dim_] := 
      Module[{i, m = Table[0, {2  dim}, {2  dim}]}, 
        Do[m[[i, i]] = 1, {i, 2  dim}]; m]
mD[len_] := mD[len, 2]
mQ[kl_] := mQ[kl, 2]

Computation of Twiss parameters from a matrix for 1,2, or 3 degrees of freedom.

twiss2[m_] := 
      Module[{μ, sμ, β, γ, α}, μ = 
          ArcCos[(1/2) Tr[m]]; 
        If[m[[1, 2]] < 0, μ = 2  π - μ]; 
        sμ = Sin[μ]; β = 
          m[[1, 2]]/sμ; γ = -m[[2, 1]]/
            sμ; α = (1/2) (m[[1, 1]] - m[[2, 2]])/
              sμ; {α, β, γ, (μ/(2  π))}]
twiss4[m_] := {twiss2[Table[m[[i, j]], {i, 1, 2}, {j, 1, 2}]], 
        twiss2[Table[
            m[[i, j]], {i, 3, 4}, {j, 3, 4}]]}
twiss6[m_] := 
      Module[{i, j, η}, η = 
          Inverse[{{1, 0}, {0, 1}} - 
                Table[m[[i, j]], {i, 2}, {j, 2}]] . {m[[1, 6]], 
              m[[2, 6]]}; {Join[
            twiss2[Table[m[[i, j]], {i, 1, 2}, {j, 1, 2}]], η], 
          Join[twiss2[Table[m[[i, j]], {i, 3, 4}, {j, 3, 4}]], {0, 
              0}]}]

Propargation of Twiss parameters trough an element that is described by its length and its transport matrix:

Off[General::"spell1", 
      General::"spell"]
prop[ele_, dim_] := 
      Module[{m, i, j, len = ele[[1]], mβ0, 
          mβ, Δμx, Δμy, η}, 
        m = Table[
            ele[[2, i, j]], {i, 4}, {j, 
              4}]; 
mβ0 = {{βxlat[[-1, 
                  2]], -αxlat[[-1, 2]], 0, 
              0}, {-αxlat[[-1, 
                    2]], γxlat[[-1, 2]], 0, 0}, {0, 
              0, βylat[[-1, 
                  2]], -αylat[[-1, 2]]}, {0, 
              0, -αylat[[-1, 2]], γylat[[-1, 
                  2]]}}; 
        mβ = 
          m . mβ0 . Transpose[m]; 
βxlat = 
          Append[βxlat, {βxlat[[-1, 1]] + len, 
              mβ[[1, 1]]}]; αxlat = 
          Append[αxlat, {αxlat[[-1, 1]] + 
                len, -mβ[[1, 2]]}]; γxlat = 
          Append[γxlat, {γxlat[[-1, 1]] + len, 
              mβ[[2, 2]]}]; βylat = 
          Append[βylat, {βylat[[-1, 1]] + len, 
              mβ[[3, 3]]}]; αylat = 
          Append[αylat, {αylat[[-1, 1]] + 
                len, -mβ[[3, 4]]}]; γylat = 
          Append[γylat, {γylat[[-1, 1]] + len, 
              mβ[[4, 
                  4]]}]; 
Δμx = 
          ArcTan[m[[1, 
                    1]] (Sqrt[βxlat[[-2, 2]]/βxlat[[-1, 2]]]) - αxlat[[-2, 2]] 
                m[[1, 
                      2]]/(Sqrt[βxlat[[-1, 2]] βxlat[[-2, 2]]]), 
            m[[1, 2]]/(Sqrt[βxlat[[-1, 2]] βxlat[[-2, 2]]])]; 
Δμy = 
          ArcTan[m[[3, 
                    3]] (Sqrt[βylat[[-2, 2]]/βylat[[-1, 2]]]) - αylat[[-2, 2]] 
                m[[3, 
                      4]]/(Sqrt[βylat[[-1, 2]] βylat[[-2, 2]]]), 
            m[[3, 4]]/(Sqrt[βylat[[-1, 2]] βylat[[-2, 2]]])]; 
ψxlat = 
          Append[ψxlat, {ψxlat[[-1, 1]] + 
                len, ψxlat[[-1, 
                    2]] + Δμ(x/(2  π))}]; 
ψylat = 
          Append[ψylat, {ψylat[[-1, 1]] + 
                len, ψylat[[-1, 
                    2]] + Δμ(y/(2  π))}]; 
If[
          dim > 2, 
η = 
            ele[[2]] . {ηxlat[[-1, 
                    2]], ηpxlat[[-1, 
                    2]], ηylat[[-1, 
                    2]], ηpylat[[-1, 2]], 0, 
                1}; 
ηxlat = 
            Append[ηxlat, {ηxlat[[-1, 1]] + 
                  len, η[[1]]}]; ηpxlat = 
            Append[ηpxlat, {ηpxlat[[-1, 1]] + 
                  len, η[[2]]}]; 
ηylat = 
            Append[ηylat, {ηylat[[-1, 1]] + 
                  len, η[[3]]}]; ηpylat = 
            Append[ηpylat, {ηpylat[[-1, 1]] + 
                  len, η[[4]]}]
]
]
On[General::"spell1", General::"spell"]

Propagation of Twiss parameters through a number of optical elements:

proper[lat_, dim_] := 
      Do[prop[lat[[j]], dim], {j, 1, Length[lat]}]

Compute a series of matrices from a lattice description:

matform[lat_, dim_] := 
      Map[lenmat[#, dim] &, lat]
lenmat[ele_, dim_] := {ele[[1]], 
        If[ele[[2]] == eleD, mD[ele[[3]], dim], 
          If[ele[[2]] == eleQ, mQ[ele[[3]], dim], 
            If[ele[[2]] == eleB, mB[ele[[3]], dim], 
              If[ele[[2]] == eleS, 
                mS[ele[[3]], dim]]]]]}

Compute a series of matrices and second order kicks from a lattice description:

kickform[lat_, dim_] := 
      Map[lenkick[#, dim] &, lat]
lenkick[ele_, dim_] := 
      Module[{res}, 
        res = {ele[[1]], 
            If[ele[[2]] == eleD, mD[ele[[3]], dim], 
              If[ele[[2]] == eleQ, mQ[ele[[3]], dim], 
                If[ele[[2]] == eleB, mB[ele[[3]], dim], 
                  If[ele[[2]] == eleS, mS[ele[[3]], dim]]]]]}; 
        If[ele[[2]] == eleS, 
          res = Append[res, 
              If[dim == 1, {ele[[3]]}, 
                If[dim ≥ 
                    2, {{{-ele[[3]] (1/2), 0}, {0, 
                        ele[[3]] (1/2)}}, {{0, -(1/2) 
                          ele[[3]]}, {-(1/2) ele[[3]], 
                        0}}}]]]]; res]

getm[lat_] := 
      Module[{m = lat[[1, 2]]}, 
        Do[m = lat[[i, 2]] . m, {i, 2, Length[lat]}]; m]

ot[n] plots a resonance diagram up to order n:

Off[General::"spell1"]
range = {{0, 1}, {0, 1}};
$TextStyle  =  {FontFamily  ->  "Times",  
          FontSize  ->  15};
lin[n_, m_] := {Line[{{0, n}, {1, n - m}}]}
lint[n_, m_] := {Line[{{n, 0}, {n - m, 1}}]}
\(ot[ord_] := 
      Module[{g = {}}, 

        Do[Do[Do[
              If[m < 0, g = Append[g, RGBColor[0, 1, 0]], 
                g = Append[g, RGBColor[0, 0, 0]]]; 
              g = Append[g, lin[n\/l, m\/l]], {n, Min[m, 0], 
                Max[l, m + l]}], {m, \(-ord\) + l, ord - l}]; 
          Do[g = Append[g, Line[{{n\/l, 0}, {n\/l, 1}}]], {n, 1, l}]; 
          Do[g = Append[g, Line[{{0, n\/l}, {1, n\/l}}]], {n, 1, l}], {l, 1, 
            ord}]; 
Show[Graphics[g], PlotRange -> range, 
          ImageSize -> {400, 400}, AspectRatio -> 1, 
          Frame -> True]]\)
On[General::"spell1"]

The fit routine searches for an N dimensional set of parameters lets the N dimensional function fct have the values of the N dimensional goal vector. The search will start with the initial parameters para and no more than nmax iterations will be performed .

Off[General::"spell1"]
\(fit[fct_, goal_, para_, nmax_] := 
      Module[{resi, Δp, mα, p, res, i, j, pi = para, 
          dim = Length[para], sfac = 0.5}, 

        resi = fct[pi]; 
Δp = 
          Table[If[pi[\([\)\(i\)\(]\)] == 0, 0.01, 
              pi[\([\)\(i\)\(]\)]*0.01], {i, dim}]; 

        Do[
mα = 
            Table[p = pi; 
              p[\([\)\(i\)\(]\)] = 
                pi[\([\)\(i\)\(]\)] + Δp[\([\)\(i\)\(]\)]; 
              res = fct[p]; 
              Table[\(res[\([\)\(j\)\(]\)] - resi[\([\)\(j\)\(]\)]\)\/Δp[\([\)\(i\)\(]\)], {j, dim}], {i, 
                dim}]; 
p = pi + Δp; 
          res = fct[p]; Δp = 
            sfac\ Inverse[Transpose[mα]] . \((goal - res)\); pi = p; 
          resi = res; 
If[\((res - goal)\) . \((res - 
                    goal)\) < 10\^\(-8\), Break[]], 
{i, 
            nmax}]; 
Print["Accuracy: ", 
          res - goal]; 
Print["Best set of parameters"]; p]\)
On[General::"spell1"]

Split elements in a series of optical elements in smaller parts for improved resolution:

Off[General::"spell1"]
split[lat_, n_] := 
      Module[{res = {}, i}, 
        Do[If[lat[[i, 1]] == 0, res = Append[res, lat[[i]]], 
            Do[res = 
                Append[res, {lat[[i, 1]]/n, lat[[i, 2]], 
                    lat[[i, 3]]/n}], {n}]], {i, Length[lat]}]; 
        res]
On[General::"spell1"]

Off[General::"spell1"]
eletrack[ele_, z_] := 
      Module[{dim = Length[ele]}, 
        mtrack[ele[[2]], 
          If[dim == 2, z, 
            If[dim == 3, 
              z + {0, ele[[3, 1]] (((z[[1]]))^2)}]]]]
lattrack[lat_, z_] := 
      Module[{z0 = z}, Do[z0 = eletrack[lat[[i]], z0], {i, Length[lat]}]; 
        z0]
On[General::"spell1"]

FODO

Definition of a FODO:

fodo = {{0, eleQ, (kl/2)}, {len, eleD, len}, {0, eleQ, -kl}, {len, 
            eleD, len}, {0, eleQ, (kl/2)}};

kl = 1. ; len = 1. ;

α, β, γ, ψ/2π in x and y of a FODO and initiallization for the propagation of Twiss parameters:

twiss0 = twiss4[getm[matform[fodo, 2]]]
βxlat = {{0, twiss0[[1, 2]]}};
βylat = {{0, twiss0[[2, 2]]}};
αxlat = {{0, 
            twiss0[[1, 1]]}};
αylat = {{0, 
            twiss0[[2, 1]]}};
γxlat = {{0, 
            twiss0[[1, 3]]}};
γylat = {{0, 
            twiss0[[2, 3]]}};
ψxlat = {{0, 0}};
ψylat = {{0, 0}};

proper[matform[fodo, 2], 2]

βxlat

Show[
      Graphics[{RGBColor[1, 0, 0], Line[βxlat], RGBColor[0, 1, 0], 
          Line[βylat]}], Frame -> True, Axes -> False, 
      PlotRange -> {All, {0, 4}}]

Split each long element for higher resolution:

twiss0 = twiss4[getm[matform[fodo, 2]]]
βxlat = {{0, twiss0[[1, 2]]}}; 
βylat = {{0, twiss0[[2, 2]]}}; 
αxlat = {{0, twiss0[[1, 1]]}}; 
αylat = {{0, twiss0[[2, 1]]}}; 
γxlat = {{0, twiss0[[1, 3]]}}; 
γylat = {{0, twiss0[[2, 3]]}}; 
ψxlat = {{0, 0}}; 
ψylat = {{0, 0}};
proper[matform[split[fodo, 10], 2], 2]
Show[
      Graphics[{RGBColor[1, 0, 0], Line[βxlat], RGBColor[0, 1, 0], 
          Line[βylat]}], Frame -> True, Axes -> False, 
      PlotRange -> {All, {0, 4}}]

Arc of Half the Ring

Definition the arc:

Off[General::"spell1"]
nfodo = 20;
fodo = {{0, eleQ, (kl/2)}, {len, eleD, len}, {0, eleQ, -kl}, {len, 
            eleD, len}, {0, eleQ, (kl/2)}};
arc = Apply[Join, Table[fodo, {i, (nfodo/2)}]]; 
    On[General::"spell1"]

kl = 1. ; len = 1. ;

Initiallization Twiss parameters for the FODO, propagate through the arc and plot:

twiss0 = twiss4[getm[matform[fodo, 2]]];
βxlat = {{0, twiss0[[1, 2]]}}; 
βylat = {{0, twiss0[[2, 2]]}}; 
αxlat = {{0, twiss0[[1, 1]]}}; 
αylat = {{0, twiss0[[2, 1]]}}; 
γxlat = {{0, twiss0[[1, 3]]}}; 
γylat = {{0, twiss0[[2, 3]]}}; 
ψxlat = {{0, 0}}; 
ψylat = {{0, 0}}; proper[matform[arc, 2], 2]
Show[
      Graphics[{RGBColor[1, 0, 0], Line[βxlat], RGBColor[0, 1, 0], 
          Line[βylat]}], Frame -> True, Axes -> False, 
      PlotRange -> {All, {0, 4}}]

Interaction Region

Clear[l1, l2, k1, k2]
ir1[l1_, l2_, k1_, 
          k2_] = {{0, eleQ, k1}, {l1, eleD, l1}, {0, eleQ, k2}, {l2, eleD, 
            l2}};

irmat[l1_, l2_, k1_, k2_] := 
      Evaluate[matform[ir1[l1, l2, k1, k2], 2]]

fodo = {{0, eleQ, (kl/2)}, {len, eleD, len}, {0, eleQ, -kl}, {len, 
            eleD, len}, {0, eleQ, (kl/2)}};

kl = 1. ; len = 1. ;

Twiss paramtere propagation through the first half of the interaction region:

βip[p_] := 
      Module[{twiss0, l1 = p[[1]], l2 = p[[2]], k1 = p[[3]], 
          k2 = p[[4]]}, 
twiss0 = 
          twiss4[getm[
              matform[fodo, 2]]]; 
βxlat = {{0, 
              twiss0[[1, 2]]}}; 
βylat = {{0, 
              twiss0[[2, 2]]}}; 
αxlat = {{0, 
              twiss0[[1, 1]]}}; 
αylat = {{0, 
              twiss0[[2, 1]]}}; 
γxlat = {{0, 
              twiss0[[1, 3]]}}; 
γylat = {{0, 
              twiss0[[2, 3]]}}; 
ψxlat = {{0, 
              0}}; 
ψylat = {{0, 
              0}}; 
proper[irmat[l1, l2, k1, k2], 
          2]; 
If[pic == 1, 
          Show[Graphics[{RGBColor[1, 0, 0], Line[βxlat], 
                RGBColor[0, 1, 0], Line[βylat]}], Frame -> True, 
            Axes -> False, 
            PlotRange -> {All, All}]]; {βxlat[[-1, 
              2]], βylat[[-1, 2]], αxlat[[-1, 
              2]], αylat[[-1, 2]]}]

pic = 1; βstart = βip[{0.5, 0.2, 0.7, -0.7}]

Fitting parameters to obtain two specified final β and α functions. For that one starts with some realizable values and then performes many steps with mall changes of the goal until the final goal is reached.

pic = 0;
fit[βip, βstart, {0.5, 0.2, 0.7, -0.7}, 
      100]

pstart = 
      fit[βip, {1. , 2. , 0, 0}, {0.5, 0.2, 0.6, -0.7}, 
        100]

pstart = fit[βip, {1. , 1.6, 0, 0}, pstart, 100]

pstart = fit[βip, {1. , 1.2, 0, 0}, pstart, 100]

pstart = fit[βip, {1.0, 0.9, 0, 0}, pstart, 100]

pstart = fit[βip, {0.5, 0.2, 0, 0}, pstart, 100]

pstart = fit[βip, {0.3, 0.1, 0, 0}, pstart, 100]

pstart = fit[βip, {0.2, 0.05, 0, 0}, pstart, 100]

Print["l1 = ", l1  = pstart[[1]]]
Print["l2 = ", 
      l2  = pstart[[2]]]
Print["k1 = ", 
      k1  = pstart[[3]]]
Print["k2 = ", k2  = pstart[[4]]]

Lattice

nfodo = 20;
len = 1. ;
kl = 1. ;
l1 = 0.1391000706229664;
l2 = 0.3598037760308997;
k1 = 3.645324218395703;
k2 = -5.175733201180464;
fodo = {{0, eleQ, (kl/2)}, {len, eleD, len}, {0, eleQ, -kl}, {len, 
            eleD, len}, {0, eleQ, (kl/2)}};
arc = 
        Apply[Join, 
          Table[fodo, {i, (nfodo/2)}]];

ir1 = {{0, eleQ, k1}, {l1, eleD, l1}, {0, eleQ, k2}, {l2, eleD, 
            l2}};

ir2 = Reverse[ir1];

lattice = Join[arc, ir1, ir2, arc];

twiss0 = twiss4[getm[matform[fodo, 2]]];
βxlat = {{0, twiss0[[1, 2]]}};
βylat = {{0, twiss0[[2, 2]]}};
αxlat = {{0, twiss0[[1, 1]]}};
αylat = {{0, twiss0[[2, 1]]}};
γxlat = {{0, twiss0[[1, 3]]}};
γylat = {{0, twiss0[[2, 3]]}};
ψxlat = {{0, 0}};
ψylat = {{0, 0}}; proper[matform[lattice, 2], 2]
Show[
      Graphics[{RGBColor[1, 0, 0], Line[βxlat], RGBColor[0, 1, 0], 
          Line[βylat]}], Frame -> True, Axes -> False, 
      PlotRange -> {All, {0, 4}}, 
      ImageSize -> {1000, 500}]

Dispersion

Here the missing magnet scheme is use which brings the dispersion and its slope to 0 at the IP:

nfodo = 20;
len = 1. ;
kl = 1. ;
l1 = 0.1391000706229664;
l2 = 0.3598037760308997;
k1 = 3.645324218395703;
k2 = -5.175733201180464;
ϕ0 = ((2  π)/(nfodo - 2));
fodo = {{0, eleQ, (kl/2)}, {(len/2), eleD, (len/2)}, {0, 
            eleB, ϕ0}, {(len/2), eleD, (len/2)}, {0, 
            eleQ, -kl}, {(len/2), eleD, (len/2)}, {0, 
            eleB, ϕ0}, {(len/2), eleD, (len/2)}, {0, eleQ, 
            (kl/2)}};
fodo1 = {{0, eleQ, (kl/2)}, {(len/2), eleD, (len/2)}, {0, eleB, 
            0}, {(len/2), eleD, (len/2)}, {0, eleQ, -kl}, {(len/2), eleD, 
            (len/2)}, {0, eleB, 0}, {(len/2), eleD, (len/2)}, {0, eleQ, 
            (kl/2)}};

arc = 
        Apply[Join, 
          Table[fodo, {i, (nfodo/2) - 1}]];

ir1 = {{0, eleQ, k1}, {l1, eleD, l1}, {0, eleQ, k2}, {l2, eleD, 
            l2}};

ir2 = Reverse[ir1];

lattice = 
        Join[arc, fodo1, fodo, ir1, ir2, fodo, fodo1, arc];

twiss0 = twiss6[getm[matform[fodo, 3]]];
βxlat = {{0, twiss0[[1, 2]]}};
βylat = {{0, twiss0[[2, 2]]}};
αxlat = {{0, twiss0[[1, 1]]}};
αylat = {{0, twiss0[[2, 1]]}};
γxlat = {{0, twiss0[[1, 3]]}};
γylat = {{0, twiss0[[2, 3]]}};
ψxlat = {{0, 0}};
ψylat = {{0, 0}};
ηxlat = {{0, twiss0[[1, 5]]}};
ηpxlat = {{0, twiss0[[1, 6]]}};
ηylat = {{0, twiss0[[2, 5]]}};
ηpylat = {{0, twiss0[[2, 6]]}};
proper[matform[Join[fodo1, fodo, ir1], 3], 3];
Show[Graphics[{RGBColor[1, 0, 0], Line[ηxlat]}], 
        Frame -> True, Axes -> False, 
        PlotRange -> {All, All}];
{ηxlat[[-1, 2]], ηpxlat[[-1, 
        2]]}

Show[
      Graphics[{RGBColor[1, 0, 0], Line[βxlat], RGBColor[0, 1, 0], 
          Line[βylat]}], Frame -> True, Axes -> False, 
      PlotRange -> {All, {0, 4}}]

twiss0 = twiss6[getm[matform[fodo, 3]]];
βxlat = {{0, twiss0[[1, 2]]}};
βylat = {{0, twiss0[[2, 2]]}};
αxlat = {{0, twiss0[[1, 1]]}};
αylat = {{0, twiss0[[2, 1]]}};
γxlat = {{0, twiss0[[1, 3]]}};
γylat = {{0, twiss0[[2, 3]]}};
ψxlat = {{0, 0}};
ψylat = {{0, 0}};
ηxlat = {{0, twiss0[[1, 5]]}};
ηpxlat = {{0, twiss0[[1, 6]]}};
ηylat = {{0, twiss0[[2, 5]]}};
ηpylat = {{0, twiss0[[2, 6]]}};
proper[matform[lattice, 3], 3];
Show[Graphics[{RGBColor[1, 0, 0], Line[ηxlat]}], 
        Frame -> True, Axes -> False, 
        PlotRange -> {All, All}];
{ηxlat[[-1, 2]], ηpxlat[[-1, 
        2]]}

Show[
      Graphics[{RGBColor[1, 0, 0], Line[βxlat], RGBColor[0, 1, 0], 
          Line[βylat]}], Frame -> True, Axes -> False, 
      PlotRange -> {All, {0, 4}}]

Chromaticity

nfodo = 20;
len = 1. ;
kl = 1. ;
l1 = 0.1391000706229664;
l2 = 0.3598037760308997;
k1 = 3.645324218395703;
k2 = -5.175733201180464;
Clear[k2la, k2lb]
ϕ0 = ((2  π)/(nfodo - 2));
fodo = {{0, eleS, (k2la/2)}, {0, eleQ, (kl/2)}, {(len/2), eleD, 
          (len/2)}, {0, eleB, ϕ0}, {(len/2), eleD, (len/2)}, {0, 
          eleQ, -(kl/2)}, {0, eleS, k2lb}, {0, 
          eleQ, -(kl/2)}, {(len/2), eleD, (len/2)}, {0, 
          eleB, ϕ0}, {(len/2), eleD, (len/2)}, {0, eleQ, (kl/2)}, {0, eleS, 
          (k2la/2)}}; 
    fodo1 = {{0, eleS, (k2la/2)}, {0, eleQ, (kl/2)}, {(len/2), eleD, 
          (len/2)}, {0, eleB, 0}, {(len/2), eleD, (len/2)}, {0, 
          eleQ, -(kl/2)}, {0, eleS, k2lb}, {0, 
          eleQ, -(kl/2)}, {(len/2), eleD, (len/2)}, {0, eleB, 
          0}, {(len/2), eleD, (len/2)}, {0, eleQ, (kl/2)}, {0, eleS, 
          (k2la/2)}};

arc = 
        Apply[Join, 
          Table[fodo, {i, (nfodo/2) - 1}]];

ir1 = {{0, eleQ, k1}, {l1, eleD, l1}, {0, eleQ, k2}, {l2, eleD, 
            l2}};

ir2 = Reverse[ir1];

lattice = 
        Join[arc, fodo1, fodo, ir1, ir2, fodo, fodo1, arc];

k2la = 0;
k2lb = 0;
twiss0 = twiss6[getm[matform[fodo, 3]]];
βxlat = {{0, twiss0[[1, 2]]}};
βylat = {{0, twiss0[[2, 2]]}};
αxlat = {{0, twiss0[[1, 1]]}};
αylat = {{0, twiss0[[2, 1]]}};
γxlat = {{0, twiss0[[1, 3]]}};
γylat = {{0, twiss0[[2, 3]]}};
ψxlat = {{0, 0}};
ψylat = {{0, 0}};
ηxlat = {{0, twiss0[[1, 5]]}};
ηpxlat = {{0, twiss0[[1, 6]]}};
ηylat = {{0, twiss0[[2, 5]]}};
ηpylat = {{0, twiss0[[2, 6]]}};
proper[matform[lattice, 3], 3];
Show[Graphics[{RGBColor[1, 0, 0], Line[ηxlat]}], 
        Frame -> True, Axes -> False, 
        PlotRange -> {All, All}];
Show[
      Graphics[{RGBColor[1, 0, 0], Line[βxlat], RGBColor[0, 1, 0], 
          Line[βylat]}], Frame -> True, Axes -> False, 
      PlotRange -> {All, {0, 4}}]

Off[General::"spell1"]
chrom[k2_]:=Module[{ξx = 0, ξy = 0},k2la = k2[[1]];k2lb = k2[[2]];Do[If[lattice[[i, 2]] == eleQ,ξx=ξx-enter your formulas for natural chromaticity;ξy=ξy+enter your formulas for natural chromaticiy];If[lattice[[i, 2]] == eleS,ξx=ξx+enter your formulas for sextupole contribution;ξy=ξy-enter your formulas for sextupole contribution],{i, Length[lattice]}];{ξx, ξy}]
On[
      General::"spell1"]

fit[chrom, {0, 0}, {0.8, -1.24}, 25]

Resonance strength

nfodo = 20;
len = 1. ;
kl = 1. ;
l1  = 0.1390925589102365;
l2  = 0.3597201861389145;
k1  = 3.6456962918177758;
k2  = -5.176530688784192;
k2la = 0.8422392259334442;
k2lb = -1.428721413036512;
ϕ0 = ((2  π)/(nfodo - 2));
fodo = {{0, eleS, (k2la/2)}, {0, eleQ, (kl/2)}, {(len/2), eleD, 
          (len/2)}, {0, eleB, ϕ0}, {(len/2), eleD, (len/2)}, {0, 
          eleQ, -(kl/2)}, {0, eleS, k2lb}, {0, 
          eleQ, -(kl/2)}, {(len/2), eleD, (len/2)}, {0, 
          eleB, ϕ0}, {(len/2), eleD, (len/2)}, {0, eleQ, (kl/2)}, {0, eleS, 
          (k2la/2)}}; 
    fodo1 = {{0, eleS, (k2la/2)}, {0, eleQ, (kl/2)}, {(len/2), eleD, 
          (len/2)}, {0, eleB, 0}, {(len/2), eleD, (len/2)}, {0, 
          eleQ, -(kl/2)}, {0, eleS, k2lb}, {0, 
          eleQ, -(kl/2)}, {(len/2), eleD, (len/2)}, {0, eleB, 
          0}, {(len/2), eleD, (len/2)}, {0, eleQ, (kl/2)}, {0, eleS, 
          (k2la/2)}};

arc = 
        Apply[Join, 
          Table[fodo, {i, (nfodo/2) - 1}]];

ir1 = {{0, eleQ, k1}, {l1, eleD, l1}, {0, eleQ, k2}, {l2, eleD, 
            l2}};

ir2 = Reverse[ir1];

lattice = 
        Join[arc, fodo1, fodo, ir1, ir2, fodo, fodo1, arc];

twiss0 = twiss6[getm[matform[fodo, 3]]];
βxlat = {{0, twiss0[[1, 2]]}};
βylat = {{0, twiss0[[2, 2]]}};
αxlat = {{0, twiss0[[1, 1]]}};
αylat = {{0, twiss0[[2, 1]]}};
γxlat = {{0, twiss0[[1, 3]]}};
γylat = {{0, twiss0[[2, 3]]}};
ψxlat = {{0, 0}};
ψylat = {{0, 0}};
ηxlat = {{0, twiss0[[1, 5]]}};
ηpxlat = {{0, twiss0[[1, 6]]}};
ηylat = {{0, twiss0[[2, 5]]}};
ηpylat = {{0, twiss0[[2, 6]]}};
proper[matform[lattice, 3], 3];
Show[Graphics[{RGBColor[1, 0, 0], Line[ηxlat]}], 
        Frame -> True, Axes -> False, 
        PlotRange -> {All, All}];
Show[
        Graphics[{RGBColor[1, 0, 0], Line[βxlat], RGBColor[0, 1, 0], 
            Line[βylat]}], Frame -> True, Axes -> False, 
        PlotRange -> {All, {0, 4}}] 
      Show[Graphics[{RGBColor[1, 0, 0], Line[ψxlat], RGBColor[0, 1, 0], 
            Line[ψylat]}], Frame -> True, Axes -> False, 
        PlotRange -> {All, All}]

Length and phase advance around the ring:

ψxlat[[-1]]
ψylat[[-1]]

Off[General::"spell"]
len = βxlat[[-1, 1]];
μx = ψxlat[[-1, 2]] 
          2  π;
μy = ψylat[[-1, 2]] 
          2  π;
On[General::"spell"]

Insert your own formulas here for resonance strength calculations.

Off[General::"spell1"]
res[] :=   ...  
On[General::"spell1"]

res[]

ot[3]

Tune

nfodo = 20;
len = 1. ;
kl = 1. ;
l1  = 0.1390925589102365;
l2  = 0.3597201861389145;
k1  = 3.6456962918177758;
k2  = -5.176530688784192;
k2la = 0.8422392259334442;
k2lb = -1.428721413036512;
ϕ0 = ((2  π)/(nfodo - 2));
fodo = {{0, eleS, (k2la/2)}, {0, eleQ, (k1a/2)}, {(len/2), eleD, 
          (len/2)}, {0, eleB, ϕ0}, {(len/2), eleD, (len/2)}, {0, eleQ, 
          (k1b/2)}, {0, eleS, k2lb}, {0, eleQ, (k1b/2)}, {(len/2), eleD, 
          (len/2)}, {0, eleB, ϕ0}, {(len/2), eleD, (len/2)}, {0, eleQ, 
          (k1a/2)}, {0, eleS, (k2la/2)}}; 
    fodo1 = {{0, eleS, (k2la/2)}, {0, eleQ, (k1a/2)}, {(len/2), eleD, 
          (len/2)}, {0, eleB, 0}, {(len/2), eleD, (len/2)}, {0, eleQ, 
          (k1b/2)}, {0, eleS, k2lb}, {0, eleQ, (k1b/2)}, {(len/2), eleD, 
          (len/2)}, {0, eleB, 0}, {(len/2), eleD, (len/2)}, {0, eleQ, 
          (k1a/2)}, {0, eleS, (k2la/2)}};

arc = 
        Apply[Join, 
          Table[fodo, {i, (nfodo/2) - 1}]];

ir1 = {{0, eleQ, k1}, {l1, eleD, l1}, {0, eleQ, k2}, {l2, eleD, 
            l2}};

ir2 = Reverse[ir1];

lattice = 
        Join[arc, fodo1, fodo, ir1, ir2, fodo, fodo1, arc];

tune[kq_] := 
      Module[{twiss}, k1a = kq[[1]]; k1b = kq[[2]]; 
        twiss = twiss6[getm[matform[lattice, 3]]]; {twiss[[1, 4]], 
          twiss[[2, 4]]}]

tune[{1. , -1. }]

fit[tune, {0.28, 0.29}, {1. , -1. }, 20]

k1a = 1.0548344347429697;
k1b = -1.0411961165546295;

twiss0 = twiss6[getm[matform[lattice, 3]]];
βxlat = {{0, twiss0[[1, 2]]}};
βylat = {{0, twiss0[[2, 2]]}};
αxlat = {{0, twiss0[[1, 1]]}};
αylat = {{0, twiss0[[2, 1]]}};
γxlat = {{0, twiss0[[1, 3]]}};
γylat = {{0, twiss0[[2, 3]]}};
ψxlat = {{0, 0}};
ψylat = {{0, 0}};
ηxlat = {{0, twiss0[[1, 5]]}};
ηpxlat = {{0, twiss0[[1, 6]]}};
ηylat = {{0, twiss0[[2, 5]]}};
ηpylat = {{0, twiss0[[2, 6]]}};
proper[matform[lattice, 3], 3];
Show[Graphics[{RGBColor[1, 0, 0], Line[ηxlat]}], 
        Frame -> True, Axes -> False, 
        PlotRange -> {All, All}];
Show[
        Graphics[{RGBColor[1, 0, 0], Line[βxlat], RGBColor[0, 1, 0], 
            Line[βylat]}], Frame -> True, Axes -> False, 
        PlotRange -> {All, {0, 4}}] 
      Show[Graphics[{RGBColor[1, 0, 0], Line[ψxlat], RGBColor[0, 1, 0], 
            Line[ψylat]}], Frame -> True, Axes -> False, 
        PlotRange -> {All, All}]

ψxlat[[-1]]

ψylat[[-1]]

len = βxlat[[-1, 1]];
μx = ψxlat[[-1, 2]] 
          2  π;
μy = ψylat[[-1, 2]] 2  π;

To be exact, now the dispersion could be rematched, the chromaticity rematched, and the resonance strength recomputed. but all these values do not change much.

Dynamic Aperture

Linear optics tracking:

mtrack[mat_, zi_] := mat . zi

zi = {0.0001, 0. };

mat = getm[matform[lattice, 1]];

tab = {zi}; 
    Do[tab = Append[tab, mtrack[mat, tab[[-1]]]], {i, 
        100}]

Show[Graphics[Map[Point[#] &, tab]], AspectRatio -> 1, 
      Frame -> True]

klat = kickform[lattice, 1];

Off[General::"spell1"]
tt = 
        Table[zi = {0.01, 0. } i; tab = {zi}; 
          Do[tab = Append[tab, lattrack[klat, tab[[-1]]]], {i, 
              100}]; tab, {i, 1, 20}];
Show[Graphics[Map[Map[Point[#] &, #] &, tt]], AspectRatio -> 1, 
      Frame -> True]
On[General::"spell1"]