MAGNETIC SHIELD

Abstract

A method of designing a magnetic shield comprising a structure enclosing a space, the structure comprising passive magnetic shielding material and a winding configured to produce a specified magnetic within the structure when current is passed through the winding, is disclosed. The method comprises determining an optimised configuration of the winding accounting for the presence of the passive magnetic shield material by implementing one or more boundary conditions at the surface of the passive magnetic shielding material.

Description

FIELD OF THE INVENTION

The present invention relates to a magnetic shield, and methods for designing a magnetic shield.

BACKGROUND

Regions of space in which magnetic fields are cancelled so as to produce a substantially magnetic field free environment are required for a variety of different purposes [1], for example laboratory experiments, hospital testing etc. Typically, this is achieved by placing a large amount of high magnetic permeability material (i.e., passive magnetic shielding material) around the region in which magnetic fields are to be cancelled. The amount of high magnetic permeability material can be partially reduced by using an active cancellation system (for example, current carrying coils). However, the presence of high magnetic permeability material warps the magnetic fields produced by the active cancellation coil system.

Optimising active field cancellation (or generation) systems in the presence of passive magnetic shielding material (e.g., high magnetic permeability material) is a long-standing challenge in electromagnetism. It is required for multiple applications, for example biomedical imaging, quantum sensing of gravity, geophysical surveys and underground mapping, and noise suppression in precision measurement. However, optimisation of arbitrary active field-generating coils, or conductor networks, in the presence of passive magnetic shielding material (for example, Mu-metal®) has not yet been implemented.

SUMMARY

According to a first aspect, there is provided a method of designing a magnetic shield (or magnetic field generating system) comprising a structure enclosing a space, the structure comprising passive magnetic shielding material, and a winding configured to produce a specified magnetic field within the structure when current is passed through the winding. The method may comprise determining an optimised configuration of the winding accounting for the presence of the passive magnetic shielding material by implementing one or more boundary conditions at the surface of the passive magnetic shielding material. The method may be computer-implemented.

Warping of a desired magnetic field to be produced by a winding may be avoided by accounting for the presence of a passive magnetic shielding material in the optimisation process. This may improve accuracy of magnetic field generation and/or magnetic field cancellation in such magnetic shields.

The one or more boundary conditions may require that the magnetic field (for example, one or more vector components of the magnetic field) produced by the winding is zero on a surface of the passive magnetic shielding material.

Accounting for the presence of the passive magnetic shielding material may comprise constructing a function, or discrete approximation thereof, for a geometry of the structure, that can be used to solve differential equations relating magnetic field to current density.

The function may be a Green's function subject to one or more Dirichlet boundary conditions. The method may further comprise implementing the method of mirror images.

Determining an optimised configuration of the winding comprises implementing an optimisation process. The optimisation process may be a least squares minimisation process. The method may comprise regularising the least squares minimisation with a penalty term. The penalty term may be or comprise at least one of a power consumption of the winding, a curvature of the winding configuration, a resistance, an inductance of the winding, a mass and/or volume of the winding, and an energy stored in the winding. The optimisation process may comprise solving one or more sets of simultaneous equations.

Determining an optimised configuration of the winding may comprise determining optimal surface currents on the structure required to produce the specified magnetic field within the structure. The surface currents may be defined by a stream function. Determining the optimal surface currents may comprise determining streamlines of the stream function. Determining streamlines of the stream function may comprise discretising the streamfunction to determine contours of the stream function.

Determining an optimised configuration of the winding may comprise defining at least one discrete winding element having at least one free parameter. Determining an optimised configuration of the winding may further comprise optimising at least one free parameter of each of the discrete winding elements to produce the specified magnetic field within the structure. The free parameters of each of the discrete winding elements may be collectively optimised to produce the specified magnetic field within the structure.

Each of the discrete winding elements may be defined by parametric curves, for example one or more of a straight line, an arc, a circular loop, or any other discrete configuration of wire. The at least one free parameter of each discrete winding element may be or comprise at least one of a size, a shape, a spacing, a placement and an orientation of the discrete winding element.

Determining an optimised configuration of the winding may comprise determining expressions for components of current density for each of the discrete winding elements. The method may further comprise substituting the expressions for components of current density into expressions for components of the magnetic field produced by an arbitrary winding configuration. The method may also comprise implementing an optimisation process to determine, for the specified field, optimal values relating to the components of current density. The optimal values may be or include an optimal free parameter, and/or an optimal relationship between free parameters, of the discrete winding elements. The optimal values may be or include an optimal relationship between radii and separations of the discrete winding elements.

The magnetic field may be one of a constant magnetic field, a linear gradient magnetic field or a higher-order magnetic field.

The structure may be a closed structure. The structure may be a hollow cylinder. The hollow cylinder may comprise an axial body and planar end surfaces. Separate boundary conditions may be implemented for an axial surface of the cylinder and planar end surfaces of the cylinder.

The method may further comprise making a winding corresponding to or approximating the optimised configuration of the winding.

According to a second aspect, there is provided a magnetic shield (or magnetic field generating system) designed using the method of the first aspect.

According to a third aspect, there is provided a magnetic shield (or magnetic field generating system). The magnetic shield may comprise a structure enclosing a space. The structure may comprise passive magnetic shielding material. The magnetic shield may also comprise a winding configured to produce a specified magnetic field within the structure when current is passed through the winding. A configuration of the winding may be optimised such that a field produced by the winding is more similar to the specified magnetic field when the winding is in the presence of the structure than when the winding is not in the presence of the structure.

A magnetic shield in accordance with the third aspect may have improved accuracy of magnetic field generation and/or magnetic field cancellation when compared to existing magnetic shields. A reduction in performance of the winding (for example, an error relative to the specified magnetic field) may be measurable when the winding is operated not in the presence of the structure comprising passive magnetic shielding material.

The structure may be a closed structure. The structure may be a hollow cylinder.

The winding may be disposed within the space enclosed by the structure. The winding may be disposed on an interior surface of the structure.

The winding, in the presence of the structure, may be configured to produce a uniform magnetic field over a volume of at least 1000 cm³ with a relative field variation of less than 1%.

BRIEF DESCRIPTION OF DRAWINGS

The invention will now be described by way of example with reference to the accompanying drawings in which:

FIG. 1 shows a method for designing a magnetic shield in accordance with the invention;

FIGS. 2A and 2B show other methods for designing a magnetic shield in accordance with the invention, using a Green's function approach;

FIG. 3 shows a conducting open cylindrical shell of length L and radius a about the origin;

FIG. 4 shows a passive magnetic shield cylinder of length {circumflex over (L)} and radius b with planar end caps located at z=±{circumflex over (L)}/2 enclosing an interior conducting open cylindrical shell of length L₁-L₂ and radius a;

FIGS. 5A to 5D show winding configurations optimised omitting interactions with passive magnetic shielding material;

FIGS. 5E to 5H show field homogeneity of magnetic fields produced by the windings shown in FIGS. 5A to 5D;

FIGS. 6A to 6D show winding configurations optimised accounting for interactions with passive magnetic shielding material in accordance with an embodiment of the invention;

FIGS. 6E and 6F show field homogeneity of magnetic fields produced by the windings shown in FIGS. 6A to 6D;

FIGS. 6G and 6H respectively show a flow of current in various regions of the winding shown in FIG. 6D, and how loops in the winding configuration of FIG. 6D are connected to produce the flow of current in FIG. 6G; and

FIGS. 7A to 7D show optimised winding configurations in accordance with an embodiment of the invention, and associated magnetic field strength distributions

FIG. 8 shows a hollow cylinder of high magnetic permeability of length L_s, radius p_s with planar end caps located at z=±L_s/2 enclosing a circular coil of radius ρ_c<ρ_s centred along the z-axis and lying in the z=z′ plane;

FIGS. 9A and 9B show planar winding configurations optimised accounting for interactions with passive magnetic shielding material in accordance with an embodiment of the invention;

FIGS. 9C and 9D show magnetic field and magnetic field homogeneity generated by the winding configurations shown in FIGS. 9A and 9B.

Features which are described in the context of separate aspects and embodiments of the invention may be used together and/or be interchangeable wherever possible. Similarly, where features are, for brevity, described in the context of a single embodiment, these may also be provided separately or in any suitable sub-combination. Features described in combination with the system may have corresponding features definable with respect to the method, and these embodiments are specifically envisaged.

DETAILED DESCRIPTION

FIG. 1 shows a method 10 for designing a magnetic shield, or other magnetic field-generating structure. The magnetic shield comprises a structure enclosing a space, the structure comprising a passive magnetic shielding material. The magnetic shield also comprises a winding configured to produce a magnetic field to cancel unwanted magnetic fields, or generate specified fields, within the structure when current is passed through the winding. The method 10, at step 12, comprises accounting for the presence of the passive magnetic shielding material. In some embodiments, accounting for the presence of the passive magnetic shielding material comprises implementing one or more boundary conditions relating to physical constraints on the magnetic field produced by the winding due to the passive magnetic shielding material. The one or more boundary conditions may require that one or more components of the magnetic field produced by the winding are zero, or close to zero if the full characteristics of the shielding material are included, on a surface of the passive magnetic shielding material. The method 10, at step 14, also comprises determining an optimised configuration of the winding by determining optimal surface currents on the structure required to produce a desired magnetic field within the structure. The term ‘passive magnetic shielding material’ used herein is intended to include all materials which are magnetisable (i.e., the material can support or alter a magnetic field within itself) in response to an applied magnetic field. It is this property which enables such materials to provide a shielding effect from external magnetic fields. For example, a passive magnetic shielding material may have a relative permeability of substantially 1000 or greater.

One approach for designing a magnetic shield, or field-generating structure, as described with respect to the method outlined in FIG. 1 (described above) is the method 20 shown in FIG. 2A.

At step 25, the method 20 comprises constructing a function (for example a Green's function), or discrete approximation thereof, for a geometry of the structure of the magnetic shield, that can be used to solve differential equations relating magnetic field and current density.

At step 30, the method 20 comprises modifying the function (for example, the Green's function) using, or subjecting the function to, one or more boundary conditions. The one or more boundary conditions represent physical constraints relating to the magnetic field intended to be produced by the optimised winding configuration. The one or more boundary conditions may be Dirichlet boundary conditions. The one or more boundary conditions may require that one or more components of the magnetic field produced by the winding is zero on a surface of the passive magnetic shielding material.

At step 35, the method 20 comprises implementing the one or more boundary conditions used to modify the function (for example, the Green's function). For example, the method of mirror images may be used to implement the one or more boundary conditions used to modify the Green's function. The method of mirror images may be used to determine expressions for components of the magnetic field produced by an arbitrary winding configuration.

Steps 25 to 35 of the method 20 may account for the presence of the passive magnetic shielding material in determining an optimised configuration of the winding.

At step 40, the method 20 comprises defining surface currents using a streamfunction. The streamfunction represents the fact that current density is restricted to flow over a surface geometry, i.e., that surface currents are used to produce a desired magnetic field. The aim of the optimisation is to determine the streamfunction, and hence, the optimal winding configuration to give a desired magnetic field (for example, a winding configuration that most closely approximates the stream function).

At step 45, the method 20 comprises determining expressions for components of current density from the streamfunction. For example, a Fourier decomposition of the streamfunction may be used to determine expressions for components of current density.

At step 50, the method 20 comprises substituting the expressions for components of current density into the expressions for components of the magnetic field produced by an arbitrary winding configuration.

At step 55, the method 20 comprises implementing an optimisation process to determine, for a desired magnetic field, optimal values relating to the components of current density. The optimisation process may be a least squares minimisation. The least squares minimisation may be regularised using a penalty term. The penalty term may comprise at least one of a power consumption of the winding, a curvature of the winding configuration, an inductance and/or resistance of the winding, and an energy stored in the winding. The optimal values may be optimal Fourier coefficients (the Fourier coefficients introduced by a Fourier decomposition of the streamfunction).

At step 60, the method 20 comprises substituting the determined optimal values relating to components of current density into the stream function to determine streamlines of the stream function (and therefore, optimal surface currents) to produce the desired magnetic field. Determining streamlines of the stream function may comprise discretising the stream function to determine contours of the stream function (for example, such that a winding configuration can be used to approximate the contours of the streamfunction).

Steps 40 to 60 of the method 20 determine an optimised configuration of the winding by determining optimal surface currents on the structure required to produce a desired magnetic field within the structure.

Another approach for designing a magnetic shield, or field-generating structure, as described with respect to the method outlined in FIG. 1 (described above) is the method 120 shown in FIG. 2B. Steps 125 to 135 of the method 120 are substantially identical to the steps 25 to 35 of the method 20 described above. Steps 125 to 135 may account for the presence of the passive magnetic shielding material in determining an optimised configuration of the winding.

At step 140, the method 120 comprises defining or specifying at least one discrete winding element having one or more free parameters. As used herein, the term ‘discrete’ is intended to mean that each of the elements may be considered as a separate and distinct element from each of the other elements (although multiple discrete winding elements may be used to form a complete winding). The discrete winding elements may be defined by parametric curves such as straight lines, arcs or circular loops, or by any other discrete configuration of wire, to form a winding. The one or more free parameters may be or include at least one of a form (for example, a size and/or a shape) a spacing, a placement and an orientation of the discrete winding elements to determine an optimal winding configuration to produce a desired magnetic field. Spacing refers to the distance between associated discrete elements (for example, the distance between a pair of coils in either a symmetric or asymmetric arrangement, depending on the desired harmonic field) of the winding. Placement refers to the location of associated discrete elements (for example, a pair of coils in either a symmetric or asymmetric arrangement) relative to other discrete elements of the winding. The free parameters of the discrete elements may be collectively optimised to generate a desired magnetic field.

At step 145, the method 120 comprises determining expressions for components of current density for each of the at least one discrete winding elements. The expressions may relate components of current density to the one or more free parameters of each of the discrete winding elements.

At step 150, the method 120 comprises substituting the expressions for components of current density into expressions for components of a magnetic field produced by an arbitrary winding configuration.

At step 155, the method 120 comprises implementing an optimisation process to determine, for a desired magnetic field, one or more optimal values relating to the components of current density. The optimisation process may be implemented by solving simultaneous equations. The one or more optimal values may be or include an optimal free parameter, or an optimal relationship between free parameters of the discrete winding elements (for example, a ratio between radii and axial separation of circular loops).

Steps 140 to 155 of the method 120 determine an optimised configuration of the winding by optimising one or more free parameters of the discrete winding elements to produce a desired magnetic field within the structure.

In the embodiment discussed below, a winding configuration is optimised for a magnetic shield, or other magnetic field-generating system, comprising a closed cylinder of passive magnetic shielding material and a winding disposed in the interior of the closed cylinder. Boundary conditions ensure that the value of the magnetic field generated by the winding attain the correct physical limit at the boundaries of the closed cylinder.

The magnetic vector potential satisfies Poisson's non-homogeneous equation,

∇²A(r)=−μ₀J(r′)   (1.1)

where J(r′) is the current density. As presented in various papers and textbooks [2], the solution to this PDE is given in integral notation by

A(r)=μ₀∫_r′dr′G(r,r′)J(r′)   (1.2)

where, in general, a Green's function with Dirichlet boundary conditions satisfies

∇²G(r,r′)=−δ(r−r′)   (1.3)

Consider the case of the Green's function in cylindrical coordinates without passive magnetic shielding material. Since the Laplacian in cylindrical coordinates is separable, it is possible to construct a Green's function of the form

G(r,r′)=R(ρ,ρ′)Φ(ϕ,ϕ′)Z(z,z′)   (1.4)

The Laplacian operator in cylindrical coordinates is given by

$\begin{array}{cc}{\nabla }^{2}G=\frac{1}{\rho }\frac{\partial }{\partial \rho }\left(\rho \frac{\partial G}{\partial \rho }\right)+\frac{1}{{\rho }^2}\frac{{\partial }^{2}G}{\partial {\varphi }^2}+\frac{{\partial }^{2}G}{\partial z^2}& \left(1.5\right)\end{array}$

Inserting 1.4 into 1.5, it is found

$\begin{array}{cc}\frac{1}{\rho }\frac{R^{\prime }}{R}+\frac{R^{″}}{R}+\frac{1}{{\rho }^2}\frac{{\Phi }^{″}}{\Phi }+\frac{Z^{″}}{Z}=0,\rho \ne 0& \left(1.6\right)\end{array}$

The separation constants associated with Z may be chosen to be either ±k². Here, −k² is chosen in order to more easily apply boundary conditions when calculating the Green's function in the presence of passive magnetic shielding material. Separating 1.6, it is found that

Z″+k ² Z=0   (1.7)

Φ″+m²Φ=0   (1.8)

Substituting 1.7 and 1.8 into 1.6, it is found that

ρ²R″+ρR′−(k²ρ+m²)R=0   (1.9)

This results in 1.9 having the general solution

$\begin{array}{cc}{R}_m\left(\rho ,{\rho }^{\prime },k\right)=\{\begin{array}{cc}{A}_m\left({\rho }^{\prime },k\right)I_m\left(k\rho \right)+B_m\left({\rho }^{\prime },k\right)K_m\left(k\rho \right)& 0<\rho <{\rho }^{\prime }<\infty \\ C_m\left({\rho }^{\prime },k\right)I_m\left(k\rho \right)+D_m\left({\rho }^{\prime },k\right)K_m\left(k\rho \right)& 0<{\rho }^{\prime }<\rho <\infty \end{array}& \left(1.1\right)\end{array}$

where I_v(z) and K_v(z) are the modified Bessel functions of the first and second kind, respectively. In cylindrical coordinates we may write

$\begin{array}{cc}\begin{array}{c}\delta \left(r-r^{\prime }\right)=\frac{1}{\rho }\delta \left(\rho -{\rho }^{\prime }\right)\delta \left(\varphi -{\varphi }^{\prime }\right)\delta \left(z-z^{\prime }\right)\\ =\frac{1}{\rho }\delta \left(\rho -{\rho }^{\prime }\right)\frac{1}{4{\pi }^2}{\sum }_{m=-\infty }^{\infty }{e}^{\mathrm{im}\left(\varphi -{\varphi }^{\prime }\right)}{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{ik}\left(z-z^{\prime }\right)}\end{array}& \left(1.11\right)\end{array}$

using the Cauchy principle value, and satisfying the requirements for the azimuthal dependence to be single valued quantised and unbounded in z. In order to avoid future problems with negative values of k, 1.11 may be rewritten as

$\begin{array}{cc}\delta \left(r-r^{\prime }\right)=\frac{1}{\rho }\delta \left(\rho -{\rho }^{\prime }\right)\frac{1}{2{\pi }^2}{\sum }_{m=-\infty }^{\infty }{e}^{\mathrm{im}\left(\varphi -{\varphi }^{\prime }\right)}{\int }_0^{\infty }\mathrm{dk}\cos \left(k\left(z-z^{\prime }\right)\right)& \left(1.12\right)\end{array}$

As for the ρ dependence, it is required that G be finite at ρ,ρ′=0 and zero at ρ,ρ′=∞ for ρ′>ρ and ρ>ρ′ respectively. From these conditions, the coefficients B_m(ρ′,k)=0 and C_m(ρ′,k)=0 giving

$\begin{array}{cc}{R}_m\left(\rho ,{\rho }^{\prime },k\right)=\{\begin{array}{cc}{A}_m\left({\rho }^{\prime },k\right)I_m\left(k\rho \right)& 0<\rho <{\rho }^{\prime }<\infty \\ D_m\left({\rho }^{\prime },k\right)K_m\left(k\rho \right)& 0<{\rho }^{\prime }<\rho <\infty \end{array}& \left(1.13\right)\end{array}$

Writing the continuity theorem we have

R _m(ρ,ρ′,k)|_ρ=ρ′+=R_m(ρ,ρ′,k)|_ρ=ρ′−  (1.14)

giving

$\begin{array}{cc}{R}_m\left(\rho ,{\rho }^{\prime },k\right)=\{\begin{array}{cc}{\mathrm{AI}}_m\left(k\rho \right)K_m\left(k{\rho }^{\prime }\right)& 0<\rho <{\rho }^{\prime }<\infty \\ {\mathrm{AI}}_m\left(k{\rho }^{\prime }\right)K_m\left(k\rho \right)& 0<{\rho }^{\prime }<\rho <\infty \end{array}& \left(1.15\right)\end{array}$

The Green's function may now be written as

$\begin{array}{cc}G\left(r,r^{\prime }\right)=\frac{A}{2{\pi }^2}{\sum }_{m=-\infty }^{\infty }{e}^{\mathrm{im}\left(\varphi -{\varphi }^{\prime }\right)}{\int }_0^{\infty }\mathrm{dk}\cos \left(k\left(z-z^{\prime }\right)\right)I_m\left(k{\rho }_{<}\right)K_m\left(k{\rho }_{>}\right)& \left(1.16\right)\end{array}$

where ρ_> is the larger of ρ and ρ′ and ρ_< is the smaller of the two quantities. Substituting G into 1.5 and integrating over the infinitesimal interval (ρ′⁻, ρ′⁺) yields

$\begin{array}{cc}{{{\int }_{\rho ={\rho }^{\prime -}}^{\rho ={\rho }^{\prime +}}d\rho \frac{{\partial }^{2}G}{\partial z^2}+\frac{1}{{\rho }^2}\frac{{\partial }^{2}G}{\partial {\varphi }^2}+\frac{1}{\rho }\frac{\partial G}{\partial \rho }+\frac{{\partial }^{2}G}{\partial {\rho }^2}=\frac{\partial G}{\partial \rho }❘}_{\rho ={\rho }^{\prime +}}-\frac{\partial G}{\partial \rho }❘}_{\rho ={\rho }^{\prime -}}& \left(1.17\right)\end{array}$

which satisfies 1.3, giving

$\begin{array}{cc}\frac{A}{2{\pi }^2}\sum _{m=-\infty }^{\infty }{e}^{\mathrm{im}\left(\varphi -{\varphi }^{\prime }\right)}{\int }_0^{\infty }\mathrm{dk}\cos \left(k\left(z-z^{\prime }\right)\right)\left(I_m\left(k{\rho }^{\prime }\right)K_m^{\prime }\left(k{\rho }^{\prime }\right)-I_m^{\prime }\left(k{\rho }^{\prime }\right)K_m\left(k{\rho }^{\prime }\right)\right)=-{\int }_{\rho ={\rho }^{\prime -}}^{\rho ={\rho }^{\prime +}}d\rho \frac{-\delta \left(\rho -{\rho }^{\prime }\right)\delta \left(\varphi -{\varphi }^{\prime }\right)\delta \left(z-z^{\prime }\right)}{\rho }& \left(1.18\right)\end{array}$

Using the identity

$\begin{array}{cc}{I}_v\left(z\right)K_v^{\prime }\left(z\right)-I_v^{\prime }\left(z\right)K_v\left(z\right)=-\frac{1}{2}& \left(1.19\right)\end{array}$

where I′_m(z) and K′_m(z) are the derivatives of I_m(z) and K_m(z) with respect to z respectively, it then follows from 1.18 that A=1, resulting in the Green's function

$\begin{array}{cc}G\left(r,r^{\prime }\right)=\frac{1}{2{\pi }^2}{\sum }_{m=-\infty }^{\infty }{e}^{\mathrm{im}\left(\varphi -{\varphi }^{\prime }\right)}{\int }_0^{\infty }\mathrm{dk}\cos \left(k\left(z-z^{\prime }\right)\right)I_m\left({\mathrm{kp}}_{<}\right)K_m\left(k{\rho }_{>}\right)& \left(1.2\right)\end{array}$

which may be finally rewritten as

$\begin{array}{cc}G\left(r,r^{\prime }\right)=\frac{1}{4{\pi }^2}{\sum }_{m=-\infty }^{\infty }{e}^{\mathrm{im}\left(\varphi -{\varphi }^{\prime }\right)}{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{ik}\left(z-z^{\prime }\right)}{I}_m\left(❘k❘{\rho }_{<}\right)K_m\left(❘k❘{\rho }_{>}\right)& \left(1.21\right)\end{array}$

Now, generating optimal surface coils, on a finite open cylinder of length L and radius a about the origin (as shown in FIG. 3), in order to generate an arbitrary specified magnetic field in the interior of the cylinder is considered. The vector potential takes the form of 1.2, where

r=ρ cos ϕi+ρ sin ϕj+zk   (1.22 )

r′=a cos ϕ′i+a sin ϕ′j+z′k   (1.23)

Following the formulation of R. Turner and R. M. Bowley [3], the three components in cylindrical coordinates may be written as

A _ρ(r)=μ₀∫_r′dr′G(r,r′)J_ϕ(r′)sin(ϕ−ϕ′)   (1.24)

A _ϕ(r)=μ₀∫_r′dr′G(r,r′)J_ϕ(r′)cos(ϕ−ϕ′)   (1.25)

A _z(r)=μ₀∫_r′dr′G(r,r′)J_z(r′)   (1.26)

Through the substitution of 1.21 into 1.24, 1.25 and 1.26 the components of A on the interior of the cylinder are

$\begin{array}{cc}{A}_{\rho }\left(\rho ,\varphi ,z\right)=-\frac{i{\mu }_{0}a}{8{\pi }^2}\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{im}\varphi }{e}^{\mathrm{ikz}}\left[I_{m-1}\left(❘k❘\rho \right)K_{m-1}\left(❘k❘a\right)-I_{m+1}\left(❘k❘\rho \right)K_{m+1}\left(❘k❘a\right)\right]{\int }_0^{2\pi }d{\varphi }^{\prime }{e}^{-\mathrm{im}{\varphi }^{\prime }}{\int }_{-\infty }^{\infty }{\mathrm{dz}}^{\prime }{e}^{-{\mathrm{ikz}}^{\prime }}{J}_{\varphi }\left({\varphi }^{\prime },z^{\prime }\right)& \left(1.27\right)\end{array}$ $\begin{array}{cc}{A}_{\varphi }\left(\rho ,\varphi ,z\right)=\frac{{\mu }_{0}a}{8{\pi }^2}\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{im}\varphi }{e}^{\mathrm{ikz}}\left[I_{m-1}\left(❘k❘\rho \right)K_{m-1}\left(❘k❘a\right)+I_{m+1}\left(❘k❘\rho \right)K_{m+1}\left(❘k❘a\right)\right]{\int }_0^{2\pi }d{\varphi }^{\prime }{e}^{-\mathrm{im}{\varphi }^{\prime }}{\int }_{-\infty }^{\infty }{\mathrm{dz}}^{\prime }{e}^{-{\mathrm{kz}}^{\prime }}{J}_{\varphi }\left({\varphi }^{\prime },z^{\prime }\right)& \left(1.28\right)\end{array}$ $\begin{array}{cc}{A}_z\left(\rho ,\varphi ,z\right)=\frac{{\mu }_{0}a}{4{\pi }^2}\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{im}\varphi }{e}^{\mathrm{ikz}}{I}_m\left(❘k❘\rho \right)K_m\left(❘k❘a\right)\times {\int }_0^{2\pi }d{\varphi }^{\prime }{e}^{-\mathrm{im}{\varphi }^{\prime }}{\int }_{-\infty }^{\infty }{\mathrm{dz}}^{\prime }{e}^{-{\mathrm{ikz}}^{\prime }}{J}_z\left({\varphi }^{\prime },z^{\prime }\right)& \left(1.29\right)\end{array}$

At this point it becomes useful to define the Fourier transforms of the current densities,

$\begin{array}{cc}{J}_{\varphi }^m\left(k\right)=\frac{1}{2\pi }{\int }_0^{2\pi }d{\varphi }^{\prime }{e}^{-\mathrm{im}{\varphi }^{\prime }}{\int }_{-\infty }^{\infty }{\mathrm{dz}}^{\prime }{e}^{-{\mathrm{ikz}}^{\prime }}{J}_{\varphi }\left({\varphi }^{\prime },z^{\prime }\right)& \left(1.3\right)\end{array}$ $\begin{array}{cc}{J}_z^m\left(k\right)=\frac{1}{2\pi }{\int }_0^{2\pi }d{\varphi }^{\prime }{e}^{-\mathrm{im}{\varphi }^{\prime }}{\int }_{-\infty }^{\infty }{\mathrm{dz}}^{\prime }{e}^{-{\mathrm{ikz}}^{\prime }}{J}_z\left({\varphi }^{\prime },z^{\prime }\right)& \left(1.31\right)\end{array}$

such that their inverse transforms are given by

$\begin{array}{cc}{J}_{\varphi }\left({\varphi }^{\prime },z^{\prime }\right)=\frac{1}{2\pi }{\sum }_{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{im}{\varphi }^{\prime }}{e}^{{\mathrm{ikz}}^{\prime }}{J}_{\varphi }^m\left(k\right)& \left(1.32\right)\end{array}$ $\begin{array}{cc}{J}_Z\left({\varphi }^{\prime },z^{\prime }\right)=\frac{1}{2\pi }{\sum }_{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{im}{\varphi }^{\prime }}{e}^{{\mathrm{ikz}}^{\prime }}{J}_Z^m\left(k\right)& \left(1.33\right)\end{array}$

For future simplification a relation between J_z^m(k) and J_ϕ^m(k) will prove useful. As such, using the continuity equation, restricting the current to flow only on the surface of a cylinder of radius a, gives

$\begin{array}{cc}\frac{1}{a}\frac{\partial J_{\varphi }\left({\varphi }^{\prime },z^{\prime }\right)}{\partial {\varphi }^{\prime }}+\frac{\partial J_z\left({\varphi }^{\prime },z^{\prime }\right)}{\partial {\varphi }^{\prime }}=0,\mathrm{on}\rho =a& \left(1.34\right)\end{array}$

Substituting 1.32 and 1.33 into 1.34 results in

$\begin{array}{cc}\frac{i}{2\pi }{\sum }_{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathrm{dk}\left(\frac{m}{a}{J}_{\varphi }^m\left(k\right)+{\mathrm{kJ}}_z^m\left(k\right)\right)e^{\mathrm{im}{\varphi }^{\prime }}{e}^{{\mathrm{ikz}}^{\prime }}=0& \left(1.35\right)\end{array}$

yielding the relationship between J_z^m(k) and J_ϕ^m(k) to be given by

$\begin{array}{cc}{J}_z^m\left(k\right)=-\frac{m}{\mathrm{ka}}{J}_{\varphi }^m\left(k\right)& \left(1.36\right)\end{array}$

Through the use of 1.30 and 1.31, 1.27, 1.28 and 1.29 are simplified to

$\begin{array}{cc}{A}_{\rho }\left(\rho ,\varphi ,z\right)=-\frac{i{\mu }_{0}a}{4\pi }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{im}\varphi }{e}^{\mathrm{ikz}}\left[I_{m-1}\left(❘k❘\rho \right)K_{m-1}\left(❘k❘a\right)-I_{m+1}\left(❘k❘\rho \right)K_{m+1}\left(❘k❘a\right)\right]J_{\varphi }^m\left(k\right)& \left(1.37\right)\end{array}$ $\begin{array}{cc}{A}_{\varphi }\left(\rho ,\varphi ,z\right)=\frac{{\mu }_{0}a}{4\pi }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{im}\varphi }{e}^{\mathrm{ikz}}\left[I_{m-1}\left(❘k❘\rho \right)K_{m-1}\left(❘k❘a\right)-I_{m+1}\left(❘k❘\rho \right)K_{m+1}\left(❘k❘a\right)\right]J_{\varphi }^m\left(k\right)& \left(1.38\right)\end{array}$ $\begin{array}{cc}{A}_z\left(\rho ,\varphi ,z\right)=\frac{{\mu }_{0}a}{2\pi }{\sum }_{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dke}}^{\mathrm{im}\varphi }{e}^{\mathrm{ikz}}{I}_m\left(❘k❘\rho \right)K_m\left(❘k❘a\right)J_z^m\left(k\right)& \left(1.39\right)\end{array}$

From this the magnetic field generated by an arbitrary coil on the surface of a cylinder may be calculated, taking the curl of A,

$\begin{array}{cc}B=▽\wedge A=\left(\frac{1}{\rho }\frac{\partial A_z}{\partial \varphi }-\frac{\partial A_{\varphi }}{\partial z}\right)\hat{\rho }+\left(\frac{\partial A_{\rho }}{\partial z}-\frac{\partial A_z}{\partial \rho }\right)\hat{\varphi }+\frac{1}{\rho }\left(\frac{\partial }{\partial \rho }\left(\rho A_{\varphi }\right)-\frac{\partial A_{\rho }}{\partial \varphi }\right)\hat{z}& \left(1.4\right)\end{array}$

while using 1.36 and the recurrence relations for the modified Bessel functions

$\begin{array}{cc}{I}_{m\pm 1}\left(z\right)=I_m^{\prime }\left(z\right)\mp \frac{m}{z}{I}_m\left(z\right)& \left(1.41\right)\end{array}$ $\begin{array}{cc}{K}_{m\pm 1}\left(z\right)=-K_m^{\prime }\left(z\right)\pm \frac{m}{z}{K}_m\left(z\right)& \left(1.42\right)\end{array}$

where I′_m(z) and K′_m(z) are the derivatives of I_m(z) and K_m(z) with respect to z respectively, the magnetic field B components of the interior of the cylinder are given by

$\begin{array}{cc}{B}_{\rho }\left(\rho ,\varphi ,z\right)=\frac{i{\mu }_{0}a}{2\pi }{\sum }_{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathrm{dk}{\mathrm{ke}}^{\mathrm{im}\varphi }{e}^{\mathrm{ikz}}{I}_m^{\prime }\left(❘k❘\rho \right)K_m^{\prime }\left(❘k❘a\right)J_{\varphi }^m\left(k\right)& \left(1.43\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }\left(\rho ,\varphi ,z\right)=\frac{{\mu }_0{a}^2}{2\pi \rho }{\sum }_{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathrm{dk}❘k❘e^{\mathrm{im}\varphi }{e}^{\mathrm{ikz}}{I}_m\left(❘k❘\rho \right)K_m^{\prime }\left(❘k❘a\right)J_z^m\left(k\right)& \left(1.44\right)\end{array}$ $\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=-\frac{{\mu }_{0}a}{2\pi }{\sum }_{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathrm{dk}❘k❘e^{\mathrm{im}\varphi }{e}^{\mathrm{ikz}}{I}_m\left(❘k❘\rho \right)K_m^{\prime }\left(❘k❘a\right)J_{\varphi }^m\left(k\right)& \left(1.45\right)\end{array}$

in agreement with Turner [4]. As with the target field method [4], the current density required to generate a specified target field can now be explicitly calculated. However, since no explicit length of coil is specified in the target field method, the result is slightly unrealistically sized coils. Instead, a Fourier series decomposition of the current density may be employed. Since the current density has been restricted to flow over a surface geometry, through the use of the continuity equation 1.34, the current density may be defined by a single function—the streamfunction, ψ(ϕ′,z). It then follows that

$\begin{array}{cc}{J}_{\varphi }\left({\varphi }^{\prime },z^{\prime }\right)=\frac{\partial \psi \left({\varphi }^{\prime },z^{\prime }\right)}{\partial z^{\prime }},J_z\left({\varphi }^{\prime },z^{\prime }\right)=-\frac{1}{a}\frac{\partial \psi \left({\varphi }^{\prime },z^{\prime }\right)}{\partial {\varphi }^{\prime }}\mathrm{on}\rho =a& \left(1.46\right)\end{array}$

The set of values ψ(ϕ′,z)=constant denote the streamlines of the surface coil, which will be the subject of the optimisation, thereby yielding the winding pattern of the surface coil that approximates the streamfunction. The aim of the optimisation is to determine the streamfunction and, therefore, the optimal coil geometry to give a desired magnetic field. Consequently, a Fourier decomposition of the streamfunction is given by

$\begin{array}{cc}\phi \left({\varphi }^{\prime },z^{\prime }\right)=H\left(z^{\prime }+L/2\right)-H\left(z^{\prime }-L/2\right))\left[-{\sum }_{n=1}^N\frac{L}{n\pi }{P}_{n0}\cos \left(\frac{n\pi \left(z^{\prime }+L/2\right)}{L}\right)+{\sum }_{n=1}^N{\sum }_{m=1}^M\frac{L}{n\pi }\left(P_{nm}\cos \left(m{\varphi }^{\prime }\right)+Q_{nm}\sin \left(m{\varphi }^{\prime }\right)\right)\sin \left(\frac{n\pi (z^{\prime }+L/2}{L}\right)\right]& \left(1.47\right)\end{array}$

where H(x) is the Heaviside step function, resulting in an axial current density,

$\begin{array}{cc}{J}_z\left({\varphi }^{\prime },z^{\prime }\right)=\left(H\left(z^{\prime }+L/2\right)-H\left(z^{\prime }-L/2\right)\right)\left\{\sum _{n=1}^N\sum _{m=1}^M\frac{\mathrm{mL}}{n\pi a}\left[P_{nm}\sin \left(m{\varphi }^{\prime }\right)-Q_{nm}\cos \left(m{\varphi }^{\prime }\right)\right]\sin \left(\frac{n\pi (z+L/2}{L}\right)\right\}& \left(1.48\right)\end{array}$

which has zero value at z=±L/2. The azimuthal component of the current density is given by

$\begin{array}{cc}{J}_{\varphi }\left({\varphi }^{\prime },z^{\prime }\right)=\left(H\left(z^{\prime }+L/2\right)-H\left(z^{\prime }-L/2\right)\right)\left[\sum _{n=1}^N{P}_{n0}\sin \left(\frac{n\pi (z^{\prime }+L/2}{L}\right)+\sum _{n=1}^N\sum _{m=1}^M\left(P_{nm}\cos \left(m{\varphi }^{\prime }\right)+Q_{nm}\sin \left(m{\varphi }^{\prime }\right)\right)\cos \left(\frac{n\pi (z^{\prime }+L/2}{L}\right)\right]& \left(1.49\right)\end{array}$

Upon substitution of 1.49 and 1.48 into 1.43, 1.45 and 1.44 respectively, and performing the integration over the surface of the coil, it is found that

$\begin{array}{cc}{B}_{\rho }\left(\rho ,\varphi ,z\right)=\sum _{n=1}^N{P}_{n0}{F}_n\left(\rho ,z\right)+\sum _{n=1}^N\sum _{m=1}^M\left(P_{nm}{G}_{nm}^p\left(\rho ,z\right)+Q_{nm}{G}_{nm}^q\left(\rho ,z\right)\right)& \left(1.5\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }\left(\rho ,\varphi ,z\right)=\sum _{n=1}^N\sum _{m=1}^M\frac{\mathrm{mL}}{n\pi \rho }\left(P_{nm}{H}_{nm}^p\left(\rho ,z\right)+Q_{nm}{H}_{nm}^q\left(\rho ,z\right)\right)& \left(1.51\right)\end{array}$ $\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=\sum _{n=1}^N{P}_{n0}{H}_{n0}\left(\rho ,z\right)+\sum _{n=1}^N\sum _{m=1}^M\left(P_{nm}{S}_{nm}^p\left(\rho ,z\right)+Q_{nm}{S}_{nm}^q\left(\rho ,z\right)\right)& \left(1.52\right)\end{array}$

in which 1.50, 1.51 and 1.52 contain the predefined functions,

$\begin{array}{cc}{F}_n\left(\rho ,z\right)=\frac{{\mu }_{0}a}{\pi }{\int }_0^{\infty }dkk{I}_0^{\prime }\left(k\rho \right)K_0^{\prime }\left(ka\right)\frac{2Ln\pi }{k^2{L}^2-n^2{\pi }^2}\left(\begin{array}{c}\sin \left(\frac{kL}{2}\right)\cos \left(kz\right)\\ \cos \left(\frac{kL}{2}\right)\sin \left(kz\right)\end{array}\right)& \left(1.53\right)\end{array}$ $\begin{array}{cc}{G}_{nm}\left(\rho ,z\right)=\frac{{\mu }_{0}a}{\pi }{\int }_0^{\infty }dkk{I}_m^{\prime }\left(k\rho \right)K_m^{\prime }\left(ka\right)\frac{2k{L}^2}{k^2{L}^2-n^2{\pi }^2}\left(\begin{array}{c}-\sin \left(\frac{kL}{2}\right)\sin \left(kz\right)\\ \cos \left(\frac{kL}{2}\right)\cos \left(kz\right)\end{array}\right)& \left(1.54\right)\end{array}$ $\begin{array}{cc}{G}_{nm}^p\left(\rho ,z\right)=\cos \left(m\varphi \right)G_{nm}\left(\rho ,z\right),G_{nm}^q\left(\rho ,z\right)=\sin \left(m\varphi \right)G_{nm}\left(\rho ,z\right)& \left(1.55\right)\end{array}$ $\begin{array}{cc}{H}_{nm}\left(\rho ,z\right)=\frac{{\mu }_{0}a}{\pi }{\int }_0^{\infty }dkk{I}_m\left(k\rho \right)K_m^{\prime }\left(ka\right)\frac{2Ln\pi }{k^2{L}^2-n^2{\pi }^2}\left(\begin{array}{c}-\sin \left(\frac{kL}{2}\right)\sin \left(kz\right)\\ \cos \left(\frac{kL}{2}\right)\cos \left(kz\right)\end{array}\right)& \left(1.56\right)\end{array}$ $\begin{array}{cc}{H}_{nm}^p\left(\rho ,z\right)=\sin \left(m\varphi \right)H_{nm}\left(\rho ,z\right),H_{nm}^q\left(\rho ,z\right)=-\cos \left(m\varphi \right)H_{nm}\left(\rho ,z\right)& \left(1.57\right)\end{array}$ $\begin{array}{cc}{S}_{nm}\left(\rho ,z\right)=-\frac{{\mu }_{0}a}{\pi }{\int }_0^{\infty }dkk{I}_m\left(k\rho \right)K_m^{\prime }\left(ka\right)\frac{2k{L}^2}{k^2{L}^2-n^2{\pi }^2}\left(\begin{array}{c}\sin \left(\frac{kL}{2}\right)\cos \left(kz\right)\\ \cos \left(\frac{kL}{2}\right)\sin \left(kz\right)\end{array}\right)& \left(1.58\right)\end{array}$ $\begin{array}{cc}{S}_{nm}^p\left(\rho ,z\right)=\cos \left(m\varphi \right)S_{nm}\left(\rho ,z\right),S_{nm}^q\left(\rho ,z\right)=\sin \left(m\varphi \right)S_{nm}\left(\rho ,z\right)& \left(1.59\right)\end{array}$

where the function in each right-most bracket corresponds to the upper line if n is even and to the lower line if n is odd. From this, the magnetic field generated by an arbitrary coil on the surface of an open cylinder may be calculated in Cartesian coordinates,

B _x(ρ,ϕ,z)=cos ϕB_ρ(ρ,ϕ,z)−sin ϕB_ϕ(ρ,ϕ,z)   (1.60)

B _y(ρ,ϕ,z)=sin ϕB_ρ(ρ,ϕ,z)+cos ϕB_ϕ(ρ,ϕ,z)   (1.61)

Using the system of governing equations 1.50 to 1.52, solving for the unknown coefficients (P_n0,P_nm,Q_nm) is a very ill-conditioned problem due to the formulation of the vector potential in integral representation of 1.24 to 1.26. A least squares minimisation may be used, as shown in the example below, although other optimisation techniques may be employed instead. In some embodiments, a penalty term may be implemented in the optimisation to regularise the problem. The regularisation term may take many forms. For example, individual terms in the regularisation term may represent the curvature of a given coil geometry, or the energy stored in the coil. In the example described below, the penalty term relates to the overall power of the coil. However, the choice of the penalty term is somewhat arbitrary as all regularisation parameters act to achieve the same general goal. If the regularisation term is large, the result is a well-conditioned inverse problem with a more simplistic design compromising field fidelity. On the other hand, if the regularisation term is small, the result is a less well-conditioned inverse problem with a more intricate design prioritising higher field fidelity.

The power dissipation in a surface coil of thickness, t, and resistivity, , is given by

$\begin{array}{cc}P={\int }_{-L/2}^{L/2}d{z}^{\prime }{\int }_0^{2\pi }d{\varphi }^{\prime }{❘J_z\left({\varphi }^{\prime },z^{\prime }\right)❘}^2+{❘J_{\varphi }({\varphi }^{\prime },z^{\prime }❘}^2& \left(1.62\right)\end{array}$

which, when integrated over the surface of the coil, gives

$\begin{array}{cc}P=\left[\sum _{n=1}^N{P}_{n0}^2\pi L+\sum _{n=1}^N\sum _{m=1}^M\left(P_{nm}^2+Q_{nm}^2\right)\left(\frac{\pi L}{2}+\frac{m^2{L}^3}{2\pi n^2{a}^2}\right)\right]& \left(1.63\right)\end{array}$

The least squares optimisation may now be represented by the cost function

Φ=αΣ_k[B^desired(r_k)−B(r_k)]²+βP   (1.64)

where α and β are weighting parameters chosen such that the physical parameters may be adjusted to fit physical constraints, with K target field points. The minimisation is achieved by taking the differential of the functional with respect to the Fourier coefficients,

$\begin{array}{cc}\begin{array}{cccc}\frac{\partial \Phi }{\partial P_{i0}}=0,& \frac{\partial \Phi }{\partial P_{ij}}=0,& \frac{\partial \Phi }{\partial Q_{ij}}=0,& i,j>1\end{array}& \left(1.65\right)\end{array}$

allowing the optimal Fourier coefficients to be found for any given physical magnetic field through a matrix inversion. The inversion process ends with the solution of the optimal continuous value of the streamfunction defined on the surface of the specified surface. The final objective, however, is to design a coil which generates the desired magnetic field to a specified accuracy. In some embodiments, therefore, an approximate solution to the current continuum is found by discretising the streamfunction, contouring, therefore finding the streamlines. Every position (x,y,z,φ) may be seen as a hyperplane in 4-dimensional space in which the streamlines exist as the set of hypersurfaces in space (x,y,z,φ_i), which may be projected as a line segment on a 2-dimensional surface, such that

$\begin{array}{cc}\Delta \phi =\frac{\max \phi -\min \phi }{N_c-1},\mathrm{where}{\phi }_i=\min \phi +i\Delta \phi & \left(1.66\right)\end{array}$

with N_c being the number of contours.

Now, consider a hollow cylinder of radius b and length L, of infinite permeability, infinitesimal thickness and planar end caps located at z=±{circumflex over (L)}/2. An arbitrary current flows over an inner open cylinder of radius a<b and length L<{circumflex over (L)} centred about the origin, as shown in FIG. 4. A Green's function needs to be constructed subject to the physical constraints relating to the magnetic field. Physically, the high permeability material changes the magnetisation of the sub-domains in the interior of the passive magnetic shielding material in order to satisfy the condition that

{circumflex over (n)}∧H=0 on S   (1.67)

which is satisfied provided

B(r)_∥=0 on S   (1.68)

This condition requires that

B _ρ|_z=±L/2=0, B_z|_ρ=b=0, B_ϕ|_z=±L/2,ρ=b=0   (1.69)

The induced magnetisation sub-domains of the passive shield may be cast into an equivalent pseudo-current density, which is mathematically equivalent. For clarity, perpendicular surfaces will be dealt with separately. The effect from the cylindrical walls of the passive magnetic shield on the vector potential, A, in the region a<ρ<b, results in 1.37, 1.38 and 1.39 being modified such that

$\begin{array}{cc}{A}_{\rho }\left(\rho ,\varphi ,z\right)=-\frac{i{\mu }_0}{4\pi }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk{e}^{im\varphi }{e}^{ikz}\left\{a\left[I_{m-1}\left(❘k❘a\right)K_{m-1}\left(❘k❘\rho \right)-I_{m+1}\left(❘k❘a\right)K_{m+1}\left(❘k❘\rho \right)\right]J_{\varphi }^m\left(k\right)+b[I_{m-1}\left(❘k❘\rho \right)K_{m-1}\left(❘k❘b\right)-I_{m+1}\left(❘k❘\rho \right)K_{m+1}\left(❘k❘b\right)\left(k\right)\right\}& \left(1.7\right)\end{array}$ $\begin{array}{cc}{A}_{\varphi }\left(\rho ,\varphi ,z\right)=\frac{{\mu }_0}{4\pi }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk{e}^{im\varphi }{e}^{ikz}\left\{a\left[I_{m-1}\left(❘k❘a\right)K_{m-1}\left(❘k❘\rho \right)+I_{m+1}\left(❘k❘a\right)K_{m+1}\left(❘k❘\rho \right)\right]J_{\varphi }^m\left(k\right)+b[I_{m-1}\left(❘k❘\rho \right)K_{m-1}\left(❘k❘b\right)+I_{m+1}\left(❘k❘\rho \right)K_{m+1}\left(❘k❘b\right)\left(k\right)\right\}& \left(1.71\right)\end{array}$ $\begin{array}{cc}{A}_z\left(\rho ,\varphi ,z\right)=\frac{{\mu }_0}{2\pi }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk{e}^{im\varphi }{e}^{ikz}\left[a{I}_m\left(❘k❘a\right)K_m\left(❘k❘\rho \right)J_z^m\left(k\right)+b{I}_m\left(❘k❘\rho \right)K_m\left(❘k❘b\right)\left(k\right)\right]& \left(1.72\right)\end{array}$

Using 1.40 and 1.69, the boundary conditions at ρ=b may be written as

$\begin{array}{cc}{\frac{1}{\rho }\left(\frac{\partial }{\partial \rho }\left(\rho A_{\varphi }\right)-\frac{\partial A_{\rho }}{\partial \varphi }\right)❘}_{\rho =b}=0& \left(1.73\right)\end{array}$ $\begin{array}{cc}{\left(\frac{\partial A_{\rho }}{\partial z}-\frac{\partial A_z}{\partial \rho }\right)❘}_{\rho =b}=0& \left(1.74\right)\end{array}$

Using 1.70, 1.71 and 1.72, the pseudo-current density on the cylindrical wall of the passive shield is found to be

$\begin{array}{cc}\left(k\right)=-\frac{a{I}_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{b{I}_m\left(❘k❘b\right)K_m^{\prime }\left(❘k❘b\right)}{J}_z^m\left(k\right)& \left(1.75\right)\end{array}$

It follows from 1.40, 1.71, 1.71 and 1.72 that the magnetic field in the presence of a passive magnetic shielding material (e.g., a high permeability magnetic material) in the region ρ<a may be written as

$\begin{array}{cc}{B}_{\rho }\left(\rho ,\varphi ,z\right)=\frac{i{\mu }_{0}a}{2\pi }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dkk{e}^{im\varphi }{e}^{ikz}{I}_m^{\prime }\left(❘k❘\rho \right)\left[K_m^{\prime }\left(❘k❘a\right)-\frac{I_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{I_m\left(❘k❘b\right)}\right]J_{\varphi }^m\left(k\right)& \left(1.76\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }\left(\rho ,\varphi ,z\right)=\frac{{\mu }_0{a}^2}{2\mathrm{\pi \rho }}\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk❘k❘e^{im\varphi }{e}^{ikz}{I}_m\left(❘k❘\rho \right)\left[K_m^{\prime }\left(❘k❘a\right)-\frac{I_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{I_m\left(❘k❘b\right)}\right]J_z^m\left(k\right)& \left(1.77\right)\end{array}$ $\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=-\frac{{\mu }_{0}a}{2\pi }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk❘k❘e^{im\varphi }{e}^{ikz}{I}_m\left(❘k❘\rho \right)\left[K_m^{\prime }\left(❘k❘a\right)-\frac{I_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{I_m\left(❘k❘b\right)}\right]J_{\varphi }^m\left(k\right)& \left(1.78\right)\end{array}$

In the specific case of the radius of the coil being equal to the radius of the passive magnetic shield, (a=b), these expressions may be further simplified using the relation 1.19, giving

$\begin{array}{cc}{B}_{\rho }\left(\rho ,\varphi ,z\right)=\frac{i{\mu }_0}{2\pi }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk\frac{k}{❘k❘}{e}^{im\varphi }{e}^{ikz}\frac{I_m\left(❘k❘\rho \right)}{I_m\left(❘k❘a\right)}{J}_{\varphi }^m\left(k\right)& \left(1.79\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }\left(\rho ,\varphi ,z\right)=\frac{{\mu }_{0a}}{2\mathrm{\pi \rho }}\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk{e}^{im\varphi }{e}^{ikz}\frac{I_m\left(❘k❘\rho \right)}{I_m\left(❘k❘a\right)}{J}_z^m\left(k\right)& \left(1.8\right)\end{array}$ $\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=-\frac{{\mu }_0}{2\pi }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk{e}^{im\varphi }{e}^{ikz}\frac{I_m\left(❘k❘\rho \right)}{I_m\left(❘k❘a\right)}{J}_{\varphi }^m\left(k\right)& \left(1.81\right)\end{array}$

however, for generality this will be not be used (but it is noted that it may be used). Now, the effect of the planar end caps needs to be incorporated. Using 1.40 and 1.69, the boundary conditions on the magnetic field at z=±{circumflex over (L)}/2 may be written as

$\begin{array}{cc}{\left(\frac{1}{\rho }\frac{\partial A_z}{\partial \varphi }-\frac{\partial A_{\varphi }}{\partial z}\right)❘}_{z\pm \hat{L}/2}=0& \left(1.82\right)\end{array}$ $\begin{array}{cc}{\left(\frac{\partial A_{\rho }}{\partial z}-\frac{\partial A_z}{\partial \rho }\right)❘}_{z\pm \hat{L}/2}=0& \left(1.83\right)\end{array}$

In the below example, these boundary conditions are implemented through the method of mirror images. The resulting solution comprises an infinite number of reflected image currents giving a magnetic field with components of the form

$\begin{array}{cc}{B}_{\rho }\left(\rho ,\varphi ,z\right)=\frac{i{\mu }_{0}a}{2\pi }\sum _{p=-\infty }^{\infty }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dkk{e}^{im\varphi }{e}^{ikz}{I}_m^{\prime }\left(❘k❘\rho \right)\left[K_m^{\prime }\left(❘k❘a\right)-\frac{I_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{I_m\left(❘k❘b\right)}\right]J_{\varphi }^{mp}\left(k\right)& \left(1.84\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }\left(\rho ,\varphi ,z\right)=\frac{{\mu }_0{a}^2}{2\mathrm{\pi \rho }}\sum _{p=-\infty }^{\infty }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk❘k❘e^{im\varphi }{e}^{ikz}{I}_m\left(❘k❘\rho \right)\left[K_m^{\prime }\left(❘k❘a\right)-\frac{I_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{I_m\left(❘k❘b\right)}\right]J_z^{mp}\left(k\right)& \left(1.85\right)\end{array}$ $\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=-\frac{{\mu }_{0}a}{2\pi }\sum _{p=-\infty }^{\infty }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk❘k❘e^{im\varphi }{e}^{ikz}{I}_m\left(❘k❘\rho \right)\left[K_m^{\prime }\left(❘k❘a\right)-\frac{I_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{I_m\left(❘k❘b\right)}\right]J_{\varphi }^{mp}\left(k\right)& \left(1.86\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{J}_{\varphi }^{mp}\left(k\right)=\frac{1}{2\pi }{\int }_0^{2\pi }d{\varphi }^{\prime }{e}^{-im{\varphi }^{\prime }}{\int }_{-\infty }^{\infty }d{z}^{\prime }{e}^{-ik{z}^{\prime }}{J}_{\varphi }^p\left({\varphi }^{\prime },z^{\prime }\right)& \left(1.87\right)\end{array}$ $\begin{array}{cc}{J}_z^{mp}\left(k\right)=\frac{1}{2\pi }{\int }_0^{2\pi }d{\varphi }^{\prime }{e}^{-im{\varphi }^{\prime }}{\int }_{-\infty }^{\infty }d{z}^{\prime }{e}^{-ik{z}^{\prime }}{J}_z^p\left({\varphi }^{\prime },z^{\prime }\right)& \left(1.88\right)\end{array}$

As is usual for the method of mirror images, the image currents must therefore be adjusted such that any odd reflections, e.g., p=1, 3, 5, . . . etc., result in an axial current that opposes the original current source. This leads to an adjusted streamfunction of the form

$$\begin{gathered} \begin{array}{cc}{\phi }^p(\text{⁠}\left({\varphi }^{\prime },z^{\prime }\right)=\left(T_e\left(z^{\prime };L_2+p\hat{L},L_1+p\hat{L}\right)+T_o\left(z^{\prime };-L_1+p\hat{L},-L_2+p\hat{L}\right)\right)\times \text{ \\ {-(-1)p⁢∑n=1NLn⁢π⁢Pn⁢0⁢cos⁡(n⁢π⁡((-1)p⁢(z′-p⁢Lˆ)-L2)L1-L2)⁠+(-1)p⁢⁠⁠∑n=1N⁢⁠∑m=1MLn⁢π⁢(Pn⁢m⁢cos⁡(m⁢ϕ′+Qn⁢m⁢sin⁡(m⁢ϕ′))⁢sin⁡(n⁢π⁡((-1)p⁢(z′-p⁢Lˆ)-L2)L1-L2)}}& \left(1.89\right)\end{array} \end{gathered}$$ $\begin{array}{cc}(T_e\left(z^{\prime };L_2+p\hat{L},L_1+p\hat{L}\right)=\left(H\left(z^{\prime }-L_2-p\hat{L}\right)-H\left(z^{\prime }-L_1-p\hat{L}\right)\right)\left(\frac{1-{\left(-1\right)}^{p+1}}{2}\right)& \left(1.9\right)\end{array}$ $\begin{array}{cc}{T}_o\left(z^{\prime };-L_1+p\hat{L},-L_2+p\hat{L}\right))=\left(H\left(z^{\prime }-L_1-p\hat{L}\right)-H\left(z^{\prime }-L_2-p\hat{L}\right)\right)\left(\frac{1-{\left(-1\right)}^p}{2}\right)& \left(1.91\right)\end{array}$

where H(x) is the Heaviside function, resulting in an axial current density

$$\begin{gathered} \begin{array}{cc}{J}_z^p(\text{⁠}\left({\varphi }^{\prime },z^{\prime }\right)=\left(T_e\left(z^{\prime };L_2+p\hat{L},L_1+p\hat{L}\right)+T_o\left(z^{\prime };-L_1+p\hat{L},-L_2+p\hat{L}\right)\right)\times \text{ \\ ∑n=1N∑m=1M(-1)p⁢m⁢Ln⁢π⁢a[Pn⁢m⁢sin⁡(m⁢ϕ′)-Qn⁢m⁢cos⁡(m⁢ϕ′)]⁢sin⁡(n⁢π((-1)p⁢(z′-p⁢Lˆ)-L2)L1-L2)}& \left(1.92\right)\end{array} \end{gathered}$$

which has zero value at z′=±L_1,2+p{circumflex over (L)} for even and odd reflections respectively. The azimuthal component of the current density is given by

$$\begin{gathered} \begin{array}{cc}{J}_{\varphi }^p\left({\varphi }^{\prime },z^{\prime }\right)=\left(T_e\left(z^{\prime };L_2+p\hat{L},L_1+p\hat{L}\right)+T_o\left(z^{\prime };-L_1+p\hat{L},-L_2+p\hat{L}\right)\right)\times \text{ \\ {∑n=1NPn⁢0⁢sin⁡(n⁢π⁡((-1)p⁢(z′-p⁢Lˆ)-L2)L1-L2)+∑n=1N∑m=1M(Pn⁢m⁢cos⁡(m⁢ϕ′+ \\ Qn⁢m⁢sin⁡(m⁢ϕ′))⁢cos⁡(n⁢π⁡((-1)p⁢(z′-p⁢Lˆ)-L2)L1-L2)}}& \left(1.93\right)\end{array} \end{gathered}$$

Substituting 1.92 and 1.93 into 1.84, 1.85 and 1.86, it is found that the magnetic field is given by

$\begin{array}{cc}{B}_{\rho }\left(\rho ,\varphi ,z\right)=\sum _{n=1}^N{P}_{n0}\widetilde{F_n}\left(\rho ,z\right)+\sum _{n=1}^N\sum _{m=1}^M\left(P_{nm}{\tilde{G}}_{nm}^p\left(\rho ,z\right)+Q_{nm}{\tilde{G}}_{nm}^q\left(\rho ,z\right)\right)& \left(1.94\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }\left(\rho ,\varphi ,z\right)=\sum _{n=1}^N\sum _{m=1}^M\left(P_{nm}{\tilde{H}}_{nm}^p\left(\rho ,z\right)+Q_{nm}{\tilde{H}}_{nm}^q\left(\rho ,z\right)\right)& \left(1.95\right)\end{array}$ $\begin{array}{cc}{B}_Z\left(\rho ,\varphi ,z\right)=\sum _{n=1}^N{P}_{n0}{\tilde{D}}_n\left(\rho ,z\right)+\sum _{n=1}^N\sum _{m=1}^M\left(P_{nm}{\tilde{S}}_{nm}^p\left(\rho ,z\right)+Q_{nm}{\tilde{S}}_{nm}^q\left(\rho ,z\right)\right)& \left(1.96\right)\end{array}$ $\mathrm{where},$ $\begin{array}{cc}{F}_n\left(\rho ,z\right)=\sum _{p=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dkI}}_0^{\prime }\left(❘k❘\rho \right)\left(K_0^{\prime }\left(❘k❘a\right)-\frac{I_0^{\prime }\left(❘k❘a\right)K_0\left(❘k❘b\right)}{I_0\left(❘k❘b\right)}\right)C_{np}^1\left(k,z\right)& \left(1.97\right)\end{array}$ $\begin{array}{cc}{G}_{nm}\left(\rho ,z\right)=\sum _{p=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dkI}}_m^{\prime }\left(❘k❘\rho \right)\left(K_m^{\prime }\left(❘k❘a\right)-\frac{I_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{I_m\left(❘k❘b\right)}\right)C_{np}^2\left(k,z\right)& \left(1.98\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}{G}_{nm}^p\left(\rho ,z\right)=\cos \left(m\varphi \right)G_{nm}\left(\rho ,z\right),& G_{nm}^q\left(\rho ,z\right)=\sin \left(m\varphi \right)G_{nm}\left(\rho ,z\right)\end{array}& \left(1.99\right)\end{array}$ $\begin{array}{cc}{H}_{nm}\left(\rho ,z\right)=\sum _{p=-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathrm{dkm}{I}_m\left(k\rho \right)\left(K_m^{\prime }\left(❘k❘a\right)-\frac{I_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{I_m\left(❘k❘b\right)}\right)C_{np}^3\left(k,z\right)& \left(1.1\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}{H}_{nm}^p\left(\rho ,z\right)=\sin \left(m\varphi \right)H_{nm}\left(\rho ,z\right),& H_{nm}^q\left(\rho ,z\right)=-\cos \left(m\varphi \right)H_{nm}\left(\rho ,z\right)\end{array}& \left(1.101\right)\end{array}$ $\begin{array}{cc}{D}_n\left(\rho ,z\right)=\sum _{p=-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathrm{dk}{I}_0\left(k\rho \right)\left(K_0^{\prime }\left(❘k❘a\right)-\frac{I_0^{\prime }\left(❘k❘a\right)K_0\left(❘k❘b\right)}{I_0\left(❘k❘b\right)}\right)C_{np}^4\left(k,z\right)& \left(1.102\right)\end{array}$ $\begin{array}{cc}{S}_{nm}\left(\rho ,z\right)=\sum _{p=-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\mathrm{dkI}}_m\left(k\rho \right)\left(K_m^{\prime }\left(❘k❘a\right)-\frac{I_m^{\prime }\left(❘k❘a\right)K_m\left(❘k❘b\right)}{I_m\left(❘k❘b\right)}\right)C_{np}^5\left(k,z\right)& \left(1.103\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}{S}_{nm}^p\left(\rho ,z\right)=\cos \left(m\varphi \right)S_{nm}\left(\rho ,z\right),& S_{nm}^q\left(\rho ,z\right)=\sin \left(m\varphi \right)S_{nm}\left(\rho ,z\right)\end{array}& \left(1.104\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{C}_{\mathrm{np}}^1\left(k,z\right)=\frac{i{\mu }_0\mathrm{an}\left(L_1-L_2\right)}{2}{\mathrm{ke}}^{\mathrm{ik}\left(z-p\hat{L}\right)}\left(\frac{e^{-{\left(-1\right)}^P{\mathrm{ikL}}_2}+{\left(-1\right)}^{n+1}{e}^{-{\left(-1\right)}^P{\mathrm{ikL}}_1}}{n^2{\pi }^2-{\left(L_1-L_2\right)}^2{k}^2}\right)& \left(1.105\right)\end{array}$ $\begin{array}{cc}{C}_{\mathrm{np}}^2\left(k,z\right)=\frac{i{\left(-1\right)}^p\left(L_1-L_2\right)k}{n\pi }{C}_{\mathrm{np}}^1\left(k,z\right)& \left(1.106\right)\end{array}$ $\begin{array}{cc}{C}_{\mathrm{np}}^3\left(k,z\right)=\frac{i{\left(-1\right)}^p\left(L_1-L_2\right)❘k❘}{n\pi k\rho }{C}_{\mathrm{np}}^1\left(k,z\right)& \left(1.107\right)\end{array}$ $\begin{array}{cc}{C}_{\mathrm{np}}^4\left(k,z\right)=\frac{i❘k❘}{k}{C}_{\mathrm{np}}^1\left(k,z\right)& \left(1.108\right)\end{array}$ $\begin{array}{cc}{C}_{\mathrm{np}}^5\left(k,z\right)=-\frac{{\left(-1\right)}^p\left(L_1-L_2\right)❘k❘}{n\pi }{C}_{\mathrm{np}}^1\left(k,z\right)& \left(1.109\right)\end{array}$

As outlined previously, the Fourier coefficients may be solved using an optimisation process (for example a least squares minimisation, optionally incorporating a regularisation term).

Whilst the example method outlined above assumes infinite permeability of the cylinder, magnetic shields comprising a passive magnetic shielding material having a relative permeability of substantially 1000 or greater, and a winding configuration optimised in accordance with the invention, will generate a magnetic field that is accurate to within substantially 99% of the magnetic field produced by the system optimised assuming infinite permeability. Magnetic shields having a winding configuration optimised in accordance with the invention and comprising a passive magnetic shielding material having a relative permeability of less than substantially 1000 may generate magnetic fields with a greater difference from the target field, but may still be useful in providing magnetic shielding.

A similar approach to that described above may be used in order to optimise winding configurations accounting for the interactions with a passive magnetic shielding material to produce magnetic fields other than constant magnetic fields, for example linear ‘gradient’ fields or quadratic ‘curvature’ fields. This is discussed in further detail below.

It will also be appreciated that depending on the region within a cylindrical passive magnetic shield over which the optimisation of the winding is to be considered (for example, if the cylindrical passive magnetic shield has a high aspect ratio, with the optimised region located near the centre of the cylindrical passive magnetic shield with respect to its axial direction), the planar end caps of the cylindrical passive magnetic shield may be disregarded for the purposes of optimising the winding configuration. In such embodiments, the boundary conditions relating to the planar end caps need not be implemented, and the optimisation process performed accordingly. Therefore, the optimisation may alternatively be performed modelling the passive magnetic shield as an open cylindrical structure, rather than a closed cylindrical structure.

It will further be appreciated that a similar approach to that described above may be employed to optimise a winding configuration for use with a passive magnetic shield having a geometry other than a cylindrical geometry. By employing the same principle of constructing a physically relevant function (for example, a Green's function), or discrete approximation thereof, for the interior geometry of the passive magnetic shield, and implementing one or more appropriate boundary conditions (depending on the geometry of the passive magnetic shield) to ensure the value of the magnetic field attains the correct physical limit as the winding approaches the interior surface of the passive magnetic shield.

In alternative embodiments, a finite element analysis approach is employed instead of the Fourier analysis illustrated in the above example. In some embodiments, the governing equations and boundary conditions as shown in 1.84 to 1.86 are used in a numerical calculation with a finite element current distribution.

In a further alternative embodiment, rather than determining optimal surface currents as described in the above example, one or more discrete winding elements are specified a priori and one or more parameters optimised to produce a specified magnetic field. This approach can sometimes be simpler than the surface currents approach discussed above, for example depending on engineering constraints of the magnetic shield in question. This alternative approach makes use of some of the method described in detail above (up to 1.88). The boundary conditions are implemented for a finite length perfect magnetic conducting (i.e., infinite magnetic permeability) cylinder of length 2L and outer radius b enclosing an interior cylindrical region of radius a.

For simplicity, the case a=b is discussed below. The resulting field in the z-direction is given by

$\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=\frac{{\mu }_0}{2\pi }\sum _{p=-\infty }^{\infty }\sum _{m=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk{e}^{im\varphi }{e}^{ikz}\frac{I_m\left(❘k❘\rho \right)}{I_m\left(❘k❘a\right)}{J}_{\varphi }^{mp}\left(k\right)& \left(1.11\right)\end{array}$

The simplest discrete winding element is now considered—a circle. For a single circular loop with current I_i, of radius a and z location d inside the cylinder, the resulting field in the z-direction may be written as

$\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=\frac{{\mu }_0{I}_i}{2\pi }\sum _{p=-\infty }^{\infty }{\int }_{-\infty }^{\infty }dk{e}^{ik\left(z-{\left(-1\right)}^{P_{d_i-2L_p}}\right)}\frac{I_0\left(❘k❘\rho \right)}{I_0\left(❘k❘a\right)}& \left(1.111\right)\end{array}$

A simple gradient field may be generated in the z-direction using a single pair of circular loops with currents flowing in opposite direction (an anti-Helmholtz arrangement). Using 1.111, the magnetic field produced by such a pair of coils separated by a distance along the z-direction of 2d from the origin may be written as

$\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=\frac{-2i{\mu }_{0}I}{\pi }{\int }_{-\infty }^{\infty }dk{e}^{ikz}\frac{I_0\left(❘k❘\rho \right)}{I_0\left(❘k❘a\right)}\sin \left(kd\right)\left(1-2\sum _{p=1}^{\infty }{\left(-1\right)}^{p+1}\cos \left(2Lpk\right)\right)& \left(1.112\right)\end{array}$

Performing the integral over k it is found that

$\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=2{\mu }_{0}I\sum _{p=-\infty }^{\infty }\sin \left(\frac{\pi pz}{2L}\right)\sin \left(\frac{\pi \mathrm{pd}}{2L}\right)\frac{I_0\left(❘\frac{\pi p}{2L}❘\rho \right)}{I_0\left(❘\frac{\pi p}{2L}❘a\right)}\left(1-{\left(-1\right)}^p\right)& \left(1.113\right)\end{array}$

Taylor expanding the spatially varying terms using

$\begin{array}{cc}{I}_0\left(\frac{k\rho }{a}\right)=1+\frac{k^2{\rho }^2}{4a^2}+O\left({\rho }^4\right)+\dots & \left(1.114\right)\end{array}$ $\begin{array}{cc}\sin \left(\frac{\mathrm{kz}}{a}\right)=\frac{kz}{a}-\frac{k^3{z}^3}{6a^3}+O\left(z^4\right)-\dots & \left(1.115\right)\end{array}$

it is seen that

$\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=4{\mu }_0{I}_i\sum _{p=1}^{\infty }p\frac{\sin \left(\frac{\pi \mathrm{dp}}{2L}\right)}{I_0\left(\frac{\pi \mathrm{ap}}{2L}\right)}\left(1-{\left(-1\right)}^p\right)\left[\frac{\pi z}{2L}+p^2\left(\frac{{\pi }^{3}z{\rho }^2}{32L^3}-\frac{{\pi }^3{z}^3}{48L^3}\right)+p^4\dots \right]& \left(1.116\right)\end{array}$

1.116 may be simplified to

$\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=\frac{4{\mu }_0\pi I}{2L}\left(z{B}_1+\left(\frac{{\pi }^{2}z{\rho }^2}{16L^2}-\frac{{\pi }^2{z}^3}{24L^2}\right)B_3+\dots \right)& \left(1.117\right)\end{array}$

where B₁, B₃, . . . B_2n−1 are the magnitudes of the linear, cubic, etc. fields generated by the asymmetric coil set up. In order to remove cubic variations in the field it is required that B₃=0, resulting, in order to cancel cubic variations in the magnetic field, in

$\begin{array}{cc}{B}_3=\sum _{p=1}^{\infty }{p}^3\frac{\sin \left(\frac{\pi \mathrm{dp}}{2L}\right)}{I_0\left(\frac{\pi \mathrm{ap}}{2L}\right)}\left(1-{\left(-1\right)}^p\right)=0& \left(1.118\right)\end{array}$

which may be generalised for N pairs of loops resulting in the set of simultaneous equations

$\begin{array}{cc}\begin{array}{cc}\sum _{i=1}^N{I}_i\sum _{p=1}^{\infty }{p}^{2n+1}\frac{\sin \left(\frac{\pi \mathrm{dp}}{2L}\right)}{I_0\left(\frac{\pi \mathrm{ap}}{2L}\right)}\left(1-{\left(-1\right)}^P\right)=0,& n\in \mathbb{Z}:n\in \left[1,N\right]\end{array}& \left(1.119\right)\end{array}$

The ratios of a and d may now be determined in order to eliminate the leading order errors in the desired magnetic field. In the specific embodiment described above, the free parameters optimised are the spacing d and the relationship between the spacing d and the radius a of the circular loop winding elements. Other free parameters may be optimised, for example, an orientation, shape, size or placement of a discrete winding element, to generate a desired magnetic field. Additionally or alternatively, relationships between two or more free parameters of the discrete winding elements may be optimised to generate a desired magnetic field.

For an effective approach to magnetic field design, the formation of any physical field must be understood from a mathematical basis. Starting from the fundamental equations for static magnetic fields, as the magnetic field, H, and the electric field, E, may be expressed in terms of the vector potential, A, as well as the scalar potential, Φ, any transformation of these variables leads to an arbitrary choice of gauge in the way the potential is defined. As the problems considered here are in the static regime, using the Coulomb gauge is the most logical choice for simplicity, and so, unless stated, this gauge is used throughout. The Maxwell equations that describe static magnetic fields, through a volume within which no current passes, are given by Gauss' and Ampere's law respectively,

∇·H=∇∧H=0   (1.120)

with the magnetic induction, B, being related to the scalar and vector potential by

B=−∇Φ=∇∧A   (1.121)

The magnetic induction may be related to the magnetic field through the use of the magnetic permeability, μ. For free space, this relation takes the simplified form of

B=μ ₀(H+M)   (1.122)

Where M is the magnetisation, relating to the bound current density of a magnetised surface

J_b=∇∧M   (1.123)

Upon substitution of 1.111 into 1.110 the following is obtained

∇·∇Φ=∇∧∇∧A=0   (1.124)

Using vector identities with choice of gauge, ∇·A=0, it is seen that the scalar and vector potential satisfy the partial differential equation (PDE),

∇²Φ=∇²A=0   (1.125)

It then follows that

∇²H=∇²B=0   (1.126)

Laplace's may be solved by a transformation into spherical polar coordinates such that the Laplacian of a function is given by

$\begin{array}{cc}{\nabla }^{2}f=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial f}{\partial r}\right)+\frac{1}{r^2{\sin }^2\theta }\frac{{\partial }^{2}f}{\partial {\varphi }^2}+\frac{1}{r^2\sin \theta }\left(\sin \theta \frac{\partial f}{\partial \theta }\right)& \left(1.127\right)\end{array}$

Separation of variables shows that the real solutions take the form

$\begin{array}{cc}f=\sum _{n=1}^{\infty }\sum _{m=-n}^n{C}_{nm}\left(\begin{array}{c}{r}_{<}^n\\ r_{>}^{-n-1}\end{array}\right)P_{nm}\left(\cos \theta \right)\left(\begin{array}{c}\cos \left(m\varphi \right)\\ \sin \left(m\varphi \right)\end{array}\right),\begin{array}{c}m\ge 0\\ m<0\end{array}& \left(1.128\right)\end{array}$

where P_nm(x) are Ferrer's associated Legendre polynomials of the nth order and mth degree, with r_< and r_> relating to the interior or exterior region of a defined spherical topology enclosing a current source. An orthogonal basis may be generated for both the scalar and vector magnetic potential as well as the individual magnetic field components and intensities. From this, the physical orthogonal set of magnetic fields may be constructed from the relationship between the scalar potential and the magnetic field. For a magnetic field generated on the interior of a surface current (r_< case), the magnetic field may be expanded as a summation of spherical harmonic fields. The harmonics produced, up to and including n=3, may be given by

Constant “Bias” Fields

$\begin{array}{cc}n=1\to \{\begin{array}{c}B\left(x,y,z\right)=\hat{x},m=1\\ B\left(x,y,z\right)=\hat{z},m=0\\ B\left(x,y,z\right)=\hat{y},m=-1\end{array}& \left(1.129\right)\end{array}$

Linear “Gradient” Fields

$\begin{array}{cc}n=2\to \{\begin{array}{c}B\left(x,y,z\right)=x\hat{x}-y\hat{y},m=2\\ B\left(x,y,z\right)=z\hat{x}+x\hat{z},m=1\\ B\left(x,y,z\right)=-x\hat{x}-y\hat{y}+2z\hat{z},m=0\\ B\left(x,y,z\right)=z\hat{y}+y\hat{z},m=-1\\ B\left(x,y,z\right)=y\hat{x}+x\widehat{y,}m=-2\end{array}& \left(1.13\right)\end{array}$

Quadratic “Curvature” Fields

$\begin{array}{cc}n=2\to \{\begin{array}{c}B\left(x,y,z\right)=\left(x^2-y^2\right)\hat{x}-2xy\widehat{y,}m=3\\ B\left(x,y,z\right)=2xz\hat{x}-2yz\hat{y}+\left(x^2-y^2\right)\hat{z},m=2\\ B\left(x,y,z\right)=\left(4z^2-3x^2-y^2\right)\hat{x}-2xy\hat{y}+8xz\hat{z},m=1\\ B\left(x,y,z\right)=-2xz\hat{x}-2yz\hat{y}+\left(2z^2-x^2-y^2\right)\widehat{z,}m=0\\ B\left(x,y,z\right)=-2xy\hat{x}+\left(4z^2-x^2-3y^2\right)\hat{y}+8yz\hat{z},m=-1\\ B\left(x,y,z\right)=yz\hat{x}+xz\hat{y}+xy\hat{z},m=-2\\ B\left(x,y,z\right)=2xy\hat{x}-\left(x^2-y^2\right)\hat{y}m=-3\end{array}& \left(1.131\right)\end{array}$

Through the use of this orthogonal set, a cancellation system may be constructed to null magnetic fields efficiently. This may be extended to higher order fields as needed, although the relative magnitude of an external field of the same order and degree decreases rapidly as n increases. Due to this, a limited set of spherical harmonics may be used to facilitate the bulk of the cancellation. Winding configurations may be optimised, for example, to generate one of the spherical harmonic fields outlined above.

Existing magnetic shields typically comprise passive magnetic shielding material, and a winding configuration designed, or optimised, not in the presence of the passive magnetic shielding.

Constant field winding configurations designed to fit the internal cavity of a substantially cylindrical magnetic shield having internal cavity length of 16.4 cm and an internal cavity diameter of 14.4 cm are now presented. A comparison is performed using two winding configurations—one configuration optimised omitting the interaction with the passive magnetic shielding material (e.g., Mu-metal®), and one configuration optimised accounting for the interaction with the passive magnetic shielding material.

Firstly, winding configurations optimised omitting the interaction with the passive magnetic shielding material are considered. FIG. 5A shows an optimised (omitting the interaction with the passive magnetic shielding material) constant field winding 501 for the z-direction, n=1,m=0 (1.129). FIG. 5B shows a geometry of the winding 501 plotted versus angular position, ϕ, around the coil and z for clarity. FIG. 5C shows an optimised (omitting the interaction with the passive magnetic shielding material) constant field winding 502 for the y/x direction, n=1,m=1/−1 (1.129). FIG. 5D shows a geometry of the winding 502 plotted versus angular position, ϕ, around the coil and z for clarity. The windings 501, 502 are optimised over a region of length 9 cm and diameter 6 cm (i.e., a region having a volume of over 1000 cm³), represented by the cylindrical region displayed within the windings 501, 502 as shown in FIGS. 5A and 5C respectively.

The resulting relative field homogeneities of the windings 501, 502 are shown in FIGS. 5E and 5F respectively. FIG. 5E shows a plot of relative field variation (in %) as a function of distance from the central point (indicated by relative position of zero in the plot of FIG. 5E) of the space within the winding 501, along the z-axis, for the winding 501. FIG. 5F shows a plot of relative field variation (in %) as a function of distance from the central point (indicated by relative position of zero in the plot of FIG. 5F) of the space within the winding 502, along the z-axis, for the winding 502. The field homogeneities are given by 100×(1−B_i), where B_i is the unitary magnetic field in a given direction i=x,y,z (e.g., B_x=1, B_y=0, B_z=0).

The approximated relative field homogeneities of the windings 501, 502 are found when those windings 501, 502 are enclosed in a capped Mu-metal® cylinder of the same length and diameter of the windings 501, 502, and calculated using the method of mirror images (1.84 to 1.86). The approximated relative field homogeneities for the winding 501 and the winding 502 are shown in FIGS. 5G and 5H respectively. It can be seen that the effect of the passive magnetic shielding material on the relative magnetic field homogeneity for both z and y/x constant field coils decreases the relative homogeneity substantially ten-fold. In other words, the windings 501, 502 perform measurably worse when in the presence of passive magnetic shielding material. This can be seen by comparing FIG. 5E with FIG. 5G, and by comparing FIG. 5F with FIG. 5H.

Now, winding configurations optimised including the interaction with passive magnetic shielding material (e.g., Mu-metal®) are presented. FIG. 6A shows an optimised (including the interaction with passive magnetic shielding material) constant field winding 601 for the z-direction, n=1,m=0 (1.129). FIG. 6B shows a geometry of the winding 601 plotted versus angular position, ϕ, around the coil and z for clarity. FIG. 6C shows an optimised (including the interaction with passive magnetic shielding material) constant field winding 602 for the y/x direction. FIG. 6D shows a geometry of the winding 602 plotted versus angular position, ϕ, around the coil and z for clarity. In the embodiment shown, the windings 601, 602 are optimised using the method outlined above (implementing relevant boundary conditions to account for the passive magnetic shielding material, and determining optimal surface currents). The windings 601, 602 are optimised over a region of length 9 cm and diameter 6 cm (i.e., a region having a volume of over 1000 cm³), represented by the cylindrical region displayed within the windings 601, 602 as shown in FIGS. 6A and 6C respectively.

The relative field homogeneities of the windings 601, 602, when enclosed in a capped Mu-metal® cylinder of the same length and diameter of the coil, are shown in FIGS. 6E and 6F respectively. FIG. 6E shows a plot of relative field variation (in %) as a function of distance from the central point (indicated by relative position of zero in the plot of FIG. 6E) of the space within the winding 601, along the z-axis, for the winding 601. FIG. 6F shows a plot of relative field variation (in %) as a function of distance from the central point (indicated by relative position of zero in the plot of FIG. 6F) of the space within the winding 602, along the z-axis, for the winding 602. The relative field homogeneities are calculated using the method of mirror images (1.84 to 1.86). It will be noted that because the optimisation of windings 601, 602 accounts for the interaction of the windings 601, 602 with the passive magnetic shielding material, the performance (for example, the field homogeneity) of the windings 601, 602 will be noticeably or measurably worse when not in the presence of passive magnetic shielding material (similar but opposite to the behaviour of the windings 501, 502). The field produced by a winding optimised accounting for the interactions with passive magnetic shielding material (e.g., windings 601, 602) will be more similar to a desired field (e.g., the desired field the optimisation is aimed towards producing) when the winding is in the presence of passive magnetic shielding material than when the winding is not in the presence of passive magnetic shielding material. This is true for uniform fields, linear ‘gradient’ fields, quadratic ‘curvature’ fields, and higher order fields. In the particular case of a uniform field, the field produced by such a winding (e.g., windings 601, 602) will be more uniform in the presence of passive magnetic shielding material than not in the presence of passive magnetic shielding material.

FIG. 6G shows the direction of current flow in the various regions of the optimised winding configuration shown in FIG. 6D. The different loops indicated A and B in FIG. 6G indicate counter-flowing current. Current flowing in regions labelled A flows in the opposite direction to current flowing in regions labelled B. FIG. 6H shows an example of how the different loops of the optimised winding configuration can be connected together to enable counter-flowing current to flow throughout the winding to produce a desired magnetic field.

A comparison of the winding configurations optimised omitting (windings 501, 502) or including (windings 601, 602) the interaction with a passive magnetic shielding material reveals drastically different winding geometries/configurations. Moreover, the resulting relative field homogeneity obtained using a winding for which the interactions with the passive magnetic shielding material are taken into account ab initio (windings 601, 602) show a substantially ten-fold reduction of unwanted field inhomogeneity (see FIGS. 6E, 6F) when compared to a winding optimised without the interaction of a passive magnetic shielding material (windings 501, 502, see FIGS. 5G and 5H). Winding configurations optimised to account for the interaction of the winding with a passive magnetic shielding material may therefore allow for a substantial reduction (e.g., a substantially three-fold reduction) in the volume of passive magnetic shielding material used to achieve the same region of field homogeneity. This allows the length of the cylinder (e.g., the overall length of the shield) to be reduced. Consequently, the optimisation method outlined above offers prospects for a magnetic shield which achieves high-specification performance whilst reducing size, weight and cost relative to existing commercial magnetic shields.

FIG. 7A shows an adjusted anti-Helmholtz coil arrangement inside a finite cylinder of infinite magnetic permeability optimised using the discrete element parameter optimisation method discussed above. FIG. 7A shows a pair of single-turn coils each having a radius a, separated in the z-direction by a distance of 2d₁, where solving 1.118 to remove third-order variations in the B_z field yields d₁=0.737a. FIG. 7B shows a plot of absolute magnetic field generated by the coil arrangement shown in FIG. 7A. The horizontal axis indicates a distance across the diameter of the cylinder relative to the radius a of the coils. The vertical axis indicates a distance along a longitudinal length of the cylinder relative to the radius a of the coils. The colour scale in FIG. 7B indicates values of the magnetic field strength produced by the coil arrangement in FIG. 7A. The corresponding variation in magnetic field strength throughout the space enclosed by the coils clearly illustrates a substantially linear gradient field. The dashed lines 701 and 702 indicate areas inside the space enclosed by the coils in which a deviation from linearity of the magnetic gradient field are respectively ≤5% and ≤0.1%.

FIG. 7C shows another adjusted anti-Helmholtz coil arrangement inside a finite cylinder of infinite magnetic permeability optimised using the discrete element parameter optimisation method discussed above. FIG. 7C shows a pair of single-turn coils each having a radius a, separated by a distance of 2d₁, (where d₁=0.106a), and a pair of three-turn coils each having a radius a, separated by a distance of 2d₂ (where d₂=0.887a), such that third-order and fifth-order variations in the B₂, field are removed. FIG. 7D shows a plot of absolute magnetic field generated by the coil arrangement shown in FIG. 7C, similar to the plot shown in FIG. 7B. The dashed lines 801 and 802 indicate areas inside the space enclosed by the coils in which a deviation from linearity of the magnetic gradient field are respectively ≤5% and ≤0.1%.

It will be appreciated that simple circular loops are the simplest arrangement of simple discrete winding elements, and that other coil arrangements comprising other discrete winding elements (for example, straight segments, arcs, or any other discrete configuration) may be utilised to generate different magnetic fields or correct for different orders of field variations. The parameters of the specified discrete winding elements to be optimised are one or more of a spacing (for example, along a longitudinal axis of a cylindrical structure), placement and orientation in order to produce a desired magnetic field.

It will be noted that because the optimisation of the winding configurations shown in FIGS. 7A and 7C accounts for the interaction of the windings with the passive magnetic shielding material, the performance (for example, the field homogeneity) of the windings will be noticeably or measurably worse when not in the presence of passive magnetic shielding material, as described above. The field produced by a winding optimised accounting for the interactions with passive magnetic shielding material will be more similar to a desired field (e.g., the desired field the optimisation is aimed towards producing) when the winding is in the presence of passive magnetic shielding material than when the winding is not in the presence of passive magnetic shielding material.

The methods and embodiments described above in relation to a current source restricted to flow over the inner cylindrical surface of a hollow cylinder. It will be appreciated that those methods and embodiments may be adapted for planar geometries (for example, to determine the response of a hollow cylinder of high magnetic permeability to a current source on a circular plane, parallel to planar end caps of the cylinder). That may enable one or more planar current-carrying winding configurations to be used in addition to (or alternatively to) winding configurations determined using the embodiments described above. That may lead to enhanced performance of a magnetic field cancellation (or generation) system, for example extension of a cancelling field over a larger internal region of the hollow cylindrical shield, simpler manufacturing of optimal designs, further miniaturisation of magnetic field cancellation (or generation) systems, and a greater number of options for generating magnetic fields inside magnetic shielding material (e.g., a cylindrical magnetic shield).

The embodiment described below illustrates how the method described above for a current source flowing over the inner surface of a hollow cylinder may be adapted.

As described above, accurate generation of magnetic fields in the vicinity of a material of high magnetic permeability, μ_r>>1, the boundary conditions at the surface of the material must be determined. Those conditions can be constructed through the interface conditions on the magnetic field, B, and the magnetic field strength, H. If no surface currents are present along the interface of the material with the surrounding air, the magnetic field strength tangent to the shared boundary must be continuous (as described at 1.67, where ∧ is the vector cross product). High magnetic permeability materials (for example, Mu-metal®) can have μ_r, values ≥100000 times greater than that of air. The magnetisation of those materials adjusts to provide a response that abruptly changes the magnetic field at the boundary to satisfy 1.67. That is equivalent to stating that the magnetic field components parallel to the surface of the high magnetic permeability material, S, approximately vanish (see 1.68).

The response of high permeability magnetic materials can be modelled as perfect magnetic conductors (e.g., infinite permeability). That is because the magnetic field generated by a material of high magnetic permeability in response to an induced field deviates from that of a perfect magnetic conductor on the scale of μ_r⁻¹ [8], [9].

FIG. 8 shows a closed hollow cylinder of high magnetic permeability material. The cylinder has radius ρ_s, length L_s, and planar end caps located at z=±L_s/2. The cylinder is surrounded by free space. Inside the cylinder, an arbitrary current flows around a circular coil of radius ρ_c≤ρ_s which is centred along the z-axis and located in the z=z′ plane, where |z′|<L_s/2 as shown in FIG. 8. Due to the boundary conditions at the surface of the high magnetic permeability material, the magnetic field at the boundary of the cylinder is given by

$\begin{array}{cc}{B}_{\rho }{❘}_{z=\pm L_s/2}=0,B_z{❘}_{\rho ={\rho }_s}=0,B_{\varphi }{❘}_{z=\pm L_s/2,\rho ={\rho }_s}=0& \left(1.132\right)\end{array}$

As described above, the magnetic field generated by a closed high magnetic permeability cylinder can be expressed via an infinite set of pseudo-current densities induced on the surface of an infinite cylinder. For embodiments in which a current source is restricted to flow over an inner cylindrical surface of a hollow cylinder, that may be achieved, for example, using a cylindrical expanded Green's function [10], determining matching conditions through a Fourier decomposition of the magnetic field, allowing the radial boundary condition on the surface of the cylinder to be satisfied. A similar expansion for a planar expanded Green's function in cylindrical coordinates is now shown for embodiments in which a current source flows on a circular plane.

The components of the vector potential are expressed in cylindrical coordinates

A _ρ(r)=μ₀∫_r′d³r′G(r,r′)[J_ρ(r′)cos(ϕ−ϕ′)+J_Φ(r′)sin(ϕ−ϕ′)]  (1.133)

A _ϕ(r)=−μ₀∫_r′d³r′G(r,r′)[J_ρ(r′)sin(ϕ−ϕ′)+J_Φ(r′)cos(ϕ−ϕ′)]  (1.134)

A _z(r)=μ₀∫_r′d³r′G(r,r′)J_z(r′)   (1.135)

where G(r−r′) is a Green's function satisfying Laplace's equation. The form of the vector potential can be simplified using a streamfunction representation of the current density. Since there is no current flow in the z-direction, the continuity equation can be used to represent the radial and azimuthal current densities in terms of a single scalar function

$\begin{array}{cc}{J}_{\rho }\left({\rho }^{\prime },{\varphi }^{\prime }\right)=\frac{1}{{\rho }^{\prime }}\frac{\partial S\left({\rho }^{\prime },{\varphi }^{\prime }\right)}{\partial {\varphi }^{\prime }},J_{\varphi }\left({\rho }^{\prime },{\varphi }^{\prime }\right)=-\frac{\partial S\left({\rho }^{\prime },{\varphi }^{\prime }\right)}{\partial {\rho }^{\prime }}& \left(1.136\right)\end{array}$

To exploit the radial symmetries of the system the Green's function can be decomposed through the use of the Bessel function of the first kind, J_m(z)

$\begin{array}{cc}G\left(r,r^{\prime }\right)=\frac{1}{4\pi }{\sum }_{m=-\infty }^{\infty }{e}^{\mathrm{im}\left(\varphi -{\varphi }^{\prime }\right)}{\int }_0^{\infty }dk{e}^{-k❘z-z^{\prime }❘}{J}_m\left(k\rho \right)J_m\left(k{\rho }^{\prime }\right)& \left(1.137\right)\end{array}$

allowing the vector potential to be expressed in terms of cylindrical harmonics defined on a circular plane.

Using 1.133 to 1.135, 1.136 and 1.137, the vector potential is cast in the simplified form

$\begin{array}{cc}{A}_{\rho }\left(r\right)=\frac{i{\mu }_0}{2\rho }{\sum }_{m=-\infty }^{\infty }{e}^{\mathrm{im}\varphi }{\int }_0^{\infty }{\mathrm{dkme}}^{-k❘z-z^{\prime }❘}{J}_m\left(k\rho \right)S^m\left(k\right)& \left(1.138\right)\end{array}$ $\begin{array}{cc}{A}_{\varphi }\left(r\right)=-\frac{{\mu }_0}{2}{\sum }_{m=-\infty }^{\infty }{e}^{\mathrm{im}\varphi }{\int }_0^{\infty }dkk{e}^{-k❘z-z^{\prime }❘}{J}_m\left(k\rho \right)S^m\left(k\right)& \left(1.139\right)\end{array}$ $\begin{array}{cc}{A}_z\left(r\right)=0& \left(1.14\right)\end{array}$

where J′_m(z) are the derivatives of J_m(z) with respect to z, and S^m(k) is defined by the m^th order Henkel Transform of the stream function

$\begin{array}{cc}{S}^m\left(k\right)=\frac{1}{2\pi }{\int }_0^{\infty }d{\rho }^{\prime }{\int }_0^{2\pi }d\varphi e^{-\mathrm{im}{\varphi }^{\prime }}{\rho }^{\prime }{J}_m\left(k{\rho }^{\prime }\right)S\left({\rho }^{\prime },{\varphi }^{\prime }\right)& \left(1.141\right)\end{array}$

Through the introduction of a pseudo-current density induced on an infinite cylinder we can equate the shared azimuthal Fourier modes at the boundary (1.132) and find a Fourier representation of the pseudo-current density. Since the pseudo and planar current densities are spatially orthogonal to the planar end caps, any subsequent images generated by the application of the method of mirror images continue to satisfy the radial boundary condition. This means that they satisfy the boundary conditions in 1.68 over the entire surface of the finite closed high magnetic permeability cylinder.

Taking the curl of the vector potential (1.138 and 1.139), including the magnetic field generated by a cylindrical current source as described in the embodiment above, while also applying the method of mirror images, the total magnetic field in the region ρ<ρ_s is given by

$\begin{array}{cc}{B}_{\rho }\left(\rho ,\varphi ,z\right)=-\frac{{\mu }_0}{2\pi }\sum _{p=-\infty }^{\infty }\sum _{m=-\infty }^{\infty }{e}^{\mathrm{im}\varphi }\left[{\int }_0^{\infty }dk{k}^2\frac{\left(z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s\right)}{❘z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s❘}{e}^{-k❘z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s❘}{J}_m^{\prime }\left(k\rho \right)S^m\left(k\right)-\frac{i{\rho }_s}{\pi }{\int }_{-\infty }^{\infty }dkk{e}^{ikz}{I}_m^{\prime }\left(❘k❘\rho \right)K_m^{\prime }\left(❘k❘{\rho }_s\right)J_{\varphi }^{mp}\left(k\right)\right]& \left(1.142\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }\left(\rho ,\varphi ,z\right)=-\frac{i{\mu }_0}{2\rho }\sum _{p=-\infty }^{\infty }\sum _{m=-\infty }^{\infty }m{e}^{\mathrm{im}\varphi }\left[{\int }_0^{\infty }\mathrm{dkk}\frac{\left(z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s\right)}{❘z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s❘}{e}^{-k❘z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s❘}{J}_m\left(k\rho \right)S^m\left(k\right)-\frac{i{\rho }_s}{\pi }{\int }_{-\infty }^{\infty }dk\frac{❘k❘}{k}{e}^{ikz}{I}_m\left(❘k❘\rho \right)K_m^{\prime }\left(❘k❘{\rho }_s\right)J_{\varphi }^{mp}\left(k\right)\right]& \left(1.143\right)\end{array}$ $\begin{array}{cc}{B}_z\left(\rho ,\varphi ,z\right)=\frac{{\mu }_0}{2}{\sum }_{p=-\infty }^{\infty }{\sum }_{m=-\infty }^{\infty }{e}^{\mathrm{im}\varphi }\left[{\int }_0^{\infty }dk{k}^2{e}^{-k❘z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s❘}{J}_m\left(k\rho \right)S^m\left(k\right)-\frac{{\rho }_s}{\pi }{\int }_{-\infty }^{\infty }\mathrm{dk}❘k❘e^{ikz}{I}_m\left(❘k❘\rho \right)K_m^{\prime }\left(❘k❘{\rho }_s\right)J_{\varphi }^{mp}\left(k\right)\right]& \left(1.144\right)\end{array}$

where J_ϕ^mp(k) is the p^th reflected Fourier transformed azimuthal pseudo-current density induced on the cylindrical surface of the magnetic shield material. Using the boundary condition (1.132) while imposing that each azimuthal Fourier mode of the magnetic field matches the azimuthal mode in both the planar streamfunction and the azimuthal pseudo-current density, the following relation is generated

$\begin{array}{cc}{\int }_{-\infty }^{\infty }\mathrm{dk}❘k❘e^{ikz}{I}_m\left(❘k❘{\rho }_s\right)K_m^{\prime }\left(❘k❘{\rho }_s\right)J_{\varphi }^{mp}\left(k\right)=\frac{\pi }{{\rho }_s}{\int }_0^{\infty }dd{k}^2{e}^{-k❘z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s❘}{J}_m\left(k{\rho }_s\right)S^m\left(k\right)& \left(1.145\right)\end{array}$

Physically, due to the formulation of the response in terms of a pseudo-current density, there must be a unique solution that satisfies the boundary condition over the infinite domain of the cylindrical shield which is independent of the axial position. Therefore, performing an inverse Fourier transform with respect to z results in the integral representation of the p^th reflected Fourier pseudo-current density

$\begin{array}{cc}{J}_{\varphi }^{mp}\left(k\right)=\frac{e^{-ik\left({\left(-1\right)}^{p_{z^{\prime }+\mathrm{pLs}}}\right)}}{{\rho }_s❘k❘I_m\left(❘k❘{\rho }_s\right)K_m^{\prime }\left(❘k❘{\rho }_s\right)}{\int }_0^{\infty }d\tilde{k}\frac{{\stackrel{\sim }{k}}^3}{{\stackrel{\sim }{k}}^2+k^2}{J}_m\left(\tilde{k}{\rho }_s\right)S^m\left(\tilde{k}\right)& \left(1.146\right)\end{array}$

Having determined the magnetic field in terms of the Fourier transformed streamfunction, S^m(k), an appropriate expansion of the real-space streamfunction, S(ρ′,ϕ′), must be chosen. Although the choice of orthogonal basis for the streamfunction expansion is arbitrary, a choice of basis which considers the symmetries and nature of a Hankel Transform will yield a solution that is greatly simplified. On a circular plane, a Fourier-Bessel series representation of the radial component of the streamfunction may be chosen, while using the Fourier series representation for the azimuthal dependence

$\begin{array}{cc}S\left({\rho }^{\prime },{\varphi }^{\prime }\right)=\left(H\left({\rho }^{\prime }\right)-H\left(\rho -{\rho }_c\right)\right){\rho }_c{\sum }_{n=1}^N{\sum }_{m^{\prime }=0}^M{J}_{m^{\prime }}\left(\frac{{\rho }_{n{m}^{\prime }}{\rho }^{\prime }}{{\rho }_c}\right)\left(W_{n{m}^{\prime }}\cos \left(m^{\prime }{\varphi }^{\prime }\right)+Q_{n{m}^{\prime }}\sin \left(m^{\prime }{\varphi }^{\prime }\right)\right)& \left(1.147\right)\end{array}$

where (W_nm′, Q_nm′) are Fourier coefficients and ρ_nm′ is the n^th zero of the m^th Bessel function of the first kind, J_m′(ρ_nm′)=0. The total magnetic field generated by an arbitrary current flow on the surface of a circular coil inside a closed finite cylinder of high magnetic permeability is thus given by

$\begin{array}{cc}{B}_{\rho }\left(r\right)=\frac{{\mu }_0{\rho }_c^3}{2}{\sum }_{n=1}^N{\sum }_{m=0}^M{\rho }_{nm}{J}_m^{\prime }\left({\rho }_{nm}\right)\left(W_{nm}\cos \left(m\varphi \right)+Q_{nm}\sin \left(m\varphi \right)\right)B_{\rho }^{nm}\left(\rho ,z\right)& \left(1.148\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }\left(r\right)=\frac{{\mu }_0{\rho }_c^3}{2\rho }{\sum }_{n=1}^N{\sum }_{m=0}^{M}m{\rho }_{nm}{J}_m^{\prime }\left({\rho }_{nm}\right)\left(W_{nm}\sin \left(m\varphi \right)-Q_{nm}\cos \left(m\varphi \right)\right)B_{\varphi }^{nm}\left(\rho ,z\right)& \left(1.149\right)\end{array}$ $\begin{array}{cc}{B}_z\left(r\right)=\frac{{\mu }_0{\rho }_c^3}{2}{\sum }_{n=1}^N{\sum }_{m=0}^M{\rho }_{nm}{J}_m^{\prime }\left({\rho }_{nm}\right)\left(W_{nm}\cos \left(m\varphi \right)+Q_{nm}\sin \left(m\varphi \right)\right)B_z^{nm}\left(\rho ,z\right)\mathrm{where}& \left(1.15\right)\end{array}$ $\begin{array}{cc}{B}_{\rho }^{nm}\left(\rho ,z\right)=-{\int }_0^{\infty }dk{k}^2\sigma \left(k;z,z^{\prime },L_s\right)\frac{J_m^{\prime }\left(k\rho \right)J_m\left(k{\rho }_c\right)}{k^2{\rho }_c^2-{\rho }_{nm}^2}-{\sum }_{p=1}^{\infty }\tilde{p}❘\tilde{p}❘{\lambda }_p\left(z,z^{\prime },L_s\right)\frac{I_m^{\prime }\left(❘\tilde{p}❘\rho \right)I_m\left(❘\tilde{p}❘{\rho }_c\right)K_m\left(❘\tilde{p}❘{\rho }_s\right)}{I_m\left(❘\tilde{p}❘{\rho }_s\right)\left({❘\tilde{p}❘}^2{\rho }_c^2+{\rho }_{nm}^2\right)}& \left(1.151\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }^{nm}\left(\rho ,z\right)={\int }_0^{\infty }dkk\sigma \left(k;z,z^{\prime },L_s\right)\frac{J_m\left(k\rho \right)J_m\left(k{\rho }_c\right)}{k^2{\rho }_c^2-{\rho }_{nm}^2}+{\sum }_{p=1}^{\infty }\tilde{p}{\lambda }_p\left(z,z^{\prime },L_s\right)\frac{I_m\left(❘\tilde{p}❘\rho \right)I_m\left(❘\tilde{p}❘{\rho }_c\right)K_m\left(❘\tilde{p}❘{\rho }_s\right)}{I_m\left(❘\tilde{p}❘{\rho }_s\right)\left({❘\tilde{p}❘}^2{\rho }_c^2+{\rho }_{nm}^2\right)}& \left(1.152\right)\end{array}$ $\begin{array}{cc}{B}_z^{nm}\left(\rho ,z\right)={\int }_0^{\infty }dk{k}^2\gamma \left(k;z,z^{\prime },L_s\right)\frac{J_m\left(k\rho \right)J_m\left(k{\rho }_c\right)}{k^2{\rho }_c^2-{\rho }_{nm}^2}-{\sum }_{p=1}^{\infty }{\tilde{p}}^2{\tau }_p\left(z,z^{\prime },L_s\right)\frac{I_m\left(❘\tilde{p}❘\rho \right)I_m\left(❘\tilde{p}❘{\rho }_c\right)K_m\left(❘\tilde{p}❘{\rho }_s\right)}{I_m\left(❘\tilde{p}❘{\rho }_s\right)\left({❘\tilde{p}❘}^2{\rho }_c^2+{\rho }_{nm}^2\right)}\mathrm{and}& \left(1.153\right)\end{array}$ $\begin{array}{cc}\gamma \left(k;z,z^{\prime },L_s\right)=e^{-k❘z-z^{\prime }❘}+\frac{2}{e^{2k{L}_s-1}}\left[e^{k{L}_s}\cosh \left(k\left(z+z^{\prime }\right)\right)+\cosh \left(k\left(z-z^{\prime }\right)\right)\right]& \left(1.154\right)\end{array}$ $\begin{array}{cc}\sigma \left(k;z,z^{\prime },L_s\right)=\frac{\left(z-z^{\prime }\right)}{❘z-z^{\prime }❘}{e}^{-k❘z-z^{\prime }❘}-\frac{2}{e^{2k{L}_s-1}}\left[e^{k{L}_s}\sinh \left(k\left(z+z^{\prime }\right)\right)+\sinh \left(k\left(z-z^{\prime }\right)\right)\right]& \left(1.155\right)\end{array}$ $\begin{array}{cc}{\lambda }_p\left(z,z^{\prime },L_s\right)=\frac{2}{L_s}\left({\left(-1\right)}^p\sin \left(\tilde{p}\left(z+z^{\prime }\right)\right)+\sin \left(\tilde{p}\left(z-z^{\prime }\right)\right)\right)& \left(1.156\right)\end{array}$ $\begin{array}{cc}{\tau }_p\left(z,z^{\prime },L_s\right)=\frac{2}{L_s}\left({\left(-1\right)}^p\cos \left(\tilde{p}\left(z+z^{\prime }\right)\right)+\cos \left(\tilde{p}\left(z-z^{\prime }\right)\right)\right)& \left(1.157\right)\end{array}$

with {tilde over (p)}=pπ/L_s. A derivation of the above is provided in Appendix A. Solving for the unknown Fourier coefficients, (W_nm′,Q_nm′), in the system of governing equations (1.148 to 1.150) is an ill-conditioned problem due to the formulation of the vector potential via the integral representation in 1.133 to 1.135. That may be solved using an optimisation process such as a least squares minimisation [11], with the addition of a penalty term which acts as a regularisation parameter [12]. The regularisation term may take many forms, with individual terms in it representing, for example: the curvature of a given wire geometry, the energy stored, the power dissipated etc. The choice, however, is somewhat arbitrary as any regularisation parameters acts to facilitate inversion. If the regularisation term is large, then the inverse problem is well-conditioned and yields a more simplistic coil design which compromises field fidelity. On the other hand, if the regularisation term is small, then the result is a less well-conditioned inverse problem, a more intricate coil design, and higher field fidelity. For example, using the power dissipated by the coil as a regularisation term, the power dissipated by a circular coil of thickness, t, and resistivity, , is given by

$\begin{array}{cc}P=\frac{t}{}{\int }_0^{{\rho }_c}d{\rho }^{\prime }{\rho }^{\prime }{\int }_0^{2\pi }d{\varphi }^{\prime }+{❘J_p\left({\rho }^{\prime },{\varphi }^{\prime }\right)❘}^2+{❘J_{\varphi }\left({\rho }^{\prime },{\varphi }^{\prime }\right)❘}^2& \left(1.158\right)\end{array}$

Using 1.136, 1.147 and 1.158 and integrating over the planar surface, the power can be related to the unknown Fourier coefficients. For m=0 it can be found that

$\begin{array}{cc}P=\frac{t}{}\pi {\rho }_c^2{\sum }_{n=1}^N{W}_{n0}^2{\rho }_{n0}^2{J_1\left({\rho }_{n0}\right)}^2& \left(1.159\right)\end{array}$

and for mϵ⁺ it can be found that

$\begin{array}{cc}P=\frac{t}{}\pi {\rho }_c^2\sum _{n=1}^N\sum _{m=1}^M\left(W_{nm}^2+Q_{nm}^2\right){\left(\frac{{\rho }_{nm}}{2}\right)}^{2m}\frac{1}{m!\left(m-1\right)!}\left[\text{⁠}{ }_2{\tilde{F}}_3\left(m,m+\frac{1}{2};m+1,m+1,2m+1;-{\rho }_{nm}^2\right)-\frac{{\rho }_{nm}^2}{2{\left(m+1\right)}^2}{ }_3{\tilde{F}}_4\left(m+\frac{1}{2},m+1,m+1;m,m+2,m+2,2m+1;-{\rho }_{nm}^2\right)\right]& \left(1.16\right)\end{array}$

where _i{tilde over (F)}_j is the regularised hypergeometric function (see Appendix B for more details). A cost function can now be formulated

Φ=αΣ_k=1^K[B^desired(r_k)−B(r_k)]²+βP   (1.161)

where α and β are weighting parameters chosen such that the physical parameters, such as total power, may be adjusted to fit the desired field fidelity. The cost function may be minimised, for example using the method of least squares, to calculate the optimal Fourier coefficients (W_nm′, Q_nm′) to generate a target magnetic field at K target field points. The minimisation is achieved by taking the differential of the functional with respect to the Fourier coefficients,

$\begin{array}{cc}\frac{\partial \Phi }{\partial W_{ij}}=0,\frac{\partial \Phi }{\partial Q_{ij}}=0,i>1,j>0& \left(1.162\right)\end{array}$

allowing the optimal Fourier coefficients to be found for any given physical target magnetic field profile through matrix inversion. The inversion process yields the optimal continuous streamfunction defined on the surface where current flow is permitted. The final objective, however, is to design a coil which generates the desired magnetic field to a specified accuracy. To do that, a discrete approximation to the current continuum may be found by contouring the stream function into N_s values to find the streams lines (e.g., contours of the streamfunction) where a winding should be laid to generate the desired target magnetic field.

In summary, a Green's function formulation is used to determine the magnetic field generated by an arbitrary planar current source. That may be achieved by decomposing the magnetic field into azimuthal Fourier modes. The magnetisation of the high magnetic permeability cylinder is then represented by an equivalent pseudo-current density induced on an infinite cylinder. The corresponding magnetic field is then decomposed into the same azimuthal Fourier modes. The magnetic fields generated by both the planar current source and the high magnetic permeability cylinder are then matched at the radial boundary of the high magnetic permeability cylinder (by exploiting the shared azimuthal symmetries). That may enable the pseudo-current on the cylinder to be related to the static current source on the plane. Since the planar end caps are substantially parallel to one another, and substantially perpendicular to the pseudo-current density induced on the infinite cylinder, the boundary conditions on the planar ends caps and on the cylindrical body of the high magnetic permeability material can be satisfied simultaneously (for example, using the method of mirror images).

It will be appreciated that the general framework of the method described above for embodiments in which current flow is restricted to an inner cylindrical surface of a hollow cylinder may be followed for embodiments in which current flow is restricted to inner planar surfaces of the hollow cylinder, but in formulating the magnetic field generated by the planar current source, a Green's function with symmetries akin to the Fourier modes of planar geometries is used. The coupling of the planar coil to a high magnetic permeability shield is treated in a substantially identical manner in both embodiments.

Application of the methodology described above for embodiments in which current flow is restricted to planar surfaces is discussed below. The optimised generation of a uniform magnetic field, B_x, between, and normal to, two planar coils inside a closed cylindrical magnetic shield is considered. The analytical model was verified through comparison with numerical simulations made used FEM software (COMSOL Multiphysics® Version 5.3a).

To generate B_x, first the optimal continuum currents on two circular plats which minimise the cost function were determined. Planar windings 901, 902 of radius ρ_c=0.95 m were located in top and bottom planes z′=±0.475 m and placed symmetrically inside a closed high magnetic permeability cylinder of length L_s=1 m and radius ρ_s=1m. The current continuum was generated using 800 Fourier coefficients (N=200 and M=1) on each plane (see FIGS. 9A and 9B showing the top plate and bottom plate respectively). The magnetic field was optimised between ρ=±ρ_c/2 and z=±z′/2 (shown by the dashed black lines 903, 904 in FIGS. 9C and 9D). to generate a constant target field, B_x, in the x-direction perpendicular to the planes. The colour map in FIGS. 9A and 9B represents the optimal current continuum solutions on the top (FIG. 9A) and bottom (FIG. 9B) planes, showing the magnitude of the stream function (where blue represents a low value, and red represents a high value). Solid and dashed black lines represent discrete wires with opposite sense of current flow which approximate the current continuum.

The magnetic field generated by the two plates is shown in FIG. 9C. The magnetic field B_x, generated by the windings 901, 902 only (i.e., without the high magnetic permeability cylinder present) is shown by the green dashed line 905 in FIG. 9C, and is calculated using the Biot-Savart Law taking N_s=100. The magnetic field B_x generated in the presence of the high magnetic permeability cylinder is calculated as a function of axial position z from the current continuum plots in FIGS. 9A and 9B both analytically (solid red line 906) and numerically using FEM software (dotted blue line 907) are shown in FIG. 9C. FIG. 9D shows the magnetic field deviation ΔB_x=B_x−B_x^desired of the field profiles shown in FIG. 9C, as calculated both analytically (solid red line 908) and numerically (dotted blue line 909).

It can be seen from that the magnetic field generated by the two plates in the presence of the high magnetic permeability cylinder is highly homogeneous. The maximum field deviation max(ΔB) was only 0.62% over a length of 50 cm across the optimised region, when calculated using the set of equations 1.148 to 1.150 and 0.60% when calculated numerically taking the high magnetic permeability cylinder parameters to be μ_r=20000 and thickness d=1mm. The green dashed line 905 indicates that the magnetic field B_x, generated by the windings 901, 902 is measurably worse when the windings 901, 902 are not in the presence of the high magnetic permeability cylinder, similar to the windings 601, 602 described above.

From reading the present disclosure, other variations and modifications will be apparent to the skilled person. Such variations and modifications may involve equivalent and other features which are already known in the art of magnetic shielding technology, and which may be used instead of, or in addition to, features already described herein.

Although the appended claims are directed to particular combinations of features, it should be understood that the scope of the disclosure of the present invention also includes any novel feature or any novel combination of features disclosed herein either explicitly or implicitly or any generalisation thereof, whether or not it relates to the same invention as presently claimed in any claim and whether or not it mitigates any or all of the same technical problems as does the present invention.

Features which are described in the context of separate embodiments may also be provided in combination in a single embodiment. Conversely, various features which are, for brevity, described in the context of a single embodiment, may also be provided separately or in any suitable sub-combination. The applicant hereby gives notice that new claims may be formulated to such features and/or combinations of such features during the prosecution of the present application or of any further application derived therefrom.

For the sake of completeness, it is also stated that the term “comprising” does not exclude other elements or steps, the term “a” or “an” does not exclude a plurality, a single processor or other unit may fulfil the functions of several means recited in the claims and any reference signs in the claims shall not be construed as limiting the scope of the claims.

APPENDIX A

In order to derive the magnetic field in terms of the streamfunction, first both the Fourier transform of the streamfunction and the induced Fourier pseudo-current density must be determined. Substituting 1.147 into 1.141 and integrating over the azimuthal coordinate yields

$$\begin{gathered} \begin{array}{cc}{S}^m\left(k\right)=\frac{{\rho }_c}{2}{\sum }_{n=1}^N{\sum }_{m^{\prime }=0}^M\left(W_{{\mathrm{nm}}^{\prime }}\left({\delta }_{{\mathrm{mm}}^{\prime }}+{\delta }_{m,-m^{\prime }}\right)+i{Q}_{{\mathrm{nm}}^{\prime }}\left({\delta }_{{\mathrm{mm}}^{\prime }}-{\delta }_{m,-m^{\prime }}\right)\right)\times \text{ \\ ∫0ρcd⁢ρ′⁢ρ′⁢Jm(k⁢ρ′)⁢Jm′(ρnm′⁢ρ′ρc)}& \left(2.1\right)\end{array} \end{gathered}$$

Upon substituting 2.1 into the magnetic field 1.132 to 1.134 the delta functions collapse the infinite summation, resulting in every Bessel function of order m being replaced by m′. Hence, without loss of generality, using the Bessel function product identity found in [12], the following may be written

$$\begin{gathered} \begin{array}{cc}{S}^m\left(k\right)=\frac{{\rho }_c^3}{2}{\sum }_{n=1}^N{\sum }_{m^{\prime }=0}^M\left(W_{{\mathrm{nm}}^{\prime }}\left({\delta }_{{\mathrm{mm}}^{\prime }}+{\delta }_{m,-m^{\prime }}\right)+i{Q}_{{\mathrm{nm}}^{\prime }}\left({\delta }_{{\mathrm{mm}}^{\prime }}-\text{ \\ δm,-m′}\right)\right)\times \frac{{\rho }_{{\mathrm{nm}}^{\prime }}}{k^2{\rho }_c^2-{\rho }_{{\mathrm{nm}}^{\prime }}^2}{J}_{m^{\prime }}\left(k{\rho }_c\right)J_{m^{\prime }}^{\prime }\left({\rho }_{{\mathrm{nm}}^{\prime }}\right)& \left(2.2\right)\end{array} \end{gathered}$$

Similarly, the induced Fourier pseudo-current density may be now be calculated. Substituting 1.147 into 1.146 and integrating over the azimuthal coordinate yields

$$\begin{gathered} \begin{array}{cc}{J}_{\varphi }^{mp}\left(k\right)=\frac{{\rho }_c{e}^{-ik\left({\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s\right)}}{2{\rho }_s❘k❘I_m\left(❘k❘{\rho }_s\right)K_m^{\prime }\left(❘k❘{\rho }_s\right)}{\sum }_{n=1}^N{\sum }_{m^{\prime }=0}^M\left(W_{{\mathrm{nm}}^{\prime }}\left({\delta }_{{\mathrm{mm}}^{\prime }}+{\delta }_{m,-m^{\prime }}\right)+{\mathrm{iQ}}_{{\mathrm{nm}}^{\prime }}\left({\delta }_{{\mathrm{mm}}^{\prime }}-{\delta }_{m,-m^{\prime }}\right)\right)\times \text{ \\ ∫0∞d⁢k˜⁢k~3k~2+k2⁢Jm(k˜⁢ρs)⁢∫0ρcd⁢ρ′⁢ρ′⁢Jm′(k˜⁢ρ′)⁢Jm′(ρnm′⁢ρ′ρc)}& \left(2.3\right)\end{array} \end{gathered}$$

First the integral over {tilde over (k)} is considered. Separating this using partial fractions yields

$\begin{array}{cc}{\int }_0^{\infty }d\tilde{k}\frac{{\tilde{k}}^3}{{\tilde{k}}^2+k^2}{J}_m\left(\tilde{k}{\rho }_s\right)J_m\left(\tilde{k}{\rho }^{\prime }\right)={\int }_0^{\infty }d\tilde{k}\tilde{k}{J}_m\left(\tilde{k}{\rho }_s\right)J_m\left(\tilde{k}{\rho }^{\prime }\right)-k^2{\int }_0^{\infty }d\tilde{k}\frac{\tilde{k}}{{\tilde{k}}^2+k^2}{J}_m\left(\tilde{k}{\rho }_s\right)J_m\left(\tilde{k}{\rho }^{\prime }\right)& \left(2.4\right)\end{array}$

The first integral is the orthogonality relation between Bessel functions, and the second can be found in [11]. It can therefore be stated that

$\begin{array}{cc}{\int }_0^{\infty }d\tilde{k}\frac{{\tilde{k}}^3}{{\tilde{k}}^2+k^2}{J}_m\left(\tilde{k}{\rho }_s\right)J_m\left(\tilde{k}{\rho }^{\prime }\right)=\frac{\delta \left({\rho }^{\prime }-{\rho }_s\right)}{{\rho }^{\prime }}-k^2{I}_{m^{\prime }}\left(❘k❘{\rho }^{\prime }\right)K_{m^{\prime }}\left(❘k❘{\rho }_s\right)& \left(2.5\right)\end{array}$

Integrating 2.5 into 2.3 and integrating over ρ′ using the same identity as in 2.2 results in the expression

$\begin{array}{cc}{J}_{\varphi }^{mp}=\frac{{\rho }_c{k}^2{e}^{-ik\left({\left(-1\right)}^p{z}^{\prime }+\mathrm{pLs}\right)}}{2{\rho }_s❘k❘I_m\left(❘k❘{\rho }_s\right)K_m^{\prime }\left(❘k❘{\rho }_s\right)}{\sum }_{n=1}^N{\sum }_{m^{\prime }=0}^M\left(W_{{\mathrm{nm}}^{\prime }}\left({\delta }_{{\mathrm{mm}}^{\prime }}+{\delta }_{m,-m^{\prime }}\right)+i{Q}_{{\mathrm{nm}}^{\prime }}\left({\delta }_{{\mathrm{mm}}^{\prime }}-{\delta }_{m,-m^{\prime }}\right)\right)\times \frac{{\rho }_{{\mathrm{nm}}^{\prime }}}{k^2{\rho }_c^2+{\rho }_{{\mathrm{nm}}^{\prime }}^2}{J}_{m^{\prime }}^{\prime }\left({\rho }_{{\mathrm{nm}}^{\prime }}\right)I_{m^{\prime }}\left(❘k❘{\rho }_c\right)K_{m^{\prime }}\left(❘k❘{\rho }_s\right)& \left(2.6\right)\end{array}$

Substituting 2.2 and 2.6 into 1.142 to 1.144 allows the expression of the summation over exponentials associated with the infinite reflections of the planar streamfunction as

$\begin{array}{cc}{\sum }_{p=-\infty }^{\infty }{e}^{-k❘z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s❘}=e^{-k❘z-z^{\prime }❘}+\frac{2}{e^{2kL}-1}\left(e^{kL}\cosh \left(k\left(z+z^{\prime }\right)\right)+\cosh \left(k\left(z-z^{\prime }\right)\right)\right)& \left(2.7\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}{\sum }_{p=-\infty }^{\infty }\frac{\left(z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s\right)}{❘z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s❘}{e}^{-k❘z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s❘}=\frac{\left(z-z^{\prime }\right)}{❘z-z^{\prime }❘}{e}^{-k\left|z-z^{\prime }\right|}-\frac{2}{e^{2kL}-1}\left(e^{kL}\sinh \left(k\left(z+z^{\prime }\right)\right)+\sinh \left(k\left(z-z^{\prime }\right)\right)\right)& \left(2.8\right)\end{array}$

while expressing the summation over the infinite pseudo-current reflections through a Fourier series expansion

$\begin{array}{cc}{\sum }_{p=-\infty }^{\infty }{e}^{ik\left(z-{\left(-1\right)}^p{z}^{\prime }+{\mathrm{pL}}_s\right)}=\frac{\pi }{L}{\sum }_{p=-\infty }^{\infty }\delta \left(k-\frac{p\pi }{L}\right)\left(e^{ik\left(z+z^{\prime }-L\right)}+e^{ik\left(z-z^{\prime }\right)}\right)& \left(2.9\right)\end{array}$

It may finally be written

$\begin{array}{cc}{B}_{\rho }\left(r\right)=\frac{{\mu }_0{\rho }_c^3}{2\rho }{\sum }_{n=1}^N{\sum }_{m=0}^M{\rho }_{nm}{J}_m^{\prime }\left({\rho }_{nm}\right)\left(W_{nm}\cos \left(m\varphi \right)+Q_{nm}\sin \left(m\varphi \right)\right)B_{\rho }^{nm}\left(\rho ,z\right)& \left(2.1\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }\left(r\right)=\frac{{\mu }_0{\rho }_c^3}{2\rho }{\sum }_{n=1}^N{\sum }_{m=0}^{M}m{\rho }_{nm}{J}_m^{\prime }\left({\rho }_{nm}\right)\left(W_{nm}\sin \left(m\varphi \right)-Q_{nm}\cos \left(m\varphi \right)\right)B_{\varphi }^{nm}\left(\rho ,z\right)& \left(2.11\right)\end{array}$ $\begin{array}{cc}{B}_z\left(r\right)=\frac{{\mu }_0{\rho }_c^3}{2}{\sum }_{n=1}^N{\sum }_{m=0}^M{\rho }_{nm}{J}_m^{\prime }\left({\rho }_{nm}\right)\left(W_{nm}\cos \left(m\varphi \right)+Q_{nm}\sin \left(m\varphi \right)\right)B_z^{nm}\left(\rho ,z\right)& \left(2.12\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{B}_{\rho }^{nm}\left(\rho ,z\right)=-{\int }_0^{\infty }dk{k}^2\sigma \left(k;z,z^{\prime },L_s\right)\frac{J_m^{\prime }\left(k\rho \right)J_m\left(k{\rho }_c\right)}{k^2{\rho }_c^2-{\rho }_{nm}^2}-{\sum }_{p=1}^{\infty }\tilde{p}❘\tilde{p}❘{\lambda }_p\left(z,z^{\prime },L_S\right)\frac{I_m^{\prime }\left(❘\tilde{p}❘\rho \right)I_m\left(❘\tilde{p}❘{\rho }_c\right)K_m\left(❘\tilde{p}❘{\rho }_s\right)}{I_m\left(❘\tilde{p}❘{\rho }_s\right)\left({❘\tilde{p}❘}^2{\rho }_c^2+{\rho }_{nm}^2\right)}& \left(2.13\right)\end{array}$ $\begin{array}{cc}{B}_{\varphi }^{nm}\left(\rho ,z\right)={\int }_0^{\infty }dkk\sigma \left(k;z,z^{\prime },L_s\right)\frac{J_m\left(k\rho \right)J_m\left(k{\rho }_c\right)}{k^2{\rho }_c^2-{\rho }_{nm}^2}+{\sum }_{p=1}^{\infty }\tilde{p}{\lambda }_p\left(z,z^{\prime },L_s\right)\frac{I_m\left(❘\tilde{p}❘\rho \right)I_m\left(❘\tilde{p}❘{\rho }_c\right)K_m\left(❘\tilde{p}❘{\rho }_s\right)}{I_m\left(❘\tilde{p}❘{\rho }_s\right)\left({❘\tilde{p}❘}^2{\rho }_c^2+{\rho }_{nm}^2\right)}& \left(2.14\right)\end{array}$ $\begin{array}{cc}{B}_z^{nm}\left(\rho ,z\right)={\int }_0^{\infty }dk{k}^2\gamma \left(k;z,z^{\prime },L_s\right)\frac{J_m\left(k\rho \right)J_m\left(k{\rho }_c\right)}{k^2{\rho }_c^2-{\rho }_{nm}^2}-{\sum }_{p=1}^{\infty }{\tilde{p}}^2{\tau }_p\left(z,z^{\prime },L_s\right)\frac{I_m\left(❘\tilde{p}❘\rho \right)I_m\left(❘\tilde{p}❘{\rho }_c\right)K_m\left(❘\tilde{p}❘{\rho }_s\right)}{I_m\left(❘\tilde{p}❘{\rho }_s\right)\left({❘\tilde{p}❘}^2{\rho }_c^2+{\rho }_{nm}^2\right)}& \left(2.15\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}\gamma \left(k;z,z^{\prime },L_s\right)=e^{-k❘z-z^{\prime }❘}+\frac{2}{e^{2k{L}_s}-1}\left[e^{k{L}_s}\cosh \left(k\left(z+z^{\prime }\right)\right)+\cosh \left(k\left(z-z^{\prime }\right)\right)\right]& \left(2.16\right)\end{array}$ $\begin{array}{cc}\sigma \left(k;z,z^{\prime },L_s\right)=\frac{\left(z-z^{\prime }\right)}{❘z-z^{\prime }❘}{e}^{-k❘z-z^{\prime }|}-\frac{2}{e^{2k{L}_s}-1}\left[e^{k{L}_s}\sinh \left(k\left(z+z^{\prime }\right)\right)+\sinh \left(k\left(z-z^{\prime }\right)\right)\right]& \left(2.17\right)\end{array}$ $\begin{array}{cc}{\lambda }_p\left(z,z^{\prime },L_s\right)=\frac{2}{L_s}\left({\left(-1\right)}^p\sin \left(\tilde{p}\left(z+z^{\prime }\right)\right)+\sin \left(\tilde{p}\left(z-z^{\prime }\right)\right)\right)& \left(2.18\right)\end{array}$ $\begin{array}{cc}{\tau }_p\left(z,z^{\prime },L_s\right)=\frac{2}{L_s}\left({\left(-1\right)}^p\cos \left(\tilde{p}\left(z+z^{\prime }\right)\right)+\cos \left(\tilde{p}\left(z-z^{\prime }\right)\right)\right)& \left(2.19\right)\end{array}$ $\mathrm{and}\tilde{p}=p\pi /L_s.$

APPENDIX B

The power dissipated in the surface of a coil of thickness, t, and resistivity, , is given by

$\begin{array}{cc}P=\frac{r}{}{\int }_0^{{\rho }_c}d{\rho }^{\prime }{\rho }^{\prime }{\int }_0^{2\pi }d{\varphi }^{\prime }{❘J_{\rho }\left({\rho }^{\prime },{\varphi }^{\prime }\right)❘}^2+{❘J_{\varphi }({\rho }^{\prime },{\varphi }^{\prime }❘}^2& \left(2.2\right)\end{array}$

Substituting the streamfunction 1.147 into the continuity relations 1.136 it is found that

$\begin{array}{cc}{J}_{\rho }\left({\rho }^{\prime },{\varphi }^{\prime }\right)=\left(H\left({\rho }^{\prime }-{\rho }_c\right)-H\left({\rho }^{\prime }\right)\right)\frac{{\rho }_c}{{\rho }^{\prime }}{\sum }_{n=1}^N{\sum }_{m=0}^{M}m{J}_m\left(\frac{{\rho }_{nm}{\rho }^{\prime }}{{\rho }_c}\right)\left(W_{nm}\sin \left(m{\varphi }^{\prime }\right)-Q_{nm}\cos \left(m{\varphi }^{\prime }\right)\right)& \left(2.21\right)\end{array}$ $\begin{array}{cc}{J}_{\varphi }\left({\rho }^{\prime },{\varphi }^{\prime }\right)=\left(H\left({\rho }^{\prime }-{\rho }_c\right)-H\left({\rho }^{\prime }\right)\right){\sum }_{n=1}^N{\sum }_{m=0}^M{\rho }_{nn}{J}_m^{\prime }\left(\frac{{\rho }_{nm}{\rho }^{\prime }}{{\rho }_c}\right)\left(W_{nm}\cos \left(m{\varphi }^{\prime }\right)+Q_{nm}\sin \left(m{\varphi }^{\prime }\right)\right)& \left(2.22\right))\end{array}$

Inserting 2.21 and 2.22 into 2.20 and integrating over the azimuthal component results in

$\begin{array}{cc}P=\frac{r}{}\pi {\rho }_c^2{\sum }_{n=1}^N{\sum }_{m=0}^M\left[\left(1+{\delta }_{m0}\right)W_{nm}^2+\left(1-{\delta }_{m0}\right)Q_{nm}^2\right]\times {\int }_0^{1}d\tilde{\rho }\left[\frac{2m^2}{\tilde{\rho }}{J_m\left({\rho }_{nm}\tilde{\rho }\right)}^2-{\rho }_{nm}^2\tilde{\rho }{J}_{m-1}\left({\rho }_{nm}\tilde{\rho }\right)J_{m+1}\left({\rho }_{nm}\tilde{\rho }\right)\right]& \left(2.23\right)\end{array}$

for m∈⁺, where _i{tilde over (F)}_j is the regularised hypergeometric function. Some specific evaluations of 2.23 are now presented:

$\begin{array}{cc}\begin{array}{cc}m=1& P=\frac{2t}{}\pi {\rho }_c^2{\sum }_{n=1}^N\left(W_{n1}^2+Q_{n1}^2\right){\rho }_{n1}^2{J_0\left({\rho }_{n1}\right)}^2\end{array}& \left(2.24\right)\end{array}$ $\begin{array}{cc}m=2P=\frac{2t}{}\pi {\rho }_c^2{\sum }_{n=1}^N\left(W_{n2}^2+Q_{n2}^2\right)\left[\left({\rho }_{n2}^2-4\right){J_0\left({\rho }_{n2}\right)}^2-\frac{2\left({\rho }_{n2}^2-8\right)}{{\rho }_{n2}}{J}_0\left({\rho }_{n2}\right)J_1\left({\rho }_{n2}\right)+\frac{\left({\rho }_{n2}^4-16\right)}{{\rho }_{n2}^2}{J}_1{\left({\rho }_{n2}\right)}^2\right]& \left(2.25\right)\end{array}$ $$\begin{gathered} \begin{array}{cc}m=3P=\frac{2t}{}\pi {\rho }_c^2{\sum }_{n=1}^N\left(W_{n3}^2+Q_{n3}^2\right)\left[\frac{\left({\rho }_{n3}^4-96\right)}{{\rho }_{n3}^2}{J_0\left({\rho }_{n3}\right)}^2-\text{ \\ 4⁢8⁢(ρn⁢32-8)ρn⁢33⁢J0(ρn⁢3)⁢J1(ρn⁢3)+(ρn⁢36-1⁢0⁢ρn⁢34+9⁢6⁢ρn⁢32-3⁢8⁢4)ρn⁢34⁢J1(ρn⁢3)2}\right]& \left(2.26\right)\end{array} \end{gathered}$$

REFERENCES

• [1] Elena Boto et al. Moving magnetoencephalography towards real-world applications with a wearable system Nature, 555 (657-661), 2018
• [2] J. D. Jackson Classical Electrodynamics (3^rd ed.) Wiley, New York, 1998
• [3] R. Turner and R. M. Bowley Passive screening of switched magnetic field gradients Journal of Physics E: Scientific Instruments, 19(876), 1986
• [4] R. Turner A target field approach to optimal coil design Journal of Physics: Applied Physics, 19(8), 1986
• [5] K. W. Rigby Design of magnets inside cylindrical superconductive shields Review of Scientific Instruments, 59(156), 1988
• [6] C. P. Bidinosti, J. W. Martin Passive magnetic shielding in static gradient fields AIP Advances 4(047135), 2014
• [7] C.-Y. Liu, Y. Andalib, D. C. M. Ostapchuk, C. P. Bidinosti Analytic models of magnetically enclosed spherical and solenoidal coils 2019
• [8] P. Hammond Electric and magnetic images Proceedings of the IEE—Part C: Monographs, 107:306, 1960
• [9] P. Hammond Effect of finite thickness of magnetic substrate on planar inductors IEEE Transactions on

Magnetics, 26:270-275, 1990

• [10] L.-W. Li, M.-S. Leong, T.-S. Yeo, and P.-S. Kooi Electromagnetic dyadic green 's functions in spectral domain for multi-layered cylinders Journal of Electromagnetic Waves and Applications, 14:961-985, 2000
• [11] Otto Bretscher Linear Algebra with Applications Pearson, London, 4 edition, 2009
• [12] I. S. Gradshteyn and I. M. Ryzhik Table of Integrals, Series, and Products. Elsevier/Academic Press, Amsterdam, seventh edition, 2007. Translated from the Russian, Translation edited and with a preface by

Alan Jeffrey and Daniel Zwillinger

Claims

1. A method of designing a magnetic shield comprising a structure enclosing a space, the structure comprising passive magnetic shielding material, and a winding configured to produce a specified magnetic field within the structure when current is passed through the winding, the method comprising: determining an optimised configuration of the winding accounting for the presence of the passive magnetic shielding material by implementing one or more boundary conditions at the surface of the passive magnetic shielding material.

2. The method of claim 1, wherein the one or more boundary conditions require that the magnetic field produced by the winding is zero on a surface of the passive magnetic shielding material.

3. The method of claim 1 or of claim 2, wherein accounting for the presence of the passive magnetic shielding material comprises constructing a function, or discrete approximation thereof, for a geometry of the structure, that can be used to solve differential equations relating magnetic field to current density.

4. The method of claim 3, wherein the function is a Green's function subject to one or more Dirichlet boundary conditions, and further comprising implementing the method of mirror images.

5. The method of any preceding claim, wherein determining an optimised configuration of the winding comprises implementing an optimisation process.

6. The method of claim 5, wherein the optimisation process comprises a least squares minimisation process, and optionally further comprises regularising the least squares minimisation with a penalty term, and further optionally wherein the penalty term comprises at least one of a power consumption of the winding, a curvature of the winding configuration, a resistance, an inductance of the winding, a mass and/or volume of the winding, and an energy stored in the winding.

7. The method of claim 5, wherein the optimisation process comprises solving a set of simultaneous equations.

8. The method of any preceding claim, wherein determining an optimised configuration of the winding comprises determining optimal surface currents on the structure required to produce the specified magnetic field within the structure.

9. The method of claim 8, wherein the surface currents are defined by a streamfunction, and determining the optimal surface currents comprises determining streamlines of the streamfunction.

10. The method of claim 9, wherein determining streamlines of the streamfunction comprises discretising the streamfunction to determine contours of the streamfunction.

11. The method of any preceding claim, wherein determining an optimised configuration of the winding comprises: defining at least one discrete winding element having at least one free parameter; andoptimising the at least one free parameter of each of the discrete winding elements to produce the specified magnetic field within the structure.

12. The method of claim 11, wherein optimising the at least one free parameter of the discrete winding elements comprises: determining expressions for components of current density for each of the discrete winding elements;substituting the expressions for components of current density into expressions for components of the magnetic field produced by an arbitrary winding configuration; andimplementing an optimisation process to determine, for the specified field, optimal values relating to the components of current density.

13. The method of claim 12, wherein the optimal values are or include an optimal relationship between free parameters of the discrete winding elements.

14. The method of any of claims 11 to 13, wherein the at least one free parameter of each discrete winding element comprises at least one of a size, a shape, a spacing, a placement and an orientation of the discrete winding element.

15. The method of any preceding claim, wherein the magnetic field is one of a constant magnetic field, a linear gradient magnetic field and a higher-order magnetic field.

16. The method of any preceding claim, wherein the structure is a closed structure.

17. The method of any preceding claim, wherein the structure is a hollow cylinder.

18. The method of claim 17 as it depends from claim 16, wherein the hollow cylinder comprises an axial body and planar end surfaces.

19. The method of claim 18, wherein separate boundary conditions are implemented for an axial surface of the cylinder and planar end surfaces of the cylinder.

20. The method of any of claims 1 to 19, further comprising manufacturing a winding corresponding to or approximating the determined optimised configuration.

21. A magnetic shield designed using the method of any of claims 1 to 20.

22. A magnetic shield comprising: a structure enclosing a space, the structure comprising passive magnetic shielding material; anda winding configured to produce a specified magnetic field within the structure when current is passed through the winding;wherein a configuration of the winding is optimised such that a field produced by the winding is more similar to the specified magnetic field when the winding is in the presence of the structure than when the winding is not in the presence of the structure.

23. The magnetic shield of claim 22, wherein the structure is a closed structure.

24. The magnetic shield of claim 22 or of claim 23, wherein the structure is a hollow cylinder.

25. The magnetic shield of any of claims 22 to 24, wherein the winding is disposed within the space enclosed by the structure, and optionally wherein the winding is disposed on an interior surface of the structure.

26. The magnetic shield of any of claims 22 to 25, wherein the winding, in the presence of the structure is configured to produce a uniform magnetic field over a volume of at least 1000 cm3 with a relative field variation of less than 1%.