From 927bbf75a40ba1661c6cd2d7ce0bbdbc49da5268 Mon Sep 17 00:00:00 2001
From: Jonathan Matthews <jonathan.matthews@ukaea.uk>
Date: Wed, 9 Oct 2024 15:56:30 +0100
Subject: [PATCH 1/7] added Greens function to PFcoils for field calculations
 outside the coils

---
 process/pfcoil.py | 206 +++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 185 insertions(+), 21 deletions(-)

diff --git a/process/pfcoil.py b/process/pfcoil.py
index 2b62b0de63..93aea99b7e 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -24,6 +24,7 @@
 from process.fortran import tfcoil_variables as tfv
 from process.fortran import times_variables as tv
 from process.utilities.f2py_string_patch import f2py_compatible_to_string
+from scipy.special import ellipe, ellipk
 
 logger = logging.getLogger(__name__)
 
@@ -3232,8 +3233,8 @@ def bfield(rc, zc, cc, rp, zp):
         t = 1.0 - s
         a = np.log(1.0 / t)
 
-        dz = zp - zc[i]
-        zs = dz**2
+        # dz = zp - zc[i]
+        # zs = dz**2
         dr = rp - rc[i]
         sd = np.sqrt(d)
 
@@ -3260,25 +3261,28 @@ def bfield(rc, zc, cc, rp, zp):
 
         #  Radial, vertical fields
 
-        brx = (
-            RMU0
-            * cc[i]
-            * dz
-            / (2 * np.pi * rp * sd)
-            * (-xk + (rc[i] ** 2 + rp**2 + zs) / (dr**2 + zs) * xe)
-        )
-        bzx = (
-            RMU0
-            * cc[i]
-            / (2 * np.pi * sd)
-            * (xk + (rc[i] ** 2 - rp**2 - zs) / (dr**2 + zs) * xe)
-        )
-
-        #  Sum fields, flux
-
-        br += brx
-        bz += bzx
-        psi += xc[i] * cc[i]
+        # brx = (
+        #     RMU0
+        #     * cc[i]
+        #     * dz
+        #     / (2 * np.pi * rp * sd)
+        #     * (-xk + (rc[i] ** 2 + rp**2 + zs) / (dr**2 + zs) * xe)
+        # )
+        # bzx = (
+        #     RMU0
+        #     * cc[i]
+        #     / (2 * np.pi * sd)
+        #     * (xk + (rc[i] ** 2 - rp**2 - zs) / (dr**2 + zs) * xe)
+        # )
+        # get the field due to unit current using Green's function
+        g_br, g_bz, g_psi = greens(rc[i], zc[i], rp, zp)
+        # convert to actual field from unit current
+
+        #  Sum fields, flux and convert from unit current to actual
+
+        br += g_br * cc[i]
+        bz += g_bz * cc[i]
+        psi += g_psi * cc[i]
 
     return xc, br, bz, psi
 
@@ -3615,3 +3619,163 @@ def init_pfcoil_variables():
     pfv.d_cond_cst = 0.0
     pfv.r_in_cst = 0.0
     pfv.r_out_cst = 3.0e-3
+
+
+GREENS_ZERO = 1e-8
+
+
+@numba.vectorize(nopython=True, cache=True)
+def clip_nb(
+    val: float | np.ndarray, val_min: float, val_max: float
+) -> float | np.ndarray:
+    """
+    Clips (limits) val between val_min and val_max. Vectorised for speed and
+    compatibility with numba.
+
+    Parameters
+    ----------
+    val:
+        The value to be clipped.
+    val_min:
+        The minimum value.
+    val_max:
+        The maximum value.
+
+    Returns
+    -------
+    The clipped values.
+    """
+    if val < val_min:
+        return val_min
+    if val > val_max:
+        return val_max
+    return val
+
+
+@numba.vectorize(nopython=True)
+def ellipe_nb(k: float | np.ndarray) -> float | np.ndarray:
+    """
+    Vectorised scipy ellipe
+
+    Notes
+    -----
+    K, E in scipy are set as K(k^2), E(k^2)
+    """
+    return ellipe(k)
+
+
+ellipe_nb.__doc__ += ellipe.__doc__
+
+
+@numba.vectorize(nopython=True)
+def ellipk_nb(k: float | np.ndarray) -> float | np.ndarray:
+    """
+    Vectorised scipy ellipk
+
+    Notes
+    -----
+    K, E in scipy are set as K(k^2), E(k^2)
+    """
+    return ellipk(k)
+
+
+@numba.jit(nopython=True)
+def calc_a_k2(
+    xc,
+    zc,
+    x,
+    z,
+):
+    a = np.hypot((x + xc), (z - zc))
+    k2 = 4 * x * xc / a**2
+    # Avoid NaN when coil on grid point
+    k2 = clip_nb(k2, GREENS_ZERO, 1.0 - GREENS_ZERO)
+    return a, k2
+
+
+@numba.jit(nopython=True)
+def calc_e_k(
+    k2,
+):
+    return ellipe_nb(k2), ellipk_nb(k2)
+
+
+@numba.jit(nopython=True)
+def calc_i1_i2(
+    a,
+    k2,
+    e,
+    k,
+):
+    if (e is None) or (k is None):
+        e, k = calc_e_k(k2)
+    i1 = k / a
+    i2 = e / (a**3 * (1 - k2))
+    return i1, i2
+
+
+@numba.jit(nopython=True)
+def greens(rc, zc, rp, zp):
+    """
+    Green's function for psi, Bx and Bz. Output is from unit current so for actual
+    field multiply result by coil current.
+
+    rc: input real array : coil r coordinates (m)
+    zc: input real array : coil z coordinates (m)
+    rp: input real array : calculation r coordinates (m)
+    zp: input real array : calculation z coordinates (m)
+    g_psi : output real : poloidal flux at (rp,zp) (Wb)
+    g_br: output real : radial field component at (rp,zp) (T)
+    g_bz : output real : vertical field component at (rp,zp) (T)
+    """
+
+    a, k2 = calc_a_k2(rc, zc, rp, zp)
+    e, k = calc_e_k(k2)
+    i1, i2 = calc_i1_i2(a, k2, e, k)
+    i1 *= 4
+    i2 *= 4
+
+    a_part = (zp - zc) ** 2 + rp**2 + rc**2
+    b_part = -2 * rp * rc
+
+    g_br = 0.25 * RMU0 / np.pi * rc * (zp - zc) * (i1 - i2 * a_part) / b_part
+    g_bz = (
+        0.25 * RMU0 / np.pi * rc * ((rc + rp * a_part / b_part) * i2 - i1 * rp / b_part)
+    )
+    g_psi = 0.25 * RMU0 / np.pi * a * ((2 - k2) * k - 2 * e)
+
+    return g_psi, g_br, g_bz
+
+
+@numba.jit(nopython=True)
+def elliptic_inductance(radius, rc):
+    """
+    Calculate the inductance of a circular coil by elliptic integrals.
+
+    Parameters
+    ----------
+    radius: input real : The radius of the circular coil
+    rc: input real : The radius of the coil cross-section
+    xc: output real : The self-inductance of the circular coil [H]
+
+    Notes
+    -----
+    The inductance is given by
+
+    .. math::
+        L = \\mu_{0} (2 r - r_c) \\Biggl((1 - k^2 / 2)~
+        \\int_0^{\\frac{\\pi}{2}} \\frac{d\\theta}{\\sqrt{1 - k~
+        \\sin (\\theta)^2}} - \\int_0^{\\frac{\\pi}{2}}~
+        \\sqrt{1 - k \\sin (\\theta)^2} \\, d\\theta\\Biggr)
+
+    where :math:`r` is the radius, :math:`\\mu_{0}` is the vacuum
+    permeability, and
+
+    .. math::
+        k = \\max\\left(10^{-8}, \\min~
+        \\left(\\frac{4r(r - r_c)}{(2r - r_c)^2}~
+        , 1.0 - 10^{-8}\\right)\\right)
+    """
+    k = 4 * radius * (radius - rc) / (2 * radius - rc) ** 2
+    k = clip_nb(k, GREENS_ZERO, 1.0 - GREENS_ZERO)
+    return RMU0 * (2 * radius - rc) * ((1 - k**2 / 2) * ellipk_nb(k) - ellipe_nb(k))

From ea940aa8864fc7094feabac520049d637734a1f9 Mon Sep 17 00:00:00 2001
From: Jonathan Matthews <jonathan.matthews@ukaea.uk>
Date: Tue, 22 Oct 2024 16:56:11 +0100
Subject: [PATCH 2/7] semianalytic function for pfcoils

---
 process/pfcoil.py | 390 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 390 insertions(+)

diff --git a/process/pfcoil.py b/process/pfcoil.py
index 93aea99b7e..e1ef1af647 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -24,6 +24,16 @@
 from process.fortran import tfcoil_variables as tfv
 from process.fortran import times_variables as tv
 from process.utilities.f2py_string_patch import f2py_compatible_to_string
+from process import fortran as ft
+import process.superconductors as superconductors
+import math
+import numpy as np
+import numba
+from numba.types import CPointer, float64, intc
+from scipy import LowLevelCallable
+from scipy.integrate import IntegrationWarning, quad, nquad
+import logging
+from scipy import optimize
 from scipy.special import ellipe, ellipk
 
 logger = logging.getLogger(__name__)
@@ -3720,6 +3730,9 @@ def greens(rc, zc, rp, zp):
     Green's function for psi, Bx and Bz. Output is from unit current so for actual
     field multiply result by coil current.
 
+    Used for producing the radial (x) and vertical (z) components of the magnetic field,
+    and the magnetic flux at any point not inside a solenoid (PF or CS).
+
     rc: input real array : coil r coordinates (m)
     zc: input real array : coil z coordinates (m)
     rp: input real array : calculation r coordinates (m)
@@ -3747,6 +3760,383 @@ def greens(rc, zc, rp, zp):
     return g_psi, g_br, g_bz
 
 
+def floatify(x):
+    """
+    Converts the np array or float into a float by returning
+    the first element or the element itself.
+
+    Notes
+    -----
+    This function aims to avoid numpy warnings for float(x) for >0 rank scalars
+    it emulates the functionality of float conversion
+
+    Raises
+    ------
+    ValueError
+        If array like object has more than 1 element
+    TypeError
+        If object is None
+    """
+    if x is None:
+        raise TypeError("The argument cannot be None")
+    return np.asarray(x, dtype=float).item()
+
+
+def integrate(func, args, bound1, bound2):
+    """
+    Utility for integration of a function between bounds. Easier to refactor
+    integration methods.
+
+    Parameters
+    ----------
+    func:
+        The function to integrate. The integration variable should be the last
+        argument of this function.
+    args:
+        The iterable of static arguments to the function.
+    bound1:
+        The lower integration bound
+    bound2:
+        The upper integration bound
+
+    Returns
+    -------
+    The value of the integral of the function between the bounds
+
+    Raises
+    ------
+    MagnetostaticsIntegrationError
+        Integration failed
+    """
+    # warnings.filterwarnings("error", category=IntegrationWarning)
+    try:
+        result = quad(func, bound1, bound2, args=args)[0]
+    except IntegrationWarning:
+        # First attempt at fixing the integration problem
+        points = [
+            0.25 * (bound2 - bound1),
+            0.5 * (bound2 - bound1),
+            0.75 * (bound2 - bound1),
+        ]
+        try:
+            result = quad(func, bound1, bound2, args=args, points=points, limit=200)[0]
+        except IntegrationWarning:
+            raise IntegrationWarning
+
+    # warnings.filterwarnings("default", category=IntegrationWarning)
+    return result
+
+
+def n_integrate(func, args, bounds):
+    """
+    Utility for n-dimensional integration of a function between bounds. Easier
+    to refactor integration methods.
+
+    Parameters
+    ----------
+    func:
+        The function to integrate. The integration variable should be the last
+        argument of this function.
+    args:
+        The iterable of static arguments to the function.
+    bounds:
+        The list of lower and upper integration bounds applied to x[0], x[1], ..
+
+    Returns
+    -------
+    The value of the integral of the function between the bounds
+    """
+    return nquad(func, bounds, args=args)[0]
+
+
+def jit_llc5(f_integrand):
+    """
+    Decorator for 4-argument integrand function to a low-level callable.
+
+    Parameters
+    ----------
+    f_integrand:
+        The integrand function
+
+    Returns
+    -------
+    The decorated integrand function as a LowLevelCallable
+    """
+    f_jitted = numba.jit(f_integrand, nopython=True, cache=True)
+
+    @numba.cfunc(float64(intc, CPointer(float64)))
+    def wrapped(n, xx):  # noqa: ARG001
+        return f_jitted(xx[0], xx[1], xx[2], xx[3], xx[4])
+
+    return LowLevelCallable(wrapped.ctypes)
+
+
+def jit_llc3(f_integrand):
+    """
+    Decorator for 2-argument integrand function to a low-level callable.
+
+    Parameters
+    ----------
+    f_integrand:
+        The integrand function
+
+    Returns
+    -------
+    The decorated integrand function as a LowLevelCallable
+    """
+    f_jitted = numba.jit(f_integrand, nopython=True, cache=True)
+
+    @numba.cfunc(float64(intc, CPointer(float64)))
+    def wrapped(n, xx):  # noqa: ARG001
+        return f_jitted(xx[0], xx[1], xx[2])
+
+    return LowLevelCallable(wrapped.ctypes)
+
+
+def jit_llc7(f_integrand):
+    """
+    Decorator for 6-argument integrand function to a low-level callable.
+
+    Parameters
+    ----------
+    f_integrand:
+        The integrand function
+
+    Returns
+    -------
+    The decorated integrand function as a LowLevelCallable
+    """
+    f_jitted = numba.jit(f_integrand, nopython=True, cache=True)
+
+    @numba.cfunc(float64(intc, CPointer(float64)))
+    def wrapped(n, xx):  # noqa: ARG001
+        return f_jitted(xx[0], xx[1], xx[2], xx[3], xx[4], xx[5], xx[6])
+
+    return LowLevelCallable(wrapped.ctypes)
+
+
+EPS = np.finfo(float).eps
+
+
+@numba.jit(nopython=True, cache=True)
+def _partial_x_integrand(phi, rr, zz):
+    """
+    Integrand edge cases derived to constant integrals. Much faster than
+    splitting up the integrands.
+    """
+    cos_phi = np.cos(phi)
+    r0 = np.sqrt(rr**2 + 1 - 2 * rr * cos_phi + zz**2)
+
+    if abs(zz) < EPS:
+        if abs(rr - 1.0) < EPS:
+            return -1.042258937608 / np.pi
+        if rr < 1.0:
+            return cos_phi * (
+                r0 + cos_phi * np.log((r0 + 1 + rr) / (r0 + 1 - rr))
+            ) - 0.25 * (1 + np.log(4))
+
+    return (r0 + cos_phi * np.log(r0 + rr - cos_phi)) * cos_phi
+
+
+@jit_llc5
+def _full_x_integrand(phi, r1, r2, z1, z2):
+    """
+    Calculate the P_x primitive integral.
+
+    \t:math:`P_{x}(R, Z) = \\int_{0}^{\\pi}[R_{0}+cos(\\phi)ln(R_{0}+R`
+    \t:math:`-cos(\\phi))]cos(\\phi)d\\phi`
+    """
+    return (
+        _partial_x_integrand(phi, r1, z1)
+        - _partial_x_integrand(phi, r1, z2)
+        - _partial_x_integrand(phi, r2, z1)
+        + _partial_x_integrand(phi, r2, z2)
+    )
+
+
+def _partial_z_integrand_nojit(phi, rr, zz):
+    """
+    Integrand edge cases derived to constant integrals. Much faster than
+    splitting up the integrands.
+    """
+    if abs(zz) < EPS:
+        return 0.0
+
+    sin_phi = np.sin(phi)
+    cos_phi = np.cos(phi)
+    r0 = np.sqrt(rr**2 + 1 - 2 * rr * cos_phi + zz**2)
+
+    # F1
+    result = zz * np.log(r0 + rr - cos_phi) - cos_phi * np.log(r0 + zz)
+
+    # F2
+    result = result - 0.5 * rr if rr - 1 < EPS else result - 0.5 / rr
+    # F3
+    if 0.5 * np.pi * sin_phi > 1e-9:  # noqa: PLR2004
+        result -= sin_phi * np.arctan(zz * (rr - cos_phi) / (r0 * sin_phi))
+    return result
+
+
+_partial_z_integrand = numba.jit(_partial_z_integrand_nojit, nopython=True, cache=True)
+_partial_z_integrand_llc = jit_llc3(_partial_z_integrand_nojit)
+
+
+def _integrate_z_by_parts(r1, r2, z1, z2):
+    """
+    Integrate the Bz integrand by parts.
+
+    This can be used as a fallback if the full integration fails.
+    """
+    return (
+        integrate(_partial_z_integrand_llc, (r1, z1), 0, np.pi)
+        - integrate(_partial_z_integrand_llc, (r1, z2), 0, np.pi)
+        - integrate(_partial_z_integrand_llc, (r2, z1), 0, np.pi)
+        + integrate(_partial_z_integrand_llc, (r2, z2), 0, np.pi)
+    )
+
+
+@jit_llc7
+def _full_psi_integrand(rp, phi, rc, zc, zp, d_rc, d_zc):
+    """
+    Integrand for psi = xBz
+    """
+    z = zp - zc  # numba issue # noqa: PLR6104
+    r1, r2 = (rc - d_rc) / rp, (rc + d_rc) / rp
+    z1, z2 = (-d_zc - z) / rp, (d_zc - z) / rp
+    return rp**2 * (
+        _partial_z_integrand(phi, r1, z1)
+        - _partial_z_integrand(phi, r1, z2)
+        - _partial_z_integrand(phi, r2, z1)
+        + _partial_z_integrand(phi, r2, z2)
+    )
+
+
+@numba.jit(nopython=True, cache=True)
+def _get_working_coords(rc, zc, rp, zp, d_rc, d_zc):
+    """
+    Convert coil and global coordinates to working coordinates.
+    """
+    z = zp - zc
+    r1, r2 = (rc - d_rc) / rp, (rc + d_rc) / rp
+    z1, z2 = (-d_zc - z) / rp, (d_zc - z) / rp
+    j_tor = 1 / (4 * d_rc * d_zc)  # Keep current out of the equation
+    return r1, r2, z1, z2, j_tor
+
+
+def is_num(thing):
+    """
+    Determine whether or not the input is a number.
+
+    Parameters
+    ----------
+    thing:
+        The input which we need to determine is a number or not
+
+    Returns
+    -------
+    Whether or not the input is a number
+    """
+    if thing is True or thing is False:
+        return False
+    try:
+        thing = floatify(thing)
+    except (ValueError, TypeError):
+        return False
+    else:
+        return not np.isnan(thing)
+
+
+def _array_dispatcher(func):
+    """
+    Decorator for float and array handling.
+    """
+
+    def wrapper(rc, zc, rp, zp, d_rc, d_zc):
+        # Handle floats
+        if is_num(rp):
+            return func(rc, zc, rp, zp, d_rc, d_zc)
+
+        # Handle arrays
+        if len(rp.shape) == 1:
+            if not isinstance(rc, np.ndarray) or len(rc.shape) == 1:
+                result = np.zeros(len(rp))
+                for i in range(len(rp)):
+                    result[i] = floatify(func(rc, zc, rp[i], zp[i], d_rc, d_zc))
+            else:
+                result = np.zeros((len(rp), len(rc)))
+                for j in range(rc.shape[1]):
+                    for i in range(len(rp)):
+                        result[i, j] = floatify(
+                            func(
+                                rc[:, j], zc[:, j], rp[i], zp[i], d_rc[:, j], d_zc[:, j]
+                            )
+                        )
+        # 2-D arrays
+        elif not isinstance(rc, np.ndarray) or len(rc.shape) == 1:
+            result = np.zeros(rp.shape)
+            for i in range(rp.shape[0]):
+                for j in range(zp.shape[1]):
+                    result[i, j] = floatify(
+                        func(rc, zc, rp[i, j], zp[i, j], d_rc, d_zc)
+                    )
+        else:
+            result = np.zeros([*list(rp.shape), rc.shape[1]])
+            for k in range(rc.shape[1]):
+                for i in range(rp.shape[0]):
+                    for j in range(zp.shape[1]):
+                        result[i, j, ..., k] = func(
+                            rc[:, k],
+                            zc[:, k],
+                            rp[i, j],
+                            zp[i, j],
+                            d_rc[:, k],
+                            d_zc[:, k],
+                        )
+        return result
+
+    return wrapper
+
+
+@_array_dispatcher
+def semianalytic(rc, zc, rp, zp, d_rc, d_zc):
+    """
+    Function that uses a semi-analytic reduction of the Biot-Savart law to
+    calculate the magnetic fields (Bx, Bz) and magnetic flux (psi) from a
+    rectangular cross-section circular current. Output is from unit current
+    so for actual field multiply result by coil current.
+
+    Used for producing the radial (x) and vertical (z) components of the magnetic
+    field, alongside the magnetic flux at any point inside the solenoid (PF or CS)
+    that contains the source point.
+
+    rc: input real array : coil r coordinates (m)
+    zc: input real array : coil z coordinates (m)
+    rp: input real array : calculation r coordinates (m)
+    zp: input real array : calculation z coordinates (m)
+    d_rc: input real array : half-width of the coil (m)
+    d_zc: input real array : half-height of the coil (m)
+    br: output real : radial field component at (rp,zp) (T)
+    bz : output real : vertical field component at (rp,zp) (T)
+    psi : output real : psi component at (rp, zp) (Wb)
+    """
+    r1, r2, z1, z2, j_tor = _get_working_coords(rc, zc, rp, zp, d_rc, d_zc)
+    Br = integrate(
+        _full_x_integrand, tuple(np.asarray((r1, r2, z1, z2)).ravel()), 0, np.pi
+    )
+
+    Bz = _integrate_z_by_parts(r1, r2, z1, z2)
+
+    psi = n_integrate(
+        _full_psi_integrand,
+        (floatify(rc), floatify(zc), floatify(zp), floatify(d_rc), floatify(d_zc)),
+        [[0, floatify(rp)], [0, np.pi]],
+    )
+    fac_B = 2e-7 * j_tor * rp  # MU_0/(2*np.pi)
+    fac_psi = 2e-7 * j_tor
+
+    return fac_B * Br, fac_B * Bz, fac_psi * psi
+
+
 @numba.jit(nopython=True)
 def elliptic_inductance(radius, rc):
     """

From 9ca0ec9c6dc864249609af0ae332847596dcb11f Mon Sep 17 00:00:00 2001
From: Jonathan Matthews <jonathan.matthews@ukaea.uk>
Date: Wed, 16 Apr 2025 09:20:39 +0100
Subject: [PATCH 3/7] replace filament approach for CS coil and use different b
 field calculation method

---
 process/pfcoil.py | 413 ++++++++++++++++++----------------------------
 1 file changed, 165 insertions(+), 248 deletions(-)

diff --git a/process/pfcoil.py b/process/pfcoil.py
index e1ef1af647..aa0ea286d2 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -679,12 +679,9 @@ def pfcoil(self):
             iii = ii
             for ij in range(pfv.n_pf_coils_in_group[ii]):
                 # Peak field
-
                 if ij == 0:
                     # Index args +1ed
-                    bri, bro, bzi, bzo = self.peakb(
-                        i + 1, iii + 1, it
-                    )  # returns b_pf_coil_peak, bpf2
+                    bri, bro, bzi, bzo = self.peakb(i, iii, it)  # returns b_pf_coil_peak, bpf2
 
                 # Issue 1871.  MDK
                 # Allowable current density (for superconducting coils) for each coil, index i
@@ -1174,42 +1171,70 @@ def ohcalc(self):
         # Peak field at the End-Of-Flattop (EOF)
         # Occurs at inner edge of coil; bmaxoh2 and bzi are of opposite sign at EOF
 
-        # Peak field due to central Solenoid itself
-        bmaxoh2 = self.bfmax(
-            pfv.j_cs_flat_top_end,
-            pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
-            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1],
-            hohc,
+        current = pfv.coheof * pfv.areaoh
+        # Self field from CS
+        br_inner, bz_inner, psi_inner = semianalytic(
+            pfv.rpf[pfv.nohc - 1],
+            pfv.zpf[pfv.nohc - 1],
+            pfv.ra[pfv.nohc - 1],
+            pfv.zpf[pfv.nohc - 1],
+            pfv.rb[pfv.nohc - 1] - pfv.ra[pfv.nohc - 1],
+            pfv.zh[pfv.nohc - 1] - pfv.zl[pfv.nohc - 1],
+            current,
+        )
+        br_outer, bz_outer, psi_outer = semianalytic(
+            pfv.rpf[pfv.nohc - 1],
+            pfv.zpf[pfv.nohc - 1],
+            pfv.rb[pfv.nohc - 1],
+            pfv.zpf[pfv.nohc - 1],
+            pfv.rb[pfv.nohc - 1] - pfv.ra[pfv.nohc - 1],
+            pfv.zh[pfv.nohc - 1] - pfv.zl[pfv.nohc - 1],
+            current,
         )
 
         # Peak field due to other PF coils plus plasma
         timepoint = 5
-        bri, bro, bzi, bzo = self.peakb(pfv.n_cs_pf_coils, 99, timepoint)
+        bri, bro, bzi, bzo = self.peakb(pfv.n_cs_pf_coils - 1, 99 - 1, timepoint)
 
-        pfv.b_cs_peak_flat_top_end = abs(bzi - bmaxoh2)
+        pfv.b_cs_peak_flat_top_end = abs(bzi + bz_inner)
 
         # Peak field on outboard side of central Solenoid
         # (self-field is assumed to be zero - long solenoid approximation)
-        bohco = abs(bzo)
+        bohco = abs(bzo + bz_outer)
 
         # Peak field at the Beginning-Of-Pulse (BOP)
-        # Occurs at inner edge of coil; b_cs_peak_pulse_start and bzi are of same sign at BOP
-        pfv.b_cs_peak_pulse_start = self.bfmax(
-            pfv.j_cs_pulse_start,
-            pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
-            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1],
-            hohc,
+        # Occurs at inner edge of coil; bmaxoh0 and bzi are of same sign at BOP
+
+        current = pfv.cohbop * pfv.areaoh
+        # Self field
+        br_inner, bz_inner, psi_inner = semianalytic(
+            pfv.rpf[pfv.nohc - 1],
+            pfv.zpf[pfv.nohc - 1],
+            pfv.ra[pfv.nohc - 1],
+            pfv.zpf[pfv.nohc - 1],
+            pfv.rb[pfv.nohc - 1] - pfv.ra[pfv.nohc - 1],
+            pfv.zh[pfv.nohc - 1] - pfv.zl[pfv.nohc - 1],
+            current,
+        )
+        br_outer, bz_outer, psi_outer = semianalytic(
+            pfv.rpf[pfv.nohc - 1],
+            pfv.zpf[pfv.nohc - 1],
+            pfv.rb[pfv.nohc - 1],
+            pfv.zpf[pfv.nohc - 1],
+            pfv.rb[pfv.nohc - 1] - pfv.ra[pfv.nohc - 1],
+            pfv.zh[pfv.nohc - 1] - pfv.zl[pfv.nohc - 1],
+            current,
         )
         timepoint = 2
-        bri, bro, bzi, bzo = self.peakb(pfv.n_cs_pf_coils, 99, timepoint)
+        bri, bro, bzi, bzo = self.peakb(pfv.n_cs_pf_coils - 1, 99 - 1, timepoint)
 
-        pfv.b_cs_peak_pulse_start = abs(pfv.b_cs_peak_pulse_start + bzi)
+        pfv.b_cs_peak_pulse_start = abs(bz_inner + bzi)
 
         # Maximum field values
         pfv.b_pf_coil_peak[pfv.n_cs_pf_coils - 1] = max(
             pfv.b_cs_peak_flat_top_end, abs(pfv.b_cs_peak_pulse_start)
         )
-        pf.bpf2[pfv.n_cs_pf_coils - 1] = max(bohco, abs(bzo))
+        pf.bpf2[pfv.n_cs_pf_coils - 1] = max(bohco, abs(bzo + bz_outer))
 
         # Stress ==> cross-sectional area of supporting steel to use
         if pfv.i_pf_conductor == 0:
@@ -1395,33 +1420,35 @@ def peakb(self, i, ii, it):
         The calculation includes the effects from all the coils
         and the plasma.
         """
-        if bv.iohcl != 0 and i == pfv.n_cs_pf_coils:
+        if (bv.iohcl != 0) and (i == pfv.n_cs_pf_coils - 1):
             # Peak field is to be calculated at the Central Solenoid itself,
             # so exclude its own contribution; its self field is
-            # dealt with externally using routine BFMAX
+            # dealt with inside ohcalc
             kk = 0
+            rp = np.array([pfv.ra[i], pfv.rb[i]])
+            zp = np.array([pfv.zpf[i], pfv.zpf[i]])
         else:
             # Check different times for maximum current
             if (
                 abs(
-                    pfv.c_pf_cs_coil_pulse_start_ma[i - 1]
-                    - pfv.c_pf_cs_coils_peak_ma[i - 1]
+                    pfv.c_pf_cs_coil_pulse_start_ma[i]
+                    - pfv.c_pf_cs_coils_peak_ma[i]
                 )
                 < 1.0e-12
             ):
                 it = 2
             elif (
                 abs(
-                    pfv.c_pf_cs_coil_flat_top_ma[i - 1]
-                    - pfv.c_pf_cs_coils_peak_ma[i - 1]
+                    pfv.c_pf_cs_coil_flat_top_ma[i]
+                    - pfv.c_pf_cs_coils_peak_ma[i]
                 )
                 < 1.0e-12
             ):
                 it = 4
             elif (
                 abs(
-                    pfv.c_pf_cs_coil_pulse_end_ma[i - 1]
-                    - pfv.c_pf_cs_coils_peak_ma[i - 1]
+                    pfv.c_pf_cs_coil_pulse_end_ma[i]
+                    - pfv.c_pf_cs_coils_peak_ma[i]
                 )
                 < 1.0e-12
             ):
@@ -1452,70 +1479,51 @@ def peakb(self, i, ii, it):
 
                 kk = pf.nfxf
 
-        # Non-Central Solenoid coils' contributions
+        # initialise variables for loop
+        current_array = np.array([])
+        rc_array = np.array([])
+        zc_array = np.array([])
         jj = 0
+        # loops to set up arrays for magnetostatic calculations
         for iii in range(pfv.n_pf_coil_groups):
             for _jjj in range(pfv.n_pf_coils_in_group[iii]):
-                jj = jj + 1
-                # Radius, z-coordinate and current for each coil
-                if iii == ii - 1:
-                    # Self field from coil (Lyle's Method)
-                    kk = kk + 1
-
-                    dzpf = pfv.z_pf_coil_upper[jj - 1] - pfv.z_pf_coil_lower[jj - 1]
-                    pf.rfxf[kk - 1] = pfv.r_pf_coil_middle[jj - 1]
-                    pf.zfxf[kk - 1] = pfv.z_pf_coil_middle[jj - 1] + dzpf * 0.125e0
-                    pf.cfxf[kk - 1] = (
-                        pfv.c_pf_cs_coils_peak_ma[jj - 1]
-                        * pfv.waves[jj - 1, it - 1]
-                        * 0.25e6
-                    )
-                    kk = kk + 1
-                    pf.rfxf[kk - 1] = pfv.r_pf_coil_middle[jj - 1]
-                    pf.zfxf[kk - 1] = pfv.z_pf_coil_middle[jj - 1] + dzpf * 0.375e0
-                    pf.cfxf[kk - 1] = (
-                        pfv.c_pf_cs_coils_peak_ma[jj - 1]
-                        * pfv.waves[jj - 1, it - 1]
-                        * 0.25e6
-                    )
-                    kk = kk + 1
-                    pf.rfxf[kk - 1] = pfv.r_pf_coil_middle[jj - 1]
-                    pf.zfxf[kk - 1] = pfv.z_pf_coil_middle[jj - 1] - dzpf * 0.125e0
-                    pf.cfxf[kk - 1] = (
-                        pfv.c_pf_cs_coils_peak_ma[jj - 1]
-                        * pfv.waves[jj - 1, it - 1]
-                        * 0.25e6
-                    )
-                    kk = kk + 1
-                    pf.rfxf[kk - 1] = pfv.r_pf_coil_middle[jj - 1]
-                    pf.zfxf[kk - 1] = pfv.z_pf_coil_middle[jj - 1] - dzpf * 0.375e0
-                    pf.cfxf[kk - 1] = (
-                        pfv.c_pf_cs_coils_peak_ma[jj - 1]
-                        * pfv.waves[jj - 1, it - 1]
-                        * 0.25e6
-                    )
-
+                current = pfv.c_pf_cs_coils_peak_ma[jj] * pfv.waves[jj, it - 1] * 1.0e6
+                if jj == i:
+                    rp = np.array([pfv.r_pf_coil_inner[i], pfv.r_pf_coil_outer[i]])
+                    zp = np.array([pfv.z_pf_coil_middle[jj], pfv.z_pf_coil_middle[jj]])
+                    dr = pfv.r_pf_coil_outer[i] - pfv.r_pf_coil_inner[i]
+                    dz = pfv.z_pf_coil_upper[jj] - pfv.z_pf_coil_lower[jj]
+                    self_current = current
+                    self_rc = pfv.r_pf_coil_middle[jj]
+                    self_zc = pfv.z_pf_coil_middle[jj]
                 else:
-                    # Field from different coil
-                    kk = kk + 1
-                    pf.rfxf[kk - 1] = pfv.r_pf_coil_middle[jj - 1]
-                    pf.zfxf[kk - 1] = pfv.z_pf_coil_middle[jj - 1]
-                    pf.cfxf[kk - 1] = (
-                        pfv.c_pf_cs_coils_peak_ma[jj - 1]
-                        * pfv.waves[jj - 1, it - 1]
-                        * 1.0e6
-                    )
+                    current_array = np.append(current_array, current)
+                    rc_array = np.append(rc_array, pfv.r_pf_coil_middle[jj])
+                    zc_array = np.append(zc_array, pfv.z_pf_coil_middle[jj])
+                jj += 1
 
-        # Plasma contribution
+        # plasma contribution
         if it > 2:
-            kk = kk + 1
-            pf.rfxf[kk - 1] = pv.rmajor
-            pf.zfxf[kk - 1] = 0.0e0
-            pf.cfxf[kk - 1] = pv.plasma_current
-
-        # Calculate the field at the inner and outer edges
-        # of the coil of interest
-        pf.xind[:kk], bri, bzi, psi = bfield(
+            current_array = np.append(current_array, pv.plasma_current)
+            rc_array = np.append(rc_array, pv.rmajor)
+            zc_array = np.append(zc_array, 0.0e0)
+        bri, bzi, _ = bfield(rc_array, zc_array, current_array, rp[0], zp[0])
+        bro, bzo, _ = bfield(rc_array, zc_array, current_array, rp[1], zp[1])
+        if not ((bv.iohcl != 0) and (i == pfv.nohc - 1)):
+            bri_self, bzi_self, _ = semianalytic(
+                self_rc, self_zc, rp[0], zp[0], dr, dz, self_current
+            )
+            bro_self, bzo_self, _ = semianalytic(
+                self_rc, self_zc, rp[1], zp[1], dr, dz, self_current
+            )
+
+            bri += bri_self
+            bzi += bzi_self
+            bro += bro_self
+            bzo += bzo_self
+
+        # remove xind as not used?
+        pf.xind[:kk] = mutual_inductance(
             pf.rfxf[:kk],
             pf.zfxf[:kk],
             pf.cfxf[:kk],
@@ -1526,85 +1534,22 @@ def peakb(self, i, ii, it):
             pf.rfxf[:kk],
             pf.zfxf[:kk],
             pf.cfxf[:kk],
-            pfv.r_pf_coil_outer[i - 1],
-            pfv.z_pf_coil_middle[i - 1],
+            pfv.rb[i - 1],
+            pfv.zpf[i - 1],
         )
 
         # b_pf_coil_peak and bpf2 for the Central Solenoid are calculated in OHCALC
-        if (bv.iohcl != 0) and (i == pfv.n_cs_pf_coils):
+        if (bv.iohcl != 0) and (i == pfv.n_cs_pf_coils - 1):
             return bri, bro, bzi, bzo
 
         bpfin = math.sqrt(bri**2 + bzi**2)
         bpfout = math.sqrt(bro**2 + bzo**2)
-        for n in range(pfv.n_pf_coils_in_group[ii - 1]):
-            pfv.b_pf_coil_peak[i - 1 + n] = bpfin
-            pf.bpf2[i - 1 + n] = bpfout
+        for n in range(pfv.n_pf_coils_in_group[ii]):
+            pfv.b_pf_coil_peak[i + n] = bpfin
+            pf.bpf2[i + n] = bpfout
 
         return bri, bro, bzi, bzo
 
-    def bfmax(self, rj, a, b, h):
-        """Calculates the maximum field of a solenoid.
-
-        author: P J Knight, CCFE, Culham Science Centre
-        This routine calculates the peak field (T) at a solenoid's
-        inner radius, using fits taken from the figure
-        on p.22 of M. Wilson's book Superconducting Magnets,
-        Clarendon Press, Oxford, N.Y., 1983
-
-        :param rj: overall current density (A/m2)
-        :type rj: float
-        :param a: solenoid inner radius (m)
-        :type a: float
-        :param b: solenoid outer radius (m)
-        :type b: float
-        :param h: solenoid half height (m)
-        :type h: float
-        :return bfmax: maximum field of solenoid
-        :rtype: float
-        """
-        beta = h / a
-        alpha = b / a
-
-        # Fits are for 1 < alpha < 2 , and 0.5 < beta < very large
-        b0 = (
-            rj
-            * constants.rmu0
-            * h
-            * math.log(
-                (alpha + math.sqrt(alpha**2 + beta**2))
-                / (1.0 + math.sqrt(1.0 + beta**2))
-            )
-        )
-
-        if beta > 3.0:
-            b1 = constants.rmu0 * rj * (b - a)
-            f = (3.0 / beta) ** 2
-            bfmax = f * b0 * (1.007 + (alpha - 1.0) * 0.0055) + (1.0 - f) * b1
-
-        elif beta > 2.0:
-            rat = (1.025 - (beta - 2.0) * 0.018) + (alpha - 1.0) * (
-                0.01 - (beta - 2.0) * 0.0045
-            )
-            bfmax = rat * b0
-
-        elif beta > 1.0:
-            rat = (1.117 - (beta - 1.0) * 0.092) + (alpha - 1.0) * (beta - 1.0) * 0.01
-            bfmax = rat * b0
-
-        elif beta > 0.75:
-            rat = (1.30 - 0.732 * (beta - 0.75)) + (alpha - 1.0) * (
-                0.2 * (beta - 0.75) - 0.05
-            )
-            bfmax = rat * b0
-
-        else:
-            rat = (1.65 - 1.4 * (beta - 0.5)) + (alpha - 1.0) * (
-                0.6 * (beta - 0.5) - 0.20
-            )
-            bfmax = rat * b0
-
-        return bfmax
-
     def vsec(self):
         """Calculation of volt-second capability of PF system.
 
@@ -1823,13 +1768,9 @@ def induct(self, output):
         nohmax = 200
         nplas = 1
 
-        br = 0.0
-        bz = 0.0
-        psi = 0.0
         rc = np.zeros(pfv.ngc2 + nohmax)
         zc = np.zeros(pfv.ngc2 + nohmax)
         xc = np.zeros(pfv.ngc2 + nohmax)
-        cc = np.zeros(pfv.ngc2 + nohmax)
         xcin = np.zeros(pfv.ngc2 + nohmax)
         xcout = np.zeros(pfv.ngc2 + nohmax)
         rplasma = np.zeros(nplas)
@@ -1916,8 +1857,8 @@ def induct(self, output):
 
                 reqv = rp * (1.0e0 + delzoh**2 / (24.0e0 * rp**2))
 
-                xcin, br, bz, psi = bfield(rc, zc, cc, reqv - deltar, zp)
-                xcout, br, bz, psi = bfield(rc, zc, cc, reqv + deltar, zp)
+                xcin = mutual_inductance(rc, zc, reqv - deltar, zp)
+                xcout = mutual_inductance(rc, zc, reqv + deltar, zp)
 
                 for ii in range(nplas):
                     xc[ii] = 0.5e0 * (xcin[ii] + xcout[ii])
@@ -1945,7 +1886,7 @@ def induct(self, output):
             ncoils = ncoils + pfv.n_pf_coils_in_group[i]
             rp = pfv.r_pf_coil_middle[ncoils - 1]
             zp = pfv.z_pf_coil_middle[ncoils - 1]
-            xc, br, bz, psi = bfield(rc, zc, cc, rp, zp)
+            xc = mutual_inductance(rc, zc, rp, zp)
             for ii in range(nplas):
                 xpfpl = xpfpl + xc[ii]
 
@@ -1983,7 +1924,7 @@ def induct(self, output):
                 ncoils = ncoils + pfv.n_pf_coils_in_group[i]
                 rp = pfv.r_pf_coil_middle[ncoils - 1]
                 zp = pfv.z_pf_coil_middle[ncoils - 1]
-                xc, br, bz, psi = bfield(rc, zc, cc, rp, zp)
+                xc = mutual_inductance(rc, zc, rp, zp)
                 for ii in range(noh):
                     xohpf = xohpf + xc[ii]
 
@@ -2014,7 +1955,7 @@ def induct(self, output):
 
             rp = pfv.r_pf_coil_middle[i]
             zp = pfv.z_pf_coil_middle[i]
-            xc, br, bz, psi = bfield(rc, zc, cc, rp, zp)
+            xc = mutual_inductance(rc, zc, rp, zp)
             for k in range(pf.nef):
                 if k < i:
                     pfv.ind_pf_cs_plasma_mutual[i, k] = (
@@ -3182,28 +3123,17 @@ def superconpf(
         return jcritwp, j_crit_cable, j_crit_sc, tmarg
 
 
-@numba.njit(cache=True)
-def bfield(rc, zc, cc, rp, zp):
-    """Calculate the field at a point due to currents in a number
-    of circular poloidal conductor loops.
-    author: P J Knight, CCFE, Culham Science Centre
-    author: D Strickler, ORNL
-    author: J Galambos, ORNL
-    nc : input integer : number of loops
+def mutual_inductance(rc, zc, rp, zp):
+    """
+    Calcualtes the mutual inductances between the loops and a poloidal
+    filament at the (R,Z) point of interest.
+
+    nc : number of loops
     rc(nc) : input real array : R coordinates of loops (m)
     zc(nc) : input real array : Z coordinates of loops (m)
-    cc(nc) : input real array : Currents in loops (A)
     xc(nc) : output real array : Mutual inductances (H)
     rp, zp : input real : coordinates of point of interest (m)
-    br : output real : radial field component at (rp,zp) (T)
-    bz : output real : vertical field component at (rp,zp) (T)
-    psi : output real : poloidal flux at (rp,zp) (Wb)
-    This routine calculates the magnetic field components and
-    the poloidal flux at an (R,Z) point, given the locations
-    and currents of a set of conductor loops.
-    <P>The mutual inductances between the loops and a poloidal
-    filament at the (R,Z) point of interest is also found."""
-
+    """
     #  Elliptic integral coefficients
 
     a0 = 1.38629436112
@@ -3228,9 +3158,6 @@ def bfield(rc, zc, cc, rp, zp):
     nc = len(rc)
 
     xc = np.empty((nc,))
-    br = 0
-    bz = 0
-    psi = 0
 
     for i in range(nc):
         d = (rp + rc[i]) ** 2 + (zp - zc[i]) ** 2
@@ -3269,21 +3196,33 @@ def bfield(rc, zc, cc, rp, zp):
 
         xc[i] = 0.5 * RMU0 * sd * ((2.0 - s) * xk - 2.0 * xe)
 
-        #  Radial, vertical fields
-
-        # brx = (
-        #     RMU0
-        #     * cc[i]
-        #     * dz
-        #     / (2 * np.pi * rp * sd)
-        #     * (-xk + (rc[i] ** 2 + rp**2 + zs) / (dr**2 + zs) * xe)
-        # )
-        # bzx = (
-        #     RMU0
-        #     * cc[i]
-        #     / (2 * np.pi * sd)
-        #     * (xk + (rc[i] ** 2 - rp**2 - zs) / (dr**2 + zs) * xe)
-        # )
+    return xc
+
+
+@numba.njit(cache=True)
+def bfield(rc, zc, cc, rp, zp):
+    """Calculate the field at a point due to currents in a number
+    of circular poloidal conductor loops.
+    author: P J Knight, CCFE, Culham Science Centre
+    author: D Strickler, ORNL
+    author: J Galambos, ORNL
+    nc : input integer : number of loops
+    rc(nc) : input real array : R coordinates of loops (m)
+    zc(nc) : input real array : Z coordinates of loops (m)
+    cc(nc) : input real array : Currents in loops (A)
+    rp, zp : input real : coordinates of point of interest (m)
+    br : output real : radial field component at (rp,zp) (T)
+    bz : output real : vertical field component at (rp,zp) (T)
+    psi : output real : poloidal flux at (rp,zp) (Wb)
+    This routine calculates the magnetic field components and
+    the poloidal flux at an (R,Z) point, given the locations
+    and currents of a set of conductor loops."""
+
+    br = 0
+    bz = 0
+    psi = 0
+
+    for i in range(len(rc)):
         # get the field due to unit current using Green's function
         g_br, g_bz, g_psi = greens(rc[i], zc[i], rp, zp)
         # convert to actual field from unit current
@@ -3294,7 +3233,7 @@ def bfield(rc, zc, cc, rp, zp):
         bz += g_bz * cc[i]
         psi += g_psi * cc[i]
 
-    return xc, br, bz, psi
+    return br, bz, psi
 
 
 @numba.njit(cache=True)
@@ -3404,7 +3343,7 @@ def fixb(lrow1, npts, rpts, zpts, nfix, rfix, zfix, cfix):
     for i in range(npts):
         # bfield() only operates correctly on nfix slices of array
         # arguments, not entire arrays
-        _, brw, bzw, _ = bfield(rfix[:nfix], zfix[:nfix], cfix[:nfix], rpts[i], zpts[i])
+        brw, bzw, _ = bfield(rfix[:nfix], zfix[:nfix], cfix[:nfix], rpts[i], zpts[i])
         bfix[i] = brw
         bfix[npts + i] = bzw
 
@@ -3484,7 +3423,7 @@ def mtrx(
         for j in range(n_pf_coil_groups):
             nc = n_pf_coils_in_group[j]
 
-            _, gmat[i, j], gmat[i + npts, j], _ = bfield(
+            gmat[i, j], gmat[i + npts, j], _ = bfield(
                 rcls[j, :nc], zcls[j, :nc], cc[:nc], rpts[i], zpts[i]
             )
 
@@ -3737,9 +3676,9 @@ def greens(rc, zc, rp, zp):
     zc: input real array : coil z coordinates (m)
     rp: input real array : calculation r coordinates (m)
     zp: input real array : calculation z coordinates (m)
-    g_psi : output real : poloidal flux at (rp,zp) (Wb)
     g_br: output real : radial field component at (rp,zp) (T)
     g_bz : output real : vertical field component at (rp,zp) (T)
+    g_psi : output real : poloidal flux at (rp,zp) (Wb)
     """
 
     a, k2 = calc_a_k2(rc, zc, rp, zp)
@@ -3757,7 +3696,7 @@ def greens(rc, zc, rp, zp):
     )
     g_psi = 0.25 * RMU0 / np.pi * a * ((2 - k2) * k - 2 * e)
 
-    return g_psi, g_br, g_bz
+    return g_br, g_bz, g_psi
 
 
 def floatify(x):
@@ -4051,24 +3990,32 @@ def _array_dispatcher(func):
     Decorator for float and array handling.
     """
 
-    def wrapper(rc, zc, rp, zp, d_rc, d_zc):
+    def wrapper(rc, zc, rp, zp, d_rc, d_zc, current):
         # Handle floats
         if is_num(rp):
-            return func(rc, zc, rp, zp, d_rc, d_zc)
+            return func(rc, zc, rp, zp, d_rc, d_zc, current)
 
         # Handle arrays
         if len(rp.shape) == 1:
             if not isinstance(rc, np.ndarray) or len(rc.shape) == 1:
                 result = np.zeros(len(rp))
                 for i in range(len(rp)):
-                    result[i] = floatify(func(rc, zc, rp[i], zp[i], d_rc, d_zc))
+                    result[i] = floatify(
+                        func(rc, zc, rp[i], zp[i], d_rc, d_zc, current)
+                    )
             else:
                 result = np.zeros((len(rp), len(rc)))
                 for j in range(rc.shape[1]):
                     for i in range(len(rp)):
                         result[i, j] = floatify(
                             func(
-                                rc[:, j], zc[:, j], rp[i], zp[i], d_rc[:, j], d_zc[:, j]
+                                rc[:, j],
+                                zc[:, j],
+                                rp[i],
+                                zp[i],
+                                d_rc[:, j],
+                                d_zc[:, j],
+                                current[:, j],
                             )
                         )
         # 2-D arrays
@@ -4077,7 +4024,7 @@ def wrapper(rc, zc, rp, zp, d_rc, d_zc):
             for i in range(rp.shape[0]):
                 for j in range(zp.shape[1]):
                     result[i, j] = floatify(
-                        func(rc, zc, rp[i, j], zp[i, j], d_rc, d_zc)
+                        func(rc, zc, rp[i, j], zp[i, j], d_rc, d_zc, current)
                     )
         else:
             result = np.zeros([*list(rp.shape), rc.shape[1]])
@@ -4091,6 +4038,7 @@ def wrapper(rc, zc, rp, zp, d_rc, d_zc):
                             zp[i, j],
                             d_rc[:, k],
                             d_zc[:, k],
+                            current[:, k],
                         )
         return result
 
@@ -4098,7 +4046,7 @@ def wrapper(rc, zc, rp, zp, d_rc, d_zc):
 
 
 @_array_dispatcher
-def semianalytic(rc, zc, rp, zp, d_rc, d_zc):
+def semianalytic(rc, zc, rp, zp, d_rc, d_zc, current):
     """
     Function that uses a semi-analytic reduction of the Biot-Savart law to
     calculate the magnetic fields (Bx, Bz) and magnetic flux (psi) from a
@@ -4119,6 +4067,9 @@ def semianalytic(rc, zc, rp, zp, d_rc, d_zc):
     bz : output real : vertical field component at (rp,zp) (T)
     psi : output real : psi component at (rp, zp) (Wb)
     """
+    # stops rp going to 0 as this results in a divide by 0 error in "_get_working_coords"
+    if rp < 0.0001:
+        rp = 0.0001
     r1, r2, z1, z2, j_tor = _get_working_coords(rc, zc, rp, zp, d_rc, d_zc)
     Br = integrate(
         _full_x_integrand, tuple(np.asarray((r1, r2, z1, z2)).ravel()), 0, np.pi
@@ -4134,38 +4085,4 @@ def semianalytic(rc, zc, rp, zp, d_rc, d_zc):
     fac_B = 2e-7 * j_tor * rp  # MU_0/(2*np.pi)
     fac_psi = 2e-7 * j_tor
 
-    return fac_B * Br, fac_B * Bz, fac_psi * psi
-
-
-@numba.jit(nopython=True)
-def elliptic_inductance(radius, rc):
-    """
-    Calculate the inductance of a circular coil by elliptic integrals.
-
-    Parameters
-    ----------
-    radius: input real : The radius of the circular coil
-    rc: input real : The radius of the coil cross-section
-    xc: output real : The self-inductance of the circular coil [H]
-
-    Notes
-    -----
-    The inductance is given by
-
-    .. math::
-        L = \\mu_{0} (2 r - r_c) \\Biggl((1 - k^2 / 2)~
-        \\int_0^{\\frac{\\pi}{2}} \\frac{d\\theta}{\\sqrt{1 - k~
-        \\sin (\\theta)^2}} - \\int_0^{\\frac{\\pi}{2}}~
-        \\sqrt{1 - k \\sin (\\theta)^2} \\, d\\theta\\Biggr)
-
-    where :math:`r` is the radius, :math:`\\mu_{0}` is the vacuum
-    permeability, and
-
-    .. math::
-        k = \\max\\left(10^{-8}, \\min~
-        \\left(\\frac{4r(r - r_c)}{(2r - r_c)^2}~
-        , 1.0 - 10^{-8}\\right)\\right)
-    """
-    k = 4 * radius * (radius - rc) / (2 * radius - rc) ** 2
-    k = clip_nb(k, GREENS_ZERO, 1.0 - GREENS_ZERO)
-    return RMU0 * (2 * radius - rc) * ((1 - k**2 / 2) * ellipk_nb(k) - ellipe_nb(k))
+    return current * fac_B * Br, current * fac_B * Bz, current * fac_psi * psi

From 4f453b544910fb9d122fe0ea06439281c1a09fff Mon Sep 17 00:00:00 2001
From: Jonathan Matthews <jonathan.matthews@ukaea.uk>
Date: Wed, 16 Apr 2025 09:21:03 +0100
Subject: [PATCH 4/7] numba-scipy addition to setup.py

---
 setup.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index 4c33269242..691c05acf0 100644
--- a/setup.py
+++ b/setup.py
@@ -36,6 +36,7 @@
     "install_requires": [
         "numpy>=1.23,<2",
         "scipy>=1.10",
+        "numba-scipy @ git+https://github.com/numba/numba-scipy@23c3b33440ea1fe0f84d05d269fb4a3df4b92787",
         "cvxpy!=1.3.0,!=1.3.1",
         "osqp<1.0",
         "pandas>=2.0",
@@ -52,7 +53,10 @@
         "test": ["pytest>=5.4.1", "requests>=2.30", "testbook>=0.4"],
         "examples": ["pillow>=5.1.0", "jupyter==1.0.0", "pdf2image==1.16.0"],
     },
-    "entry_points": {"console_scripts": ["process=process.main:main"]},
+    "entry_points": {
+        "console_scripts": ["process=process.main:main"],
+        "numba_extensions": ["init = numba_scipy:_init_extension"],
+    },
     "extra_link_args": EXTRA_ARGS,
 }
 

From e12f430ab3beea69e7ab9c7109abc992545e51bd Mon Sep 17 00:00:00 2001
From: Jonathan Matthews <jonathan.matthews@ukaea.uk>
Date: Wed, 16 Apr 2025 09:58:39 +0100
Subject: [PATCH 5/7] param name changes from rebase

---
 process/pfcoil.py | 129 ++++++++++++++++++++--------------------------
 1 file changed, 57 insertions(+), 72 deletions(-)

diff --git a/process/pfcoil.py b/process/pfcoil.py
index aa0ea286d2..53c0e2cd5f 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -3,7 +3,9 @@
 
 import numba
 import numpy as np
-from scipy import optimize
+from numba.types import CPointer, float64, intc
+from scipy import LowLevelCallable, optimize
+from scipy.integrate import IntegrationWarning, nquad, quad
 from scipy.linalg import svd
 from scipy.special import ellipe, ellipk
 
@@ -24,17 +26,6 @@
 from process.fortran import tfcoil_variables as tfv
 from process.fortran import times_variables as tv
 from process.utilities.f2py_string_patch import f2py_compatible_to_string
-from process import fortran as ft
-import process.superconductors as superconductors
-import math
-import numpy as np
-import numba
-from numba.types import CPointer, float64, intc
-from scipy import LowLevelCallable
-from scipy.integrate import IntegrationWarning, quad, nquad
-import logging
-from scipy import optimize
-from scipy.special import ellipe, ellipk
 
 logger = logging.getLogger(__name__)
 
@@ -681,7 +672,9 @@ def pfcoil(self):
                 # Peak field
                 if ij == 0:
                     # Index args +1ed
-                    bri, bro, bzi, bzo = self.peakb(i, iii, it)  # returns b_pf_coil_peak, bpf2
+                    bri, bro, bzi, bzo = self.peakb(
+                        i, iii, it
+                    )  # returns b_pf_coil_peak, bpf2
 
                 # Issue 1871.  MDK
                 # Allowable current density (for superconducting coils) for each coil, index i
@@ -1171,24 +1164,28 @@ def ohcalc(self):
         # Peak field at the End-Of-Flattop (EOF)
         # Occurs at inner edge of coil; bmaxoh2 and bzi are of opposite sign at EOF
 
-        current = pfv.coheof * pfv.areaoh
+        current = pfv.j_cs_flat_top_end * pfv.a_cs_poloidal
         # Self field from CS
         br_inner, bz_inner, psi_inner = semianalytic(
-            pfv.rpf[pfv.nohc - 1],
-            pfv.zpf[pfv.nohc - 1],
-            pfv.ra[pfv.nohc - 1],
-            pfv.zpf[pfv.nohc - 1],
-            pfv.rb[pfv.nohc - 1] - pfv.ra[pfv.nohc - 1],
-            pfv.zh[pfv.nohc - 1] - pfv.zl[pfv.nohc - 1],
+            pfv.r_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
+            - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
+            - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1],
             current,
         )
         br_outer, bz_outer, psi_outer = semianalytic(
-            pfv.rpf[pfv.nohc - 1],
-            pfv.zpf[pfv.nohc - 1],
-            pfv.rb[pfv.nohc - 1],
-            pfv.zpf[pfv.nohc - 1],
-            pfv.rb[pfv.nohc - 1] - pfv.ra[pfv.nohc - 1],
-            pfv.zh[pfv.nohc - 1] - pfv.zl[pfv.nohc - 1],
+            pfv.r_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
+            - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
+            - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1],
             current,
         )
 
@@ -1203,26 +1200,30 @@ def ohcalc(self):
         bohco = abs(bzo + bz_outer)
 
         # Peak field at the Beginning-Of-Pulse (BOP)
-        # Occurs at inner edge of coil; bmaxoh0 and bzi are of same sign at BOP
+        # Occurs at inner edge of coil; b_cs_peak_pulse_start and bzi are of same sign at BOP
 
-        current = pfv.cohbop * pfv.areaoh
+        current = pfv.j_cs_pulse_start * pfv.a_cs_poloidal
         # Self field
         br_inner, bz_inner, psi_inner = semianalytic(
-            pfv.rpf[pfv.nohc - 1],
-            pfv.zpf[pfv.nohc - 1],
-            pfv.ra[pfv.nohc - 1],
-            pfv.zpf[pfv.nohc - 1],
-            pfv.rb[pfv.nohc - 1] - pfv.ra[pfv.nohc - 1],
-            pfv.zh[pfv.nohc - 1] - pfv.zl[pfv.nohc - 1],
+            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
+            - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
+            - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1],
             current,
         )
         br_outer, bz_outer, psi_outer = semianalytic(
-            pfv.rpf[pfv.nohc - 1],
-            pfv.zpf[pfv.nohc - 1],
-            pfv.rb[pfv.nohc - 1],
-            pfv.zpf[pfv.nohc - 1],
-            pfv.rb[pfv.nohc - 1] - pfv.ra[pfv.nohc - 1],
-            pfv.zh[pfv.nohc - 1] - pfv.zl[pfv.nohc - 1],
+            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
+            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
+            - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
+            pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
+            - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1],
             current,
         )
         timepoint = 2
@@ -1425,31 +1426,22 @@ def peakb(self, i, ii, it):
             # so exclude its own contribution; its self field is
             # dealt with inside ohcalc
             kk = 0
-            rp = np.array([pfv.ra[i], pfv.rb[i]])
-            zp = np.array([pfv.zpf[i], pfv.zpf[i]])
+            rp = np.array([pfv.r_pf_coil_inner[i], pfv.r_pf_coil_outer[i]])
+            zp = np.array([pfv.z_pf_coil_middle[i], pfv.z_pf_coil_middle[i]])
         else:
             # Check different times for maximum current
             if (
-                abs(
-                    pfv.c_pf_cs_coil_pulse_start_ma[i]
-                    - pfv.c_pf_cs_coils_peak_ma[i]
-                )
+                abs(pfv.c_pf_cs_coil_pulse_start_ma[i] - pfv.c_pf_cs_coils_peak_ma[i])
                 < 1.0e-12
             ):
                 it = 2
             elif (
-                abs(
-                    pfv.c_pf_cs_coil_flat_top_ma[i]
-                    - pfv.c_pf_cs_coils_peak_ma[i]
-                )
+                abs(pfv.c_pf_cs_coil_flat_top_ma[i] - pfv.c_pf_cs_coils_peak_ma[i])
                 < 1.0e-12
             ):
                 it = 4
             elif (
-                abs(
-                    pfv.c_pf_cs_coil_pulse_end_ma[i]
-                    - pfv.c_pf_cs_coils_peak_ma[i]
-                )
+                abs(pfv.c_pf_cs_coil_pulse_end_ma[i] - pfv.c_pf_cs_coils_peak_ma[i])
                 < 1.0e-12
             ):
                 it = 5
@@ -1509,7 +1501,7 @@ def peakb(self, i, ii, it):
             zc_array = np.append(zc_array, 0.0e0)
         bri, bzi, _ = bfield(rc_array, zc_array, current_array, rp[0], zp[0])
         bro, bzo, _ = bfield(rc_array, zc_array, current_array, rp[1], zp[1])
-        if not ((bv.iohcl != 0) and (i == pfv.nohc - 1)):
+        if not ((bv.iohcl != 0) and (i == pfv.n_cs_pf_coils - 1)):
             bri_self, bzi_self, _ = semianalytic(
                 self_rc, self_zc, rp[0], zp[0], dr, dz, self_current
             )
@@ -1527,15 +1519,8 @@ def peakb(self, i, ii, it):
             pf.rfxf[:kk],
             pf.zfxf[:kk],
             pf.cfxf[:kk],
-            pfv.r_pf_coil_inner[i - 1],
-            pfv.z_pf_coil_middle[i - 1],
-        )
-        pf.xind[:kk], bro, bzo, psi = bfield(
-            pf.rfxf[:kk],
-            pf.zfxf[:kk],
-            pf.cfxf[:kk],
-            pfv.rb[i - 1],
-            pfv.zpf[i - 1],
+            pfv.r_pf_coil_inner[i],
+            pfv.z_pf_coil_middle[i],
         )
 
         # b_pf_coil_peak and bpf2 for the Central Solenoid are calculated in OHCALC
@@ -3759,8 +3744,8 @@ def integrate(func, args, bound1, bound2):
         ]
         try:
             result = quad(func, bound1, bound2, args=args, points=points, limit=200)[0]
-        except IntegrationWarning:
-            raise IntegrationWarning
+        except IntegrationWarning as err:
+            raise IntegrationWarning from err
 
     # warnings.filterwarnings("default", category=IntegrationWarning)
     return result
@@ -3911,7 +3896,7 @@ def _partial_z_integrand_nojit(phi, rr, zz):
     # F2
     result = result - 0.5 * rr if rr - 1 < EPS else result - 0.5 / rr
     # F3
-    if 0.5 * np.pi * sin_phi > 1e-9:  # noqa: PLR2004
+    if 0.5 * np.pi * sin_phi > 1e-9:
         result -= sin_phi * np.arctan(zz * (rr - cos_phi) / (r0 * sin_phi))
     return result
 
@@ -3939,7 +3924,7 @@ def _full_psi_integrand(rp, phi, rc, zc, zp, d_rc, d_zc):
     """
     Integrand for psi = xBz
     """
-    z = zp - zc  # numba issue # noqa: PLR6104
+    z = zp - zc  # numba issue
     r1, r2 = (rc - d_rc) / rp, (rc + d_rc) / rp
     z1, z2 = (-d_zc - z) / rp, (d_zc - z) / rp
     return rp**2 * (
@@ -4071,18 +4056,18 @@ def semianalytic(rc, zc, rp, zp, d_rc, d_zc, current):
     if rp < 0.0001:
         rp = 0.0001
     r1, r2, z1, z2, j_tor = _get_working_coords(rc, zc, rp, zp, d_rc, d_zc)
-    Br = integrate(
+    br = integrate(
         _full_x_integrand, tuple(np.asarray((r1, r2, z1, z2)).ravel()), 0, np.pi
     )
 
-    Bz = _integrate_z_by_parts(r1, r2, z1, z2)
+    bz = _integrate_z_by_parts(r1, r2, z1, z2)
 
     psi = n_integrate(
         _full_psi_integrand,
         (floatify(rc), floatify(zc), floatify(zp), floatify(d_rc), floatify(d_zc)),
         [[0, floatify(rp)], [0, np.pi]],
     )
-    fac_B = 2e-7 * j_tor * rp  # MU_0/(2*np.pi)
+    fac_b = 2e-7 * j_tor * rp  # MU_0/(2*np.pi)
     fac_psi = 2e-7 * j_tor
 
-    return current * fac_B * Br, current * fac_B * Bz, current * fac_psi * psi
+    return current * fac_b * br, current * fac_b * bz, current * fac_psi * psi

From ff959c0556cb983aab0605f3ab1b07816c464957 Mon Sep 17 00:00:00 2001
From: Jonathan Matthews <jonathan.matthews@ukaea.uk>
Date: Wed, 16 Apr 2025 10:03:11 +0100
Subject: [PATCH 6/7] missed variable to remove

---
 process/pfcoil.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/process/pfcoil.py b/process/pfcoil.py
index 53c0e2cd5f..42eb1daff1 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -1518,7 +1518,6 @@ def peakb(self, i, ii, it):
         pf.xind[:kk] = mutual_inductance(
             pf.rfxf[:kk],
             pf.zfxf[:kk],
-            pf.cfxf[:kk],
             pfv.r_pf_coil_inner[i],
             pfv.z_pf_coil_middle[i],
         )

From 55b0d7d2c96abfad40e0a7f819fa12dfa9054c4a Mon Sep 17 00:00:00 2001
From: Jonathan Matthews <jonathan.matthews@ukaea.uk>
Date: Wed, 16 Apr 2025 11:35:30 +0100
Subject: [PATCH 7/7] made the half height inputs to semianalytic function
 actually half height not full

---
 process/pfcoil.py | 65 ++++++++++++++++++++++++++++++-----------------
 1 file changed, 42 insertions(+), 23 deletions(-)

diff --git a/process/pfcoil.py b/process/pfcoil.py
index 42eb1daff1..adfedc7c29 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -1171,10 +1171,16 @@ def ohcalc(self):
             pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
             pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
             pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
-            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
-            - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
-            pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
-            - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1],
+            0.5
+            * (
+                pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
+                - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1]
+            ),
+            0.5
+            * (
+                pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
+                - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1]
+            ),
             current,
         )
         br_outer, bz_outer, psi_outer = semianalytic(
@@ -1182,17 +1188,22 @@ def ohcalc(self):
             pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
             pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1],
             pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
-            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
-            - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
-            pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
-            - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1],
+            0.5
+            * (
+                pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
+                - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1]
+            ),
+            0.5
+            * (
+                pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
+                - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1]
+            ),
             current,
         )
 
         # Peak field due to other PF coils plus plasma
         timepoint = 5
         bri, bro, bzi, bzo = self.peakb(pfv.n_cs_pf_coils - 1, 99 - 1, timepoint)
-
         pfv.b_cs_peak_flat_top_end = abs(bzi + bz_inner)
 
         # Peak field on outboard side of central Solenoid
@@ -1209,10 +1220,16 @@ def ohcalc(self):
             pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
             pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
             pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
-            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
-            - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
-            pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
-            - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1],
+            0.5
+            * (
+                pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
+                - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1]
+            ),
+            0.5
+            * (
+                pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
+                - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1]
+            ),
             current,
         )
         br_outer, bz_outer, psi_outer = semianalytic(
@@ -1220,15 +1237,20 @@ def ohcalc(self):
             pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
             pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1],
             pfv.z_pf_coil_middle[pfv.n_cs_pf_coils - 1],
-            pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
-            - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1],
-            pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
-            - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1],
+            0.5
+            * (
+                pfv.r_pf_coil_outer[pfv.n_cs_pf_coils - 1]
+                - pfv.r_pf_coil_inner[pfv.n_cs_pf_coils - 1]
+            ),
+            0.5
+            * (
+                pfv.z_pf_coil_upper[pfv.n_cs_pf_coils - 1]
+                - pfv.z_pf_coil_lower[pfv.n_cs_pf_coils - 1]
+            ),
             current,
         )
         timepoint = 2
         bri, bro, bzi, bzo = self.peakb(pfv.n_cs_pf_coils - 1, 99 - 1, timepoint)
-
         pfv.b_cs_peak_pulse_start = abs(bz_inner + bzi)
 
         # Maximum field values
@@ -1332,7 +1354,6 @@ def ohcalc(self):
         if pfv.i_pf_conductor == 0:
             # Allowable coil overall current density at EOF
             # (superconducting coils only)
-
             (
                 jcritwp,
                 pfv.jcableoh_eof,
@@ -1363,7 +1384,6 @@ def ohcalc(self):
             pfv.j_cs_critical_flat_top_end = jcritwp * pfv.awpoh / pfv.a_cs_poloidal
 
             # Allowable coil overall current density at BOP
-
             (
                 jcritwp,
                 pfv.jcableoh_bop,
@@ -1483,8 +1503,8 @@ def peakb(self, i, ii, it):
                 if jj == i:
                     rp = np.array([pfv.r_pf_coil_inner[i], pfv.r_pf_coil_outer[i]])
                     zp = np.array([pfv.z_pf_coil_middle[jj], pfv.z_pf_coil_middle[jj]])
-                    dr = pfv.r_pf_coil_outer[i] - pfv.r_pf_coil_inner[i]
-                    dz = pfv.z_pf_coil_upper[jj] - pfv.z_pf_coil_lower[jj]
+                    dr = 0.5 * (pfv.r_pf_coil_outer[i] - pfv.r_pf_coil_inner[i])
+                    dz = 0.5 * (pfv.z_pf_coil_upper[jj] - pfv.z_pf_coil_lower[jj])
                     self_current = current
                     self_rc = pfv.r_pf_coil_middle[jj]
                     self_zc = pfv.z_pf_coil_middle[jj]
@@ -3087,7 +3107,6 @@ def superconpf(
                 arguments = (isumat, jsc, bmax, strain, bc20m, tc0m, c0)
             else:
                 arguments = (isumat, jsc, bmax, strain, bc20m, tc0m)
-
             another_estimate = 2 * thelium
             t_zero_margin, root_result = optimize.newton(
                 superconductors.current_density_margin,