\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins,hyperref}

\newcommand{\rem}[1]{}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\ddh}[2]{\mathrm d #1/\mathrm d #2}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\ddhp}[3]{\mathrm d^{#3} #1/{\mathrm d #2}^{#3}}
\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{\vecthree}[3]{\left[\begin{array}{c}#1\\#2\\#3\end{array}\right]}
\newcommand{\mxtwo}[1]{\left[\begin{array}{cc}#1\end{array}\right]}
\newcommand{\mxthree}[1]{\left[\begin{array}{ccc}#1\end{array}\right]}

\newcommand{\sq}[1]{\left(#1\right)^2}

\newcommand{\ii}{\mathrm i}
\newcommand{\e}{\mathrm e}

\newcommand{\B}{\mathbf B}
\newcommand{\C}{\mathbf C}
\newcommand{\F}{\mathbf F}
\newcommand{\x}{\mathbf x}
\newcommand{\p}{\mathbf p}
\newcommand{\vv}{\mathbf v}
\newcommand{\uv}{\mathbf u}

\begin{document}

\title{Extrapolation of Magnetic Fields from a Curved Surface}
\author{Stephen Brooks}
\maketitle

\section{Flat Surface}
Maxwell's equations in free space, $\nabla\cdot\B=0$ and $\nabla\times\B=\mathbf0$, can be rearranged to give
\setlength{\arraycolsep}{3pt}
\[\partial_y\B = \mxthree{0&\partial_x&0\\-\partial_x&0&-\partial_z\\0&\partial_z&0}\B, \]
which can be used in a Taylor expansion of the form
\[ \B(x,y,z) = \sum_{n=0}^\infty \frac{y^n}{n!}\partial_y^n\B(x,0,z) \]
to obtain fields for $y\ne0$ in terms of the derivatives of $\B$ on the $y=0$ plane.  The repeated derivatives $\partial_y^n$ must be expressed as combinations of $\partial_x$ and $\partial_z$ in order for the two-dimensional function $\B(x,0,z)$ to hold all the relevant information for extrapolating the field.  This is more easily expressed by breaking $\B$ into $B_y$ and $\B_{x,z}=(B_x,B_z)$ so that
\begin{eqnarray*}
&&\partial_y B_y = -\nabla_{x,z}^T\B_{x,z} \qquad\mbox{and}\qquad \partial_y \B_{x,z} = \nabla_{x,z}B_y \\
\Rightarrow&&\partial^2_y B_y = -\nabla_{x,z}^T\partial_y\B_{x,z} = -\nabla_{x,z}^T\nabla_{x,z}B_y = -\nabla^2_{x,z}B_y,
\end{eqnarray*}
where $\nabla^2_{x,z} = \nabla_{x,z}^T\nabla_{x,z} = \partial^2_x+\partial^2_z$.  Continuing to use the fact that all partial derivatives commute, this may be repeated to get
\begin{eqnarray*}
\partial_y^{2n} B_y &=& (-\nabla^2_{x,z})^n B_y \\
\partial_y^{2n+1} B_y &=& -(-\nabla^2_{x,z})^n \nabla_{x,z}^T\B_{x,z} \\ 
\partial_y^{2n+1} \B_{x,z} &=& \nabla_{x,z}(-\nabla^2_{x,z})^n B_y \\
\partial_y^{2n+2} \B_{x,z} &=& -\nabla_{x,z}(-\nabla^2_{x,z})^n \nabla_{x,z}^T\B_{x,z}
\end{eqnarray*}
for $n\ge0$.  Reassembling into matrix form gives
\begin{eqnarray*}
\partial^{2n+1}_y \B &=& (-\nabla^2_{x,z})^n \mxthree{0&\partial_x&0\\-\partial_x&0&-\partial_z\\0&\partial_z&0}\B \qquad \mbox{and} \\
\partial^{2n+2}_y \B &=& (-\nabla^2_{x,z})^n \mxthree{-\partial_x^2&0&-\partial_x\partial_z\\0&-\nabla^2_{x,z}&0\\-\partial_x\partial_z&0&-\partial_z^2}\B.
\end{eqnarray*}
More explicit formulae can be obtained by using the binomial expansion of $(\partial^2_x+\partial^2_z)^n$:
\[ (-\nabla^2_{x,z})^n = (-1)^n \sum_{i=0}^n {n \choose i} \partial^{2i}_x\partial^{2(n-i)}_z. \]

\section{Curved Surface with $Y$ Displacement Function\label{-dispy}}
Instead of tackling fields specified on a fully-general curved surface, this section considers the case where the surface can be represented as the graph of a two-argument function $y=Y(x,z)$.  It will help to define a transformed field function
\[ \C(x,y,z) = \B(x,y+Y(x,z),z) \qquad \Leftrightarrow \qquad \B(x,y,z) = \C(x,y-Y(x,z),z) \]
so that the initial conditions are the two-dimensional function $\C(x,0,z) = \B(x,Y(x,z),z)$.  The strategy will be to evaluate fields via Taylor expansion as before:
\[ \B(x,y,z) = \C(x,y-Y(x,z),z) = \sum_{n=0}^\infty \frac{(y-Y(x,z))^n}{n!}\partial_y^n\C(x,0,z), \]
where now the $y$ derivatives of $\C$ need to be related to the $x$ and $z$ derivatives of $\C(x,0,z)$ via Maxwell's equations, which are normally stated in terms of $\B$.

\subsection{Transformation of Partial Derivative Operators}
The chain rule gives the derivative rules below, which can be applied to any component of $\B$ and $\C$ as well as the full vectors:
\begin{eqnarray*}
\partial_x \B(x,y,z) &=& \partial_x(\C(x,y-Y(x,z),z)) = \partial_x\C - \pdd{Y}{x}\partial_y\C = \left(\partial_x - \pdd{Y}{x}\partial_y\right)\C; \\
\partial_y \B &=& \partial_y \C; \\
\partial_z \B &=&  \left(\partial_z - \pdd{Y}{z}\partial_y\right)\C.
\end{eqnarray*}
Henceforth $\B$ and $\C$ are evaluated at their respective transformed locations unless otherwise specified.  Note that the formula for $\C$ in terms of $\B$ is identical to the formula going the other way except for having $+Y$ instead of $-Y$, so the reverse transformations are simply
\begin{eqnarray*}
\partial_x \C(x,y,z) &=& \left(\left(\partial_x + \pdd{Y}{x}\partial_y\right)\B\right)(x,y+Y(x,z),z); \\
\partial_y \C &=& \partial_y \B; \\
\partial_z \C &=&  \left(\partial_z + \pdd{Y}{z}\partial_y\right)\B.
\end{eqnarray*}

\subsection{Maxwell's Equations in Transformed Field}
The purpose of this section is to substitute the above rules to express the free-space Maxwell's equations in terms of derivatives of $\C$'s components and then rearrange to get $\partial_y\C$ in terms of $\C$'s $x$ and $z$ derivatives.  For compactness, the notation $Y_x = \pdd{Y}{x}$ and $Y_z = \pdd{Y}{z}$ will be used.  First, the substitution:
\[ 0 = \nabla\cdot\B = \partial_x B_x + \partial_y B_y + \partial_z B_z
= \left(\partial_x - Y_x\partial_y\right)C_x + \partial_y C_y + \left(\partial_z - Y_z\partial_y\right)C_z \]
\[ \Rightarrow \qquad \partial_y\left( - Y_xC_x + C_y - Y_zC_z\right) = -\partial_x C_x - \partial_z C_z. \]
The equation $\nabla\times\B=\mathbf0$ has three components:
\[ \partial_y B_x = \partial_x B_y
\quad \Rightarrow \quad 
\partial_y C_x = \left(\partial_x - Y_x\partial_y\right) C_y
\quad \Rightarrow \quad 
\partial_y \left(C_x + Y_xC_y\right) = \partial_x C_y; \]
\[ \partial_y B_z = \partial_z B_y
\quad \Rightarrow \quad 
\partial_y C_z = \left(\partial_z - Y_z\partial_y\right) C_y
\quad \Rightarrow \quad 
\partial_y \left(C_z + Y_zC_y\right) = \partial_z C_y. \]
The third equation component enforces a consistency condition within the initial conditions:
\begin{eqnarray*}
\partial_x B_z = \partial_z B_x
&\quad \Rightarrow \quad&
\left(\partial_x - Y_x\partial_y\right) C_z = \left(\partial_z - Y_z\partial_y\right) C_x \\ 
&\quad \Rightarrow \quad&
\partial_y \left(Y_z C_x - Y_x C_z\right) = \partial_z C_x - \partial_x C_z,
\end{eqnarray*}
or at least it would do in the case $Y=0$; here the coordinate change has mixed the equations.  The four equations can be written in matrix form:
\[
\partial_y \mxthree{
 -Y_x & 1 & -Y_z \\
1 & Y_x & 0 \\
0 & Y_z & 1 \\
Y_z & 0 & -Y_x} \C = 
\mxthree{
-\partial_x & 0 & -\partial_z \\
0 & \partial_x & 0 \\
0 & \partial_z & 0 \\
\partial_z & 0 & -\partial_x} \C.
\]
Taking $-Y_z$ times the second row, $Y_x$ times the third row and adding the fourth row gives the true consistency condition
\[ 0 = \partial_z C_x - Y_z \partial_x C_y + Y_x \partial_z C_y - \partial_x C_z, \]
which involves only $\partial_x$ and $\partial_z$ of $\C$.  This condition must be checked on $\C(x,0,z)$ before extrapolating otherwise the field cannot be Maxwellian.

Deleting the fourth row and swapping the first and second gives a familiar-looking relation
\[
\mxthree{
1 & Y_x & 0 \\
 -Y_x & 1 & -Y_z \\
0 & Y_z & 1} \partial_y \C = 
\mxthree{
0 & \partial_x & 0 \\
-\partial_x & 0 & -\partial_z \\
0 & \partial_z & 0} \C,
\]
the only complicating factor being the matrix premultiplied on the left-hand side.  It does have an inverse, so
\[
\partial_y \C = 
\frac{1}{1+Y_x^2+Y_z^2}
\mxthree{
1+Y_z^2&-Y_x&-Y_xY_z \\
Y_x&1&Y_z \\
-Y_xY_z&-Y_z&1+Y_x^2}
\mxthree{
0 & \partial_x & 0 \\
-\partial_x & 0 & -\partial_z \\
0 & \partial_z & 0} \C \]
but this is about as much as formal algebra can do.

\subsection{Computation of Repeated Derivatives of $\C$}
The equation in the last section can be written 
\[ \partial_y \C = M(A \partial_x + B \partial_z) \C \]
for a $3\times3$ variable matrix $M(x,z)$ and two constant matrices $A=\tiny\mxthree{0&1&0\\-1&0&0\\0&0&0},B=\tiny\mxthree{0&0&0\\0&0&-1\\0&1&0}$.  Repeated application of this rule produces a combinatorial explosion of terms that is best handled by a computer.  The inputs for the algebraic part of the computation are the repeated derivatives of $\C$ and $M$; let
\[ \C_{ni} = \partial_x^n\partial_z^i\C \qquad\mbox{and}\qquad M_{ni} = \partial_x^n\partial_z^i M. \]
These will be calculated in a problem-specific way depending on the form of the initial condition and surface.  Note that $M_{ni}$ are many matrices, not elements of a single matrix.

The starting point for the computation is
\[ \partial_y^0 \C = \C = \C_{00}. \]
The first step gives
\[ \partial_y^1 \C = MA\C_{10} + MB\C_{01} = M_{00}A\C_{10} + M_{00}B\C_{01}, \]
where subscripts are used on $M_{00}$ since it will soon be subject to differentiation.  The next step is
\begin{eqnarray*}
 \partial_y^2 \C &=& \left(M_{00}A\partial_x + M_{00}B\partial_z\right)\partial_y^1 \C \\
&=& M_{00}AM_{10}A\C_{10} + M_{00}AM_{00}A\C_{20} + M_{00}AM_{10}B\C_{01} + M_{00}AM_{00}B\C_{11} \\
&+& M_{00}BM_{01}A\C_{10} + M_{00}BM_{00}A\C_{11} + M_{00}BM_{01}B\C_{01} + M_{00}BM_{00}B\C_{02}.
\end{eqnarray*}
It becomes clear that the normal form of a term in $\partial_y^n \C$ is $n$ repeats of $M_{ni}\{A,B\}$ followed by a $\C_{ni}$.  If there is no collection of terms, the expression has $2^n$ branches depending on whether $M_{00}A\partial_x$ or $M_{00}B\partial_z$ is applied at each step.  However, each branch has $n!$ terms since each term in $\partial_y^{n-1} \C$ has $n$ targets for the derivative operator of a product to act on ($n-1$ $M$'s and one $\C$).  This gives a total of $2^nn!$ terms, which equals 1, 2, 8, 48, 384, 3840 for $n=0$ to $5$.  The action of the derivative operator $\partial_x$ is simply to add one to the first index of each $M$ or $\C$ in the term in turn, summing the results.  $\partial_z$ adds one to the second index of each.

Although the number of algebraic terms explodes, the number of inputs is more well-behaved.  For $n\ge1$, $\partial_y^n\C$ requires $M_{ij}$ for $0\le i+j\le n-1$ and $\C_{ij}$ for $1\le i+j\le n$, which increase in number only quadratically.

\subsection{$M_{ni}$ Matrices in Terms of Derivatives of $Y$}
The $M_{ni}$ matrices are not directly the input to the algorithm, instead they are made of repeated derivatives of $Y(x,z)$.  Adopting the notation
\[ Y_{ni} = \partial_x^n\partial_z^i Y \qquad\Rightarrow\qquad Y_x=Y_{10}, \quad Y_z=Y_{01}, \]
the first matrix is:
\[
M_{00} = 
\frac{1}{1+Y_{10}^2+Y_{01}^2}
\mxthree{
1+Y_{01}^2&-Y_{10}&-Y_{10}Y_{01} \\
Y_{10}&1&Y_{01} \\
-Y_{10}Y_{01}&-Y_{01}&1+Y_{10}^2}.
\]
Further partial derivatives can be done component-wise in the matrix and term-wise in any sum, so the problem reduces to evaluating derivatives of the general term
\[ T = \frac{\prod_i Y_{a_ib_i}}{\left(1+Y_{10}^2+Y_{01}^2\right)^n}. \]
Using the chain and product rules,
\begin{eqnarray*}
\pdd{T}{x} &=& \frac{-n}{\left(1+Y_{10}^2+Y_{01}^2\right)^{n+1}}\left(2Y_{10}Y_{20}+2Y_{01}Y_{11}\right)\prod_i Y_{a_ib_i} + \frac{\sum_i Y_{a_i+1,b_i} \prod_{j\ne i} Y_{a_jb_j}}{\left(1+Y_{10}^2+Y_{01}^2\right)^n} \\
&=& \frac{-2n\left(Y_{10}Y_{20}+Y_{01}Y_{11}\right)\prod_i Y_{a_ib_i} + \left(1+Y_{10}^2+Y_{01}^2\right) \sum_i Y_{a_i+1,b_i} \prod_{j\ne i} Y_{a_jb_j}}{\left(1+Y_{10}^2+Y_{01}^2\right)^{n+1}}
\end{eqnarray*}
and similarly
\[ \pdd{T}{z} = \frac{-2n\left(Y_{10}Y_{11}+Y_{01}Y_{02}\right)\prod_i Y_{a_ib_i} + \left(1+Y_{10}^2+Y_{01}^2\right) \sum_i Y_{a_i,b_i+1} \prod_{j\ne i} Y_{a_jb_j}}{\left(1+Y_{10}^2+Y_{01}^2\right)^{n+1}}. \]


\subsection{Evaluation Strategy}
If a magnetic field is specified by initial conditions $\B(x,Y(x,z),z) = \C(x,z)$ on some surface $y=Y(x,z)$, then subject to series convergence,
\[ \B(x,y,z) = \sum_{n=0}^\infty \frac{(y-Y(x,z))^n}{n!}\partial_y^n\C(x,0,z). \]
This section considers a situation where the infinite sum is truncated at order $N$, so the expression below is evaluated:
\[ \B_N(x,y,z) = \sum_{n=0}^N \frac{(y-Y(x,z))^n}{n!}\partial_y^n\C(x,0,z), \]
assuming the derivatives $Y_{ij}, \C_{ij}$ for $i+j\le N$ can be calculated at $(x,z)$ using a subroutine.

\subsubsection{Once-Only Calculations}
\begin{itemize}
\item The algebraic form of the 6 distinct elements of $M_{ij}$ (it is symmetric) for $i+j\le N-1$, in terms of the derivatives of $Y$.
\item The algebraic form of $\partial_y^n\C$ in terms of the $M$ matrices and the derivatives of $\C$ (and the constant matrices $A,B$), for $n\le N$.
\end{itemize}

\subsubsection{Per-Point Calculations\label{-dispy-eval-point}}
\begin{itemize}
\item $Y_{ij}, \C_{ij}$ for $i+j\le N$ at $(x,z)$.
\item $M_{ij}$ for $i+j\le N-1$ using the precalculated algebraic form with $Y_{ij}$.
\item $\partial_y^n\C$ using the precalculated algebraic form with $M_{ij}, \C_{ij}$.
\item $\B_N(x,y,z)$ using the finite Taylor sum formula with $\partial_y^n\C$ at this point.
\end{itemize}

\subsubsection{Additional Condition}
For $\B$ to be Maxwellian, the input function $\C(x,z)$ must satisfy 
\[ \partial_z C_x - \partial_x C_z =  Y_z \partial_x C_y - Y_x \partial_z C_y. \]


\section{Examples and Applications}
This section presents some useful specific cases with formulae for generating their $Y_{ij}$ and $\C_{ij}$ derivatives.

\subsection{Purely Vertical Field at Surface\label{-exai-vert}}
The extrapolation from fields on a flat surface is much simplified if only the $B_y$ component is present on the mid-plane.  In this case, the field is symmetrical about the $y=0$ plane and the Maxwell consistency condition $\partial_x B_z - \partial_z B_x =0$ is always satisfied since those components of $\B$ are zero within the plane.

For a curved surface, setting $C_x=C_z=0$ in the consistency condition gives
\[ 0 =  Y_z \partial_x C_y - Y_x \partial_z C_y = \nabla_{x,z}Y \cdot \vectwo{-\partial_z C_y}{\partial_x C_y} = \nabla_{x,z}Y \cdot R_{90^\circ}\nabla_{x,z}C_y. \]
This says the gradient of $Y$ in the mid-plane is perpendicular to the gradient of $C_y$ rotated by 90 degrees; in other words, the gradients of $Y$ and $C_y$ are parallel (or anti-parallel).  This may be achieved in simple cases such as when the gradients both point towards or away from the origin: that is, when $Y$ and $C_y$ are both functions of radius.  It also happens whenever the functions share the same contours, for example when $C_y=f(Y)$ for some function $f$.

\subsection{Functions of Radius\label{-exai-rad}}
It may be that $Y$ or $\C$ or some of its elements are only functions of radius $r=\sqrt{x^2+z^2}$ from the $y$ axis.  If the scalar involved can be expressed as $f(r^2)$, then the first partial derivative is $\partial_xf(r^2) = 2xf'(r^2)$, suggesting the general form 
\[ \partial_x^n f(r^2) = \sum_i p_{ni}(x) f^{(i)}(r^2) \]
for some polynomials $p_{ni}$.  The general term differentiates as
\[ \partial_x\left(p_{ni}(x) f^{(i)}(r^2)\right) = p'_{ni}(x) f^{(i)}(r^2) + 2x p_{ni}(x) f^{(i+1)}(r^2), \]
so the recurrence relation for these polynomials is:
\[ p_{n+1,i}(x) = p'_{ni}(x) + 2x p_{n,i-1}(x) \]
with the initial conditions $p_{00}=1$, $p_{0i}=0$ for $i>0$ and trivial boundary condition $p_{ni}=0$ for $i<0$.

It turns out the polynomials are all actually monomials of the form $p_{ni}(x) = a_{ni}x^{2i-n}$.  This is made clear by substitution into the recurrence relation:
\begin{eqnarray*}
 a_{n+1,i}x^{2i-n-1} &=& (2i-n)a_{ni}x^{2i-n-1} + 2x a_{n,i-1}x^{2i-n-2} \\
\Rightarrow\qquad
a_{n+1,i} &=& (2i-n)a_{ni} + 2a_{n,i-1}.
\end{eqnarray*}
The array of coefficients $a_{ni}$ is now easy to compute and taking account of the areas where they are zero,
\newcommand{\ceil}[1]{\left\lceil{#1}\right\rceil }
\[ \partial_x^n f(r^2) = \sum_{i=\ceil{n/2}}^n a_{ni}x^{2i-n} f^{(i)}(r^2). \]

Since $\partial_z r^2 = 2z$ analogously to $\partial_x r^2 = 2x$, a similar formula holds for $z$ derivatives of an arbitrary function of $r^2$:
\[ \partial_z^n f(r^2) = \sum_{i=\ceil{n/2}}^n a_{ni}z^{2i-n} f^{(i)}(r^2). \]
Using both of these formulae enables arbitrary mixed derivatives of $f(r^2)$ to be calculated:
\begin{eqnarray*}
f_{nm} \equiv \partial_x^n \partial_z^m f(r^2)
&=& \partial_z^m \sum_{i=\ceil{n/2}}^n a_{ni}x^{2i-n} f^{(i)}(r^2) \\
&=& \sum_{i=\ceil{n/2}}^n a_{ni}x^{2i-n} \partial_z^m f^{(i)}(r^2) \\
&=& \sum_{i=\ceil{n/2}}^n a_{ni}x^{2i-n} \sum_{j=\ceil{m/2}}^m a_{mj}z^{2j-m} f^{(i+j)}(r^2). 
\end{eqnarray*}

\subsection{Functions with Azimuthal Flutter\label{-exai-azi}}
Often a function of the form $f(r^2)$ is not complex enough but the field model is still motivated by polar coordinates.  Simple models of cyclotrons with sectors may be obtained using fields of the form $C_y=f(r^2)g(\theta)$.  The first step in evaluating the repeated derivatives of this is binomial expansion of the product rule
\begin{eqnarray*}
\phi_{nm} \equiv \partial_x^n \partial_z^m \left(f(r^2)g(\theta)\right)
&=& \partial_z^m \sum_{i=0}^n {n \choose i} \partial_x^i f(r^2)\partial_x^{n-i} g(\theta) \\
&=& \sum_{i=0}^n {n \choose i} \partial_z^m\left(\partial_x^i f(r^2)\partial_x^{n-i} g(\theta)\right) \\
&=& \sum_{i=0}^n {n \choose i} \sum_{j=0}^m {m \choose j} \partial_x^i\partial_z^j f(r^2)\partial_x^{n-i}\partial_z^{m-j} g(\theta) \\
&=& \sum_{i=0}^n \sum_{j=0}^m {n \choose i} {m \choose j} f_{ij}\partial_x^{n-i}\partial_z^{m-j} g(\theta),
\end{eqnarray*}
where $f_{ij}$ is the result from the previous section.  The next step is to find what partial derivatives of $g(\theta)$ look like:
\[ \partial_x g(\theta) = g'(\theta)\pdd{\theta}{x} = \frac{z}{r^2}g'(\theta),
\qquad\qquad
\partial_z g(\theta) = g'(\theta)\pdd{\theta}{z} = \frac{-x}{r^2}g'(\theta). \]
This suggests a general term $T=x^az^b(r^2)^{-c}g^{(d)}(\theta)$, which differentiates as
\begin{eqnarray*}
\pdd{T}{x} &=& ax^{a-1}z^b(r^2)^{-c}g^{(d)}(\theta) + x^az^b(-c)(r^2)^{-c-1}(2x)g^{(d)}(\theta) + x^az^b(r^2)^{-c}\frac{z}{r^2}g^{(d+1)}(\theta) \\ 
&=& ax^{a-1}z^b(r^2)^{-c}g^{(d)}(\theta) - 2cx^{a+1}z^b(r^2)^{-c-1}g^{(d)}(\theta) + x^az^{b+1}(r^2)^{-c-1}g^{(d+1)}(\theta), \\
\pdd{T}{z} &=& bx^az^{b-1}(r^2)^{-c}g^{(d)}(\theta) - 2cx^az^{b+1}(r^2)^{-c-1}g^{(d)}(\theta) - x^{a+1}z^b(r^2)^{-c-1}g^{(d+1)}(\theta). 
\end{eqnarray*}
Computational algebra can now express $\partial_x^{n-i}\partial_z^{m-j} g(\theta)$ as the sum of terms of the above form with $1\le d\le n-i+m-j$ (or $d=0$ in the trivial no-derivative case).  The algebraic forms of the combined derivatives of $g$ only need to be calculated and stored once.  The order of evaluation per-point is: $f^{(n)}$ giving $f_{ij}$ as in the previous section, $g^{(n)}$ giving $\partial_x^i\partial_z^j g$ and then using both of these to get $\phi_{ij}$, which is an input to the calculations in section \ref{-dispy-eval-point}.

\subsubsection{Spiral Sectors\label{-exai-azi-spiral}}
\newcommand{\dr}{\mathrm dr}
Instead of the radial lines given by $\theta = \mbox{constant}$, some machines use spiralled azimuthal variation where the constant lines of a new variable $\eta$ make an `edge angle' $\epsilon$ with the radial lines.  This means that in a radial increment $\dr$, the line of constant $\eta$ must move a distance $(\tan\epsilon)\dr$ in the positive $\theta$ direction: that is, an angle of $\frac{\tan\epsilon}{r}\dr$.  This can be achieved if $\pdd{\eta}{r} = -\frac{\tan\epsilon}{r}$, where the partial derivative keeps $\theta$ constant, assuming $\eta$ has the normal angular dependence ($\pdd{\eta}{\theta}=1$).  An explicit solution to these conditions is $\eta = \theta - (\tan\epsilon)\ln r$.

Defining $\tau = \tan\epsilon$ for convenience, $\eta = \theta - \tau\ln r$.  Evaluating mixed partial derivatives of $g(\eta)$ will proceed similarly to the last section except for the factors
\[ \pdd{\eta}{x} = \pdd{\eta}{\theta}\pdd{\theta}{x} + \pdd{\eta}{r}\pdd{r}{x} = \frac{z}{r^2} -\frac{\tau}{r}\frac{x}{r} = \frac{z-\tau x}{r^2} \]
and
\[ \pdd{\eta}{z} = \pdd{\eta}{\theta}\pdd{\theta}{z} + \pdd{\eta}{r}\pdd{r}{z} = \frac{-x}{r^2} -\frac{\tau}{r}\frac{z}{r} = \frac{-x-\tau z}{r^2}. \]
The general term $T$ will have the same form as before with $\theta$ replaced by $\eta$ and the derivatives will differ by $-\tau x^{a+1}z^b(r^2)^{-c-1}g^{(d+1)}(\eta)$ added to $\pdd{T}{x}$ and $-\tau  x^az^{b+1}(r^2)^{-c-1}g^{(d+1)}(\eta)$ added to $\pdd{T}{z}$.

\subsection{General Maxwellian Consistency}
If a desired $Y$ displacement function and vertical field $C_y$ do not satisfy the conditions of the simple case in section \ref{-exai-vert}, can the situation be salvaged?  One way would be to add horizontal field components $C_x$ and $C_z$ such that the consistency condition 
\[ \partial_z C_x - \partial_x C_z = Y_z \partial_x C_y - Y_x \partial_z C_y \equiv F \]
holds, noting that the right-hand side is now completely known.  In fact, this is easy to do because there are two arbitrary functions $C_x$ and $C_z$ but only one constraint.  That means fairly trivial solutions exist, such as:
\[ C_x=0; \qquad C_z=\int_0^x F(X,z) \mathrm dX. \]
However, it is desirable that these additional components exhibit the same symmetry or at least follow the same general shape as the main $C_y$ component, since the components are all produced from the same magnet (typically the extra components are parts of the fringe field, where $C_y$ changes).  A similar technique was used in \cite{VFFAGedges}.

\subsubsection{Radial Integration\label{-exai-max-rad}}
For the case where $C_y$ is produced from radial sector magnets (section \ref{-exai-azi}), the natural direction of integration is radially outwards.  As the consistency condition only involves the derivatives of $C_x$ and $C_z$, a valid and symmetrical initial condition is to set $C_x=C_z=0$ at $r=0$.

Changing to cylindrical polar coordinates, $\partial_z C_x - \partial_x C_z = (\nabla\times\C)_y = \frac1r C_\theta + \partial_rC_\theta - \frac1r\partial_\theta C_r$.  If the integration is to proceed radially, it makes sense to set $C_r=0$; this also removes any transverse field that could deflect a circulating beam travelling in the $\theta$ direction.  For a constant value of $\theta$, $C_\theta$ is now the solution of the ordinary differential equation
\[ \frac1r C_\theta + \partial_rC_\theta = F \]
with $C_\theta=0$ at $r=0$.  Multiplying both sides by $r$ gives
\[ C_\theta + r\partial_rC_\theta = \partial_r(rC_\theta) = rF \]
and therefore
\[ rC_\theta = \int_0^r RF \mathrm dR
\qquad\Rightarrow\qquad
C_\theta(r,\theta) = \frac1r \int_0^r RF(R,\theta) \mathrm dR.\]
This can be converted back to Cartesian coordinates with $C_x=zC_\theta/r$ and $C_z=-xC_\theta/r$.

\subsubsection{Radial Sectors with $Y$ a Function of Radius\label{-exai-max-sectors}}
\newcommand{\avg}[1]{\left\langle{#1}\right\rangle}
Further simplification can take place if the form $C_y=f(r^2)g(\theta)$ from section \ref{-exai-azi} is combined with a radial variation of height $Y=Y(r^2)$ as per section \ref{-exai-rad}.  These functions have various partial derivatives $Y_x=2xY'$, $Y_z=2zY'$ and 
\[ \partial_xC_y = 2xf'g + f\frac{z}{r^2}g', \qquad \partial_zC_y = 2zf'g + f\frac{-x}{r^2}g', \]
meaning that
\begin{eqnarray*}
F &=& 2zY'\left(2xf'g + f\frac{z}{r^2}g'\right) - 2xY'\left(2zf'g + f\frac{-x}{r^2}g'\right) \\
&=& 2zY'f\frac{z}{r^2}g' - 2xY'f\frac{-x}{r^2}g' \\
&=& 2Y'fg'\frac{z^2+x^2}{r^2} \\
&=& 2Y'fg'.
\end{eqnarray*}
Now using this in the previous section's result gives
\[ C_\theta(r,\theta) = \frac1r \int_0^r R2Y'(R^2)f(R^2)g'(\theta) \mathrm dR
 = \frac{2g'(\theta)}r \int_0^r RY'(R^2)f(R^2) \mathrm dR.\]
Notice that $2RY'(R^2) = \dd{R^2}{R}\dd{Y}{R^2} = \dd{Y}{R}$ so if the functions $Y$ and $f$ are re-expressed as functions of $r$ instead of $r^2$ (just for the following one formula), then
\[ C_\theta(r,\theta) = \frac{g'(\theta)}r \int_0^r Y'(R)f(R) \mathrm dR = \avg{Y'f}_{[0,r]}g'(\theta).\]

The additional components ($C_\theta$ or $C_x,C_z$) will also have to be differentiated repeatedly to find the off-surface field, so it is helpful to express $C_\theta=f_\theta(r^2)g_\theta(\theta)$, where $f_\theta(r^2) = \avg{2RY'(R^2)f(R^2)}_{[0,r]}$ and $g_\theta=g'$.  The Cartesian components also preserve this form via $C_x=C_\theta\cos\theta$ and $C_z=-C_\theta\sin\theta$, so $f_x=f_z=f_\theta$, $g_x(\theta)=g'(\theta)\cos\theta$ and $g_z(\theta)=-g'(\theta)\sin\theta$.

\subsubsection{Spiral Integration}
Starting with the polar form of the consistency condition $\frac1r C_\theta + \partial_rC_\theta - \frac1r\partial_\theta C_r=F$, consider working in coordinates of $(r,\eta)$ instead of $(r,\theta)$.  The partial derivatives in the $(r,\eta)$ system are $\partial_\eta$, which changes $\eta$ while keeping $r$ constant and $\partial_{r|\eta}$, which changes $r$ while keeping $\eta$ constant, differing from $\partial_r$ that keeps $\theta$ constant.  Because $\pdd{\eta}{\theta}=1$ and they both keep $r$ constant, $\partial_\eta=\partial_\theta$.  However, to keep $\eta$ constant, $\partial_{r|\eta} = \partial_r - \pdd{\eta}{r}\partial_\eta = \partial_r + \frac{\tau}{r}\partial_\theta$.  Substituting this into the consistency condition gives:
\begin{eqnarray*}
F &=& \frac1r C_\theta + \left(\partial_{r|\eta} - \frac{\tau}{r}\partial_\theta\right)C_\theta - \frac1r\partial_\theta C_r \\
&=& \frac1r C_\theta + \partial_{r|\eta}C_\theta - \frac1r\partial_\theta(C_r+\tau C_\theta). 
\end{eqnarray*}
Now the logical choice for integration is to set $C_r+\tau C_\theta=0$ by choosing $C_r=-\tau C_\theta$ and then solving
\[  \frac1r C_\theta + \partial_{r|\eta}C_\theta = F \]
analogously to section \ref{-exai-max-rad}, yielding
\[ C_\theta(r,\eta) = \frac1r \int_0^r RF(R,\eta) \mathrm dR. \]
The conversions back to Cartesian coordinates are now $C_x=(zC_\theta+xC_r)/r=(z-\tau x)C_\theta/r$ and $C_z=(-xC_\theta+zC_r)/r=(-x-\tau z)C_\theta/r$.

\subsubsection{Spiral Sectors with $Y$ a Function of Radius}
Following section \ref{-exai-max-sectors}, let $C_y=f(r^2)g(\eta)$ and $Y=Y(r^2)$; partial derivatives of $C_y$ differ:
\[ \partial_xC_y = 2xf'g + f\frac{z-\tau x}{r^2}g', \qquad \partial_zC_y = 2zf'g + f\frac{-x-\tau z}{r^2}g', \]
giving
\begin{eqnarray*}
F &=& 2zY'\left(2xf'g + f\frac{z-\tau x}{r^2}g'\right) - 2xY'\left(2zf'g + f\frac{-x-\tau z}{r^2}g'\right) \\
&=& 2zY'f\frac{z-\tau x}{r^2}g' - 2xY'f\frac{-x-\tau z}{r^2}g' \\
&=& 2Y'fg'\frac{z^2-\tau xz+x^2+\tau xz}{r^2} \\
&=& 2Y'fg'.
\end{eqnarray*}
This is exactly the same as last time so all the formulae for $C_\theta$ are the same with $\theta$ dependencies replaced by $\eta$.  In particular, $C_\theta=f_\theta(r^2)g_\theta(\eta)$ where $f_\theta(r^2) = \frac{2}r \int_0^r RY'(R^2)f(R^2) \mathrm dR$ and $g_\theta=g'$.

The Cartesian components can be expressed as $C_x=C_\theta(\cos\theta - \tau\sin\theta)$, $C_z=C_\theta(-\sin\theta-\tau\cos\theta)$.  This looks problematic because $\theta=\eta+\tau\ln r$ but the trigonometric angle-sum identities give:
\begin{eqnarray*}
C_x &=& C_\theta(\cos\eta\cos(\tau\ln r) - \sin\eta\sin(\tau\ln r) - \tau\sin\eta\cos(\tau\ln r) - \tau\cos\eta\sin(\tau\ln r)) \\
&=& C_\theta(\cos(\tau\ln r) - \tau\sin(\tau\ln r))\cos\eta + C_\theta(-\sin(\tau\ln r)-\tau\cos(\tau\ln r))\sin\eta \\
&=& C_\theta f_c(r)\cos\eta + C_\theta f_s(r)\sin\eta \\
C_z &=& C_\theta(-\sin\eta\cos(\tau\ln r)-\cos\eta\sin(\tau\ln r) - \tau\cos\eta\cos(\tau\ln r) + \tau\sin\eta\sin(\tau\ln r)) \\
&=& C_\theta(-\sin(\tau\ln r)-\tau\cos(\tau\ln r))\cos\eta + C_\theta(-\cos(\tau\ln r)+\tau\sin(\tau\ln r))\sin\eta \\
&=& C_\theta f_s(r)\cos\eta - C_\theta f_c(r)\sin\eta.
\end{eqnarray*}
This means that for example $C_x=f_{x1}(r^2)g_{x1}(\eta)+f_{x2}(r^2)g_{x2}(\eta)$ where $g_{x1}(\eta)=g'(\eta)\cos\eta$, $f_{x1}(r^2)=f_\theta(r^2)f_c(r)$ and similarly for the other term with $g_{x2}(\eta)=g'(\eta)\sin\eta$.  Differentation is a linear operation so repeated derivatives of $C_x$ and $C_z$ can be found with two applications of the algorithm in sections \ref{-exai-azi} and \ref{-exai-azi-spiral}.

A cleaner way of deriving this result is to consider the complex function $Z=C_z+\ii C_x=C_\theta(\ii\e^{\ii\theta}-\tau\e^{\ii\theta})=C_\theta(\ii-\tau)\e^{\ii(\eta+\tau\ln r)}=f_\theta(r^2)(\ii-\tau)\e^{\ii\tau\ln r}g_\theta(\eta)e^{\ii\eta}$.  $Z$ is a complex product so each of its components is the sum of two real products.

\subsubsection{Evaluation Strategy for Maxwellian Spiral Sectors\label{-exai-max-eval}}
The initial condition vector including fringe fields breaks down into product components as
\begin{eqnarray*}
C_x &=& f_1(r^2)g_1(\eta)+f_2(r^2)g_2(\eta) \\
C_y &=& f(r^2)g(\eta) \\
C_z &=& f_2(r^2)g_1(\eta)-f_1(r^2)g_2(\eta),
\end{eqnarray*}
bearing in mind the algorithm in section \ref{-exai-azi} will want repeated derivatives of all the single-parameter functions.  It is assumed that repeated derivatives of $f$ and $g$ are provided by the user (as well as those of $Y(r^2)$).  The new single parameter functions are
\begin{eqnarray*}
f_1(r^2)=f_\theta(r^2)f_c(r^2) &\qquad& g_1(\eta)=g'(\eta)\cos\eta \\
f_2(r^2)=f_\theta(r^2)f_s(r^2) &\qquad& g_2(\eta)=g'(\eta)\sin\eta
\end{eqnarray*}
where
\begin{eqnarray*}
f_\theta(r^2) &=& \frac{1}r \int_0^r 2RY'(R^2)f(R^2) \mathrm dR \\
f_c(r^2) &=& \cos(0.5\tau\ln r^2) - \tau\sin(0.5\tau\ln r^2) \\
f_s(r^2) &=& -\sin(0.5\tau\ln r^2)-\tau\cos(0.5\tau\ln r^2).
\end{eqnarray*}
Functions that are products can be dealt with by calculating the derivatives of the factors first and then using the rule $(fg)^{(n)}=\sum_{i=0}^n {n \choose i} f^{(i)}g^{(n-i)}$.

Assuming numerical integration is not used, the user will also have to supply a function to calculate $I(r^2)=\int_0^r 2RY'(R^2)f(R^2) \mathrm dR$, although not its derivatives since
\[ I'(r^2) = \dd{r}{r^2}\dd{I}{r} = \frac1{2r}2rY'(r^2)f(r^2) = Y'(r^2)f(r^2). \]

The trigonometric functions of logarithms have repeated derivatives related by the formula
\[ \dd{}{x}\left(\frac1{x^n}\left(a\cos(k\ln x)+b\sin(k\ln x)\right)\right) \]
\[ = \frac{-n}{x^{n+1}}\left(a\cos(k\ln x)+b\sin(k\ln x)\right)+\frac1{x^n}\left(-a\sin(k\ln x)\frac{k}{x}+b\cos(k\ln x)\frac{k}{x}\right) \]
\[ = \frac1{x^{n+1}}\left((-na+kb)\cos(k\ln x)+(-nb-ka)\sin(k\ln x)\right). \]
Writing $z=a+\ii b$, the `new' $z$ can be written $(-na+kb)+\ii(-nb-ka) = (-n-\ii k)z$.  So if 
\begin{eqnarray*}
f_c^{(n)}(r^2) &=& \frac1{r^{2n}}\left(\mathrm{Re\,}z_n\cos(0.5\tau\ln r^2) + \mathrm{Im\,}z_n\sin(0.5\tau\ln r^2)\right) \\
f_s^{(n)}(r^2) &=& \frac1{r^{2n}}\left(\mathrm{Re\,}w_n\cos(0.5\tau\ln r^2) + \mathrm{Im\,}w_n\sin(0.5\tau\ln r^2)\right),
\end{eqnarray*}
then the recurrence starts with $z_0=1-\ii\tau$ and $w_0=-\tau-\ii=-\ii z_0$, continuing with $z_{n+1}=(-n-0.5\ii\tau)z_n$, where by linearity $w_n=-\ii z_n$.


\subsection{The Isochronous Field}
In an isochronous machine, path length must be proportional to velocity, so the mean orbit radius must satisfy $r=\beta R$ for some $R$.  The average field required for the orbit to close is $\avg{B_y}=p/qr=m\beta\gamma c/(q\beta R)=\gamma B_0$ where $B_0=mc/qR$.  This is familiar as the field relation for an isochronous cyclotron but it will also work for other machines whose mean field is a function of radius.

In the notation of the previous sections, we can choose
\[ f(r^2) = \gamma B_0 = \frac{B_0}{\sqrt{1-\beta^2}} = \frac{B_0}{\sqrt{1-\frac{1}{R^2}r^2}}, \]
which will be multiplied by some azimuthal variation $g(\theta)$ or $g(\eta)$ to give the field $C_y$.  The functions in the previous sections need $f$'s repeated derivatives with respect to $r^2$.  These are
\[ f^{(n)}(r^2) = B_0\frac{(2n)!}{4^nn!R^{2n}}\left(1-\frac{1}{R^2}r^2\right)^{-\frac12-n}. \]

\subsubsection{Integrals}
The algorithm in section \ref{-exai-max-eval} requires the integral of $Y'(r^2)f(r^2)$ with respect to $r^2$.  If $Y(r^2)$ can be approximated by a polynomial, it suffices to find integrals $I_n$ such that $\dd{I_n}{r^2}=r^{2n}f(r^2)$.  Consider the derivative
\begin{eqnarray*}
 \dd{}{r^2}\left(r^{2n}\sqrt{1-\frac{1}{R^2}r^2}\right) &=& 
nr^{2n-2}\sqrt{1-\frac{1}{R^2}r^2} +
r^{2n}\frac{-\frac{1}{R^2}}{2\sqrt{1-\frac{1}{R^2}r^2}} \\
&=& \frac{nr^{2n-2}\left(1-\frac{1}{R^2}r^2\right)-\frac1{2R^2}r^{2n}}{\sqrt{1-\frac{1}{R^2}r^2}} \\
&=& \frac{nr^{2n-2}-(n+\frac12)\frac1{R^2}r^{2n}}{\sqrt{1-\frac{1}{R^2}r^2}} \\
&=& \left(nr^{2n-2}-\left(n+\frac12\right)\frac1{R^2}r^{2n}\right)\frac{f(r^2)}{B_0} \\
&=& \frac1{B_0}\left(n\dd{I_{n-1}}{r^2}-\left(n+\frac12\right)\frac1{R^2}\dd{I_n}{r^2}\right) \\
\Rightarrow \qquad r^{2n}\sqrt{1-\frac{1}{R^2}r^2} &=& \frac1{B_0}\left(nI_{n-1}-\left(n+\frac12\right)\frac1{R^2}I_n\right) + \mbox{const.}.
\end{eqnarray*}
Setting $n=0$ gives
\[ \sqrt{1-\frac{1}{R^2}r^2} = \frac1{B_0}\left(-\frac1{2R^2}I_0\right) + \mbox{const.}
\qquad \Rightarrow \qquad 
I_0 = -2B_0R^2\sqrt{1-\frac{1}{R^2}r^2} + \mbox{const.}, \]
where the constant is set by noting the definite integral starts at $r^2=0$ so requires $I_0(0)=0$; thus
\[ I_0 = -2B_0R^2\sqrt{1-\frac{1}{R^2}r^2} + 2B_0R^2. \]
The cases for $n>0$ form a recurrence
\begin{eqnarray*}
B_0r^{2n}\sqrt{1-\frac{1}{R^2}r^2} - nI_{n-1} &=& -\left(n+\frac12\right)\frac1{R^2}I_n \\
\Rightarrow \qquad I_n &=& \frac{R^2}{n+\frac12}\left(nI_{n-1}-B_0r^{2n}\sqrt{1-\frac{1}{R^2}r^2}\right),
\end{eqnarray*}
which satisfy $I_n(0)=0$ without the addition of a constant, by induction.

\subsection{Flutter Models}
A simple periodic variation with $n$ sectors can be produced with $g(\theta)=1+a\sin n\theta$.  This is easy to differentiate repeatedly, however it is not very realistic as magnets commonly have harder edges and do not fade out gradually like a sine wave.

Another simple function worth considering is $g(\theta)=1+(\sin n\theta)^{1/(2m+1)}$ for integer $m\ge0$.  This has mean value $1$ and increasing $m$ makes the magnet edges sharper and the field within the sectors flatter.  It is limited to sectors and drifts of equal length and its repeated derivatives do not have simple form, though.

It is possible to construct a function with arbitrary fringe field lengths by summing scaled copies of the $\tanh$ function and this may be differentiated repeatedly using the method in \cite{VFFAGedges} section 5.  Care must be taken that the function wraps properly from $\theta=2\pi$ to $0$, although for short fringe lengths all but a few of the $\tanh$ functions may be approximated by a constant at any given point.

\subsubsection{Flutter Varying with Radius}
The field form $C_y=f(r^2)g(\theta)$ used in most of the preceding sections is not complex enough to allow the proportion of field flutter to vary as a function of radius.  One way of incorporating this variation is 
\[ C_y=f(r^2)(1+h(r^2)g(\theta)) = f(r^2) + f(r^2)h(r^2)g(\theta), \]
where $g(\theta)$ now oscillates about zero rather than one and $h(r^2)$ controls the amplitude of those oscillations.  Notice that this is the sum of two terms, both of which have the form $f(r^2)g(\theta)$: namely, $f.1$ and $fh.g$.  Because Maxwell's equations are linear, the fields corresponding to each of the two terms can be calculated using the method in section \ref{-dispy} and then added together.

% Better flutter functions (not just sin)?

% Generalised operations on sets of derivatives?  (+ is trivial, * already written out.  what about f(g(x,z))?)

\begin{thebibliography}{9}   % Use for  1-9  references, {99} for 10-99 references
\bibitem{VFFAGedges} \emph{Fringe Fields for VFFAG Magnets with Edge Angles}, S.J. Brooks note 2011-6, available from \url{http://stephenbrooks.org/ap/report/2011-6/skewVFFAGfringes.pdf}.
\end{thebibliography}

\end{document}

