\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins,amsmath,amssymb,bm,hyperref,graphicx,listings}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\pdd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdds}[2]{\frac{\partial^2 #1}{\partial{#2}^2}}
\newcommand{\pddd}[3]{\frac{\partial^2 #1}{\partial #2\partial #3}}

\newcommand{\rowfour}[4]{\left[\begin{array}{cccc}#1&#2&#3&#4\end{array}\right]}
\newcommand{\vectwo}[2]{\left[\begin{array}{c}#1\\#2\end{array}\right]}
\newcommand{\mxtwo}[4]{\left[\begin{array}{cc}#1&#2\\#3&#4\end{array}\right]}

\newcommand{\sq}[1]{\left(#1\right)^2}
\newcommand{\floor}[1]{\left\lfloor#1\right\rfloor}
\newcommand{\abs}[1]{\left|#1\right|}

\newcommand{\rem}[1]{}

\newcommand{\A}{\mathbf A}
\newcommand{\E}{\mathbf E}
\newcommand{\B}{\mathbf B}
\newcommand{\J}{\mathbf J}
\newcommand{\0}{\mathbf 0}

\newcommand{\f}{\mathbf f}
\newcommand{\uv}{\mathbf u}
\newcommand{\vv}{\mathbf v}
\newcommand{\w}{\mathbf w}
\newcommand{\x}{\mathbf x}
\newcommand{\y}{\mathbf y}
\newcommand{\z}{\mathbf z}

\newcommand{\s}{\mathbf s}
\newcommand{\rv}{\mathbf r}
\newcommand{\cv}{\mathbf c}

\newcommand{\e}{\mathbf e}
\newcommand{\m}{\mathbf m}
\newcommand{\M}{\mathbf M}

\newcommand{\de}{\,\mathrm d}

\newcommand{\ii}{\mathrm i}

\begin{document}

\title{Benchmark Midplane Field for Medical FFA}
\author{Stephen Brooks, Dejan Trbojevic}
\date{April 2, 2024}
\maketitle

\section{Introduction}
This is a lattice cell for a fixed-field medical accelerator, transmitting protons from 10--250\,MeV.  The goal is to have fixed tunes throughout this energy range.  The cell is based on the lattice \texttt{Feb7\_2024\_ProtonAccelerator.f90} from Dejan Trbojevic and had the magnet field profiles re-optimised to attain fixed tunes with realistic magnet fringe fields.

\subsection{Coordinate Convention}
$z$ is the forward longitudinal coordinate, $x$ is transverse (orbit excursion) and $y$ is vertical.  The curvature of the machine is anticlockwise as seen from above (positive $y$ side).  $x$ points towards the machine centre, so the high energy orbits are at more negative $x$.

\section{Cell Geometry}
\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{cell_axes}
   \caption{Closed orbits in \texttt{Muon1} on top of the cell geometry and 1\,cm grid.}
   \label{fig:orbits}
\end{figure}

The cell as drawn in Figure~\ref{fig:orbits} has element start and end coordinates given in Table~\ref{tab:coords}.

\begin{table}[!hbt]
   \centering
   \caption{Element start and end positions in absolute coordinates.}
   \begin{tabular}{lllccc}
       \hline
\textbf{Element name} & \textbf{Length (m)} & \textbf{Angle (rad)} & \\
       \hline
HO & 0.0045 & 0 \\
F & 0.212 & 0.0639198516350064 \\
O & 0.009 & 0 \\
D & 0.113611696486384 & 0.0342549187896746 \\
HO & 0.0045 & 0 \\
       \hline
       \hline
\textbf{Element name} & \textbf{Start $x$ (m)} & \textbf{$y$} & \textbf{$z$} \\
       \hline
HO & 0 & 0 & 0 \\
F & -0.00112915456451474 & 0 & 0.00451804630577078 \\
O & 0.00677319767087102 & 0 & 0.216355666412717 \\
D & 0.00702461901507897 & 0 & 0.22536077394586 \\
HO & 0.0165454801670591 & 0 & 0.33857051247452 \\
       \hline
       \hline
\textbf{Element name} & \textbf{End $x$ (m)} & \textbf{$y$} & \textbf{$z$} \\
       \hline
HO & 0 & 0 & 0.0045 \\
F & 0.00564519630940407 & 0 & 0.216409783214592 \\
O & 0.00734808467505134 & 0 & 0.225337286808416 \\
D & 0.0162224641982652 & 0 & 0.338599535960023 \\
HO & 0.0169865572985422 & 0 & 0.343048843744545 \\
       \hline
   \end{tabular}
   \label{tab:coords}
\end{table}

The cell starts at the beginning of the first half drift \texttt{HO} at $(0,0,0)$.  It ends at the end of the second half drift \texttt{HO}, with total cell rotation angle 0.09817477\,rad = 5.625$^\circ$ = 360$^\circ$/64. 

The ``Length'' and ``Angle'' values define a sequence of circular arcs and straight lines that form the reference axis for the cell.  However, the rectangular magnets are only aligned `on average' on the curved parts, so the start and end coordinates for the magnets do not lie on the reference axis.

\clearpage

\section{Midplane ($y=0$) Magnetic Field Model}
The field of each magnet is defined in terms of local Cartesian coordinates $z_{loc}$ and $x_{loc}$, where $(z_{loc},x_{loc})=(0,0)$ is the start coordinate of the magnet and $(z_{loc},x_{loc})=(L,0)$ is the end.

The midplane field is vertical and the product of three functions:
\[B_y(x,0,z) = g(x_{loc})f(z_{loc})f(L-z_{loc}), \qquad B_x(x,0,z) = B_z(x,0,z) = 0, \]
where $g$ defines the body field and $f$ are fringe field sigmoid-type functions.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{Dejan_medicalfixedtune_LDRD_2024-Feb-07_fourier8atan_opt_83500}
   \caption{Colour plot of the vertical field component in the cell, colour contours are 0.1\,T.}
   \label{fig:fields}
\end{figure}

Figure \ref{fig:fields} shows the field in the cell.

\subsection{Body Field Function}
To avoid ill-conditioned optimisation with high-order polynomials, the body field function is defined using a Fourier series:
\[ g(x) = c_0 + \sum_{n=1}^8 s_n\sin(nkx)+c_n\cos(nkx),  \]
where $k$ is a transverse scale factor and $c_0...c_8$ and $s_1...s_8$ are the coefficients defining the field (units Tesla).

\begin{table}[!hbt]
   \centering
   \caption{Fourier coefficients for the magnet body fields.}
   \begin{tabular}{lcc}
       \hline
\textbf{Magnet F} & \textbf{$k = 60$ (m$^{-1}$)} & \\
\textbf{$n$} & \textbf{$c_n$ (T)} & \textbf{$s_n$ (T)} \\
       \hline
0 & -0.69802703489829 & --- \\
1 & 0.0756049884526512 & 1.73518865753323 \\
2 & 0.305088336939702 & -0.442262271952788 \\
3 & -0.139524713341402 & 0.0826714806226137 \\
4 & 0.0276793694098833 & -0.0208114430730033 \\
5 & -0.00196329427184875 & 0.0191643586396017 \\
6 & 0.00341848475778073 & -0.0112392171957729 \\
7 & -0.00264986640664568 & 0.00271130619879153 \\
8 & 0.000542396377167934 & -6.52682264979307E-05 \\
       \hline
       \hline
\textbf{Magnet D} & \textbf{$k = 60$ (m$^{-1}$)} & \\
\textbf{$n$} & \textbf{$c_n$ (T)} & \textbf{$s_n$ (T)} \\
       \hline
0 & 0.989342248880683 & --- \\
1 & -1.91337116068355 & -2.46074568747532 \\
2 & 0.397559249919083 & 0.493240962238742 \\
3 & -0.0885340223614818 & -0.134280872283561 \\
4 & 0.0116131686129637 & 0.0706611962587838 \\
5 & -0.010779760306873 & 0.00724018936522096 \\
6 & 0.0207629971432377 & -0.0480848747522068 \\
7 & -0.014783458990518 & 0.0306916335109993 \\
8 & 0.00390264678161641 & -0.00685949055711114 \\
       \hline
   \end{tabular}
   \label{tab:fouriercoeffs}
\end{table}

The values for these coefficients of the two magnets are given in Table~\ref{tab:fouriercoeffs}.

\subsection{Fringe Field Function}
It was found that functions of the form
\[ f(z) = \frac{\frac1a\arctan(z/a)-\frac1b\arctan(z/b)}{\left(\frac1a-\frac1b\right)\pi}+\frac12 \]
are a good fit to real permanent magnet fringe fields.

Both magnets in this cell use the values $a=0.008$\,m and $b=0.04$\,m.

The use of arc tangent is physically motivated from the field of two half-infinite current sheets with the beam emerging between them.  This would have a long-range decay as $1/z$, which is too slow, so the difference of two arc tangents is used, cancelling the $1/z$ term and leaving the remainder decaying as $1/z^3$.  This is the field of two current sheets that are infinite transversely but not vertically.

This fringe field is symmetrical so magnets can be split into longitudinal segments while keeping the total field the same. 

\appendix
\renewcommand{\thesection}{Appendix \Alph{section}}

\clearpage

\section{\texttt{Muon1} Lattice Input File}
\lstset{basicstyle=\scriptsize\ttfamily}
\begin{lstlisting}[breaklines]
{Dejan's medical fixed-tune NSFFA for LDRD, design 2024-Feb-07
Proper Fourier approximation (order 5), by least squares fit, using k=60m^-1, zero extended to order 8
But with TwoAtan fringe fields fitted to the 3D magnet model, increase FringeLimit to 50cm
Optimised values including max |Bmidl|}

#// Dejan/medicalfixedtune_LDRD_2024-Feb-07 - Composite cell of 17 elements:

{MatchScan AllowUnstable Egoal Estart Species}
Match 1 1.602176634e-012 4.005441585e-011 'Proton'
{Drift Length}
HO 0.0045
{Magnet3Rect AlignMode Angle FringeFunc FringeFunca FringeFuncb FringeLimit func3 func3c0 func3s1 func3c1 func3s2 func3c2 func3s3 func3c3 func3s4 func3c4 func3s5 func3c5 func3s6 func3c6 func3s7 func3c7 func3s8 func3c8 func3order func3k HalfHeight HalfWidth Length}
F 'Integral' 0.0639198516350064 'TwoAtan' 0.008 0.04 0.5 'Fourier' -0.69802703489829 1.73518865753323 0.0756049884526512 -0.442262271952788 0.305088336939702 0.0826714806226137 -0.139524713341402 -0.0208114430730033 0.0276793694098833 0.0191643586396017 -0.00196329427184875 -0.0112392171957729 0.00341848475778073 0.00271130619879153 -0.00264986640664568 -6.52682264979307E-05 0.000542396377167934 8 60 0.01 0.07 0.212
{Drift Length}
O 0.009
{Magnet3Rect AlignMode Angle FringeFunc FringeFunca FringeFuncb FringeLimit func3 func3c0 func3s1 func3c1 func3s2 func3c2 func3s3 func3c3 func3s4 func3c4 func3s5 func3c5 func3s6 func3c6 func3s7 func3c7 func3s8 func3c8 func3order func3k HalfHeight HalfWidth Length}
D 'Integral' 0.0342549187896746 'TwoAtan' 0.008 0.04 0.5 'Fourier' 0.989342248880683 -2.46074568747532 -1.91337116068355 0.493240962238742 0.397559249919083 -0.134280872283561 -0.0885340223614818 0.0706611962587838 0.0116131686129637 0.00724018936522096 -0.010779760306873 -0.0480848747522068 0.0207629971432377 0.0306916335109993 -0.014783458990518 -0.00685949055711114 0.00390264678161641 8 60 0.01 0.07 0.113611696486384
{Drift Length}
HO 0.0045
{Match-Aperture}
MatchEnd

Dejan_medicalfixedtune_LDRD_2024_Feb_07: HO,F,O,D,HO,Match,HO,F,O,D,HO,MatchEnd,HO,F,O,D,HO;
\end{lstlisting}


\end{document}