! This module handles the loading, interpolation and evaluation of the
! MHD equilibrium data (poloidal current function, safety factor,
! poloidal flux, magnetic field, etc.)
!
module gray_equil
  use const_and_precisions, only : wp_, comp_huge
  use splines,              only : spline_simple, spline_1d, spline_2d, linear_1d
  use types,                only : contour

  implicit none

  ! macro for suppressing unused variable warnings
# define unused(x) associate(x => x); end associate

  type, abstract :: abstract_equil
    ! Generic equilibrium interface
    real(wp_) :: psi_a         = 0                        ! Poloidal flux at the edge minus flux on axis, ψ_a
    real(wp_) :: phi_a         = 0                        ! Toroidal flux at the edge (r=a), Φ_a
    real(wp_) :: b_axis        = 0                        ! Value of B_φ at the magnetic axis (used in J_cd def)
    real(wp_) :: b_centre      = 0                        ! Value of B_φ at R_centre (used in Jcd_astra def)
    real(wp_) :: r_centre      = 1                        ! Alternative reference radius for B_φ
    real(wp_) :: sgn_bphi      = 0                        ! Sign of B_φ (>0 means counter-clockwise)
    real(wp_) :: axis(2)       = [0, 0]                   ! Magnetic axis position (R₀, z₀)
    real(wp_) :: r_range(2)    = [-comp_huge, comp_huge]  ! R range of the equilibrium domain
    real(wp_) :: z_range(2)    = [-comp_huge, comp_huge]  ! z range of the equilibrium domain
    real(wp_) :: z_boundary(2) = [0, 0]                   ! z range of the plasma boundary
    contains
      procedure(pol_flux_sub),     deferred :: pol_flux
      procedure(pol_curr_sub),     deferred :: pol_curr
      procedure(safety_fun),       deferred :: safety
      procedure(rho_conv_fun),     deferred :: pol2tor, tor2pol
      procedure(flux_contour_sub), deferred :: flux_contour
      procedure                             :: b_field
      procedure                             :: tor_curr
  end type

  abstract interface
    pure subroutine pol_flux_sub(self, R, z, psi_n, dpsidr, dpsidz, &
                            ddpsidrr, ddpsidzz, ddpsidrz)
      ! Computes the normalised poloidal flux ψ_n and its
      ! derivatives wrt (R, z) up to the second order.
      !
      ! Note: all output arguments are optional.
      import :: abstract_equil, wp_
      class(abstract_equil), intent(in)            :: self
      real(wp_),             intent(in)            :: R, z
      real(wp_),             intent(out), optional :: &
        psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz
    end subroutine pol_flux_sub

    pure subroutine pol_curr_sub(self, psi_n, fpol, dfpol)
      ! Computes the poloidal current function F(ψ_n)
      ! and (optionally) its derivative dF/dψ_n given ψ_n.
      import :: abstract_equil, wp_
      class(abstract_equil), intent(in)            :: self
      real(wp_),             intent(in)            :: psi_n  ! normalised poloidal flux
      real(wp_),             intent(out)           :: fpol   ! poloidal current
      real(wp_),             intent(out), optional :: dfpol  ! derivative
    end subroutine pol_curr_sub

    pure function safety_fun(self, psi_n) result(q)
      ! Computes the safety factor q as a function of the
      ! normalised poloidal flux ψ_n.
      !
      ! Note: this returns the absolute value of q.
      import :: abstract_equil, wp_
      class(abstract_equil), intent(in) :: self
      real(wp_),             intent(in) :: psi_n
      real(wp_)                         :: q
    end function safety_fun

    pure function rho_conv_fun(self, rho_in) result(rho_out)
      ! Converts between poloidal (ρ_p) and toroidal (ρ_t) normalised radius
      import :: abstract_equil, wp_
      class(abstract_equil), intent(in) :: self
      real(wp_),             intent(in) :: rho_in
      real(wp_)                         :: rho_out
    end function rho_conv_fun

    subroutine flux_contour_sub(self, psi0, R_min, R, z, &
                                R_hi, z_hi, R_lo, z_lo)
      ! Computes a contour of the ψ(R,z)=ψ₀ flux surface.
      ! Notes:
      !  - R,z are the contour points
      !  - R_min is a value such that R>R_min for any contour point
      !  - (R,z)_hi and (R,z)_lo are a guess for the higher and lower
      !    horizontal-tangent points of the contour. These variables
      !    will be updated with the exact value on success.
      import :: abstract_equil, wp_
      class(abstract_equil), intent(in)    :: self
      real(wp_),             intent(in)    :: psi0
      real(wp_),             intent(in)    :: R_min
      real(wp_),             intent(out)   :: R(:), z(:)
      real(wp_),             intent(inout) :: R_hi, z_hi, R_lo, z_lo
    end subroutine flux_contour_sub
  end interface


  type, extends(abstract_equil) :: analytic_equil
    ! Analytical equilibrium
    private
    real(wp_) :: q0    ! Safety factor at the magnetic axis
    real(wp_) :: q1    ! Safety factor at the edge
    real(wp_) :: alpha ! Exponent for the q(ρ_p) power law
    real(wp_) :: R0    ! R of the magnetic axis (m)
    real(wp_) :: z0    ! z of the magnetic axis (m)
    real(wp_) ::  a    ! Minor radius (m)
    real(wp_) :: B0    ! Magnetic field at the magnetic axis (T)
  contains
    procedure :: pol_flux     => analytic_pol_flux
    procedure :: pol_curr     => analytic_pol_curr
    procedure :: safety       => analytic_safety
    procedure :: tor2pol      => analytic_tor2pol
    procedure :: pol2tor      => analytic_pol2tor
    procedure :: flux_contour => analytic_flux_contour
  end type


  type, extends(abstract_equil) :: vacuum
    ! Vacuum
    contains
      procedure :: pol_flux     => vacuum_pol_flux
      procedure :: pol_curr     => vacuum_pol_curr
      procedure :: safety       => vacuum_safety
      procedure :: tor2pol      => vacuum_conv
      procedure :: pol2tor      => vacuum_conv
      procedure :: flux_contour => vacuum_flux_contour
  end type


  type, extends(abstract_equil) :: numeric_equil
    ! Numerical equilibrium
    private
    real(wp_) :: fpol_a  ! Poloidal current at the edge (r=a), F_a
    ! Splines
    type(spline_2d)     :: psi_spline
    type(contour)       :: psi_domain
    type(spline_1d)     :: fpol_spline
    type(spline_simple) :: q_spline
    type(linear_1d)     :: rhop_spline, rhot_spline
    contains
      procedure :: pol_flux      => numeric_pol_flux
      procedure :: pol_curr      => numeric_pol_curr
      procedure :: safety        => numeric_safety
      procedure :: tor2pol       => numeric_tor2pol
      procedure :: pol2tor       => numeric_pol2tor
      procedure :: flux_contour  => numeric_flux_contour
      procedure :: init          => numeric_init
      procedure :: find_ox_point
      procedure :: find_htg_point
  end type


  type eqdsk_data
    ! MHD equilibrium data from G-EQDSK file format
    real(wp_), allocatable :: grid_r(:)     ! R values of the uniform grid
    real(wp_), allocatable :: grid_z(:)     ! z values of the uniform grid
    real(wp_), allocatable :: fpol(:)       ! Poloidal current function, F(ψ_n)
    real(wp_), allocatable :: q(:)          ! Safety factor, q(ψ_n)
    real(wp_), allocatable :: psi(:)        ! Normalised poloidal flux, ψ_n(R)
    real(wp_), allocatable :: psi_map(:,:)  ! Normalised poloidal flux 2D map, ψ(R, z)
    real(wp_)              :: psi_a         ! Poloidal flux at the edge minus flux on axis, ψ_a
    real(wp_)              :: r_ref         ! Reference R₀ (B = B₀R₀/R without the plasma)
    real(wp_)              :: axis(2)       ! Magnetic axis position (R₀, z₀)
    type(contour)          :: limiter       ! limiter contour (wall)
    type(contour)          :: boundary      ! boundary contour (plasma)
  end type


  private
  public abstract_equil                         ! The abstract equilibrium object
  public analytic_equil, numeric_equil, vacuum  ! Implementations
  public load_equil                             ! To load equilibrium from file
  public eqdsk_data                             ! G-EQDSK data structure
  public contour                                ! re-export contours eqdsk_data

contains

  pure subroutine b_field(self, R, z, B_R, B_z, B_phi)
    ! Computes the magnetic field as a function of
    ! (R, z) in cylindrical coordinates
    !
    ! Note: all output arguments are optional.

    ! subroutine arguments
    class(abstract_equil), intent(in)            :: self
    real(wp_),             intent(in)            :: R, z
    real(wp_),             intent(out), optional :: B_R, B_z, B_phi

    ! local variables
    real(wp_) :: psi_n, fpol, dpsidr, dpsidz

    call self%pol_flux(R, z, psi_n, dpsidr, dpsidz)
    call self%pol_curr(psi_n, fpol)

    ! The field in cocos=3 is given by
    !
    !   B = F(ψ)∇φ + ∇ψ×∇φ.
    !
    ! Writing the gradient of ψ=ψ(R,z) as
    !
    !   ∇ψ = ∂ψ/∂R ∇R + ∂ψ/∂z ∇z,
    !
    ! and carrying out the cross products:
    !
    !   B = F(ψ)∇φ - ∂ψ/∂z ∇R/R + ∂ψ/∂R ∇z/R
    !
    if (present(B_R))   B_R = - 1/R * dpsidz * self%psi_a
    if (present(B_z))   B_z = + 1/R * dpsidr * self%psi_a
    if (present(B_phi)) B_phi = fpol / R
  end subroutine b_field


  pure function tor_curr(self, R, z) result(J_phi)
    ! Computes the toroidal current J_φ as a function of (R, z)
    use const_and_precisions, only : mu0_

    ! function arguments
    class(abstract_equil), intent(in) :: self
    real(wp_),             intent(in) :: R, z
    real(wp_)                         :: J_phi

    ! local variables
    real(wp_) :: dB_Rdz, dB_zdR              ! derivatives of B_R, B_z
    real(wp_) :: dpsidr, ddpsidrr, ddpsidzz  ! derivatives of ψ_n

    call self%pol_flux(R, z, dpsidr=dpsidr, ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz)

    ! In the usual MHD limit we have ∇×B = μ₀J. Using the
    ! curl in cylindrical coords the toroidal current is
    !
    !   J_φ = 1/μ₀ (∇×B)_φ = 1/μ₀ [∂B_R/∂z - ∂B_z/∂R].
    !
    ! Finally, from B = F(ψ)∇φ + ∇ψ×∇φ we find:
    !
    !   B_R = - 1/R ∂ψ/∂z,
    !   B_z = + 1/R ∂ψ/∂R,
    !
    ! from which:
    !
    !   ∂B_R/∂z = - 1/R ∂²ψ/∂z²
    !   ∂B_z/∂R = + 1/R ∂²ψ/∂R² - 1/R² ∂ψ/∂R.
    !
    dB_Rdz = - 1/R * ddpsidzz * self%psi_a
    dB_zdR = + 1/R * (ddpsidrr - 1/R * dpsidr) * self%psi_a
    J_phi  = 1/mu0_ * (dB_Rdz - dB_zdR)
  end function tor_curr


  !
  ! Analytical model
  !

  pure subroutine analytic_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
                          ddpsidrr, ddpsidzz, ddpsidrz)
    use const_and_precisions, only : pi

    ! function arguments
    class(analytic_equil), intent(in)            :: self
    real(wp_),             intent(in)            :: R, z
    real(wp_),             intent(out), optional :: &
      psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz

    ! local variables
    real(wp_) :: r_g, rho_t, rho_p         ! geometric radius, √Φ_n, √ψ_n
    real(wp_) :: gamma                     ! γ = 1/√(1 - r²/R₀²)
    real(wp_) :: dpsidphi                  !  (∂ψ_n/∂Φ_n)
    real(wp_) :: ddpsidphidr, ddpsidphidz  ! ∇(∂ψ_n/∂Φ_n)
    real(wp_) :: phi_n                     !   Φ_n
    real(wp_) :: dphidr, dphidz            !  ∇Φ_n
    real(wp_) :: ddphidrdr, ddphidzdz      ! ∇∇Φ_n
    real(wp_) :: ddphidrdz                 ! ∂²Φ_n/∂R∂z
    real(wp_) :: q, dq                     ! q(ρ_p), Δq=(q₁-q₀)/(α/2 + 1)
    real(wp_) :: dqdr, dqdz                ! ∇q
    real(wp_) :: dphidr2, ddphidr2dr2      ! dΦ_n/d(r²), d²Φ_n/d(r²)²

    ! The normalised poloidal flux ψ_n(R, z) is computed as follows:
    !   1. ψ_n = ρ_p²
    !   2. ρ_p = ρ_p(ρ_t), using `self%tor2pol`, which in turns uses q(ψ)
    !   3. ρ_t = √Φ_n
    !   4. Φ_n = Φ(r_g)/Φ(a), where Φ(r) is the flux of B_φ=B₀R₀/R
    !      through a circular surface
    !   5. r_g = √[(R-R₀)²+(z-z₀)²] is the geometric minor radius
    r_g = hypot(R - self%R0, z - self%z0)

    ! The exact flux of the toroidal field B_φ = B₀R₀/R is:
    !
    !   Φ(r) = B₀πr² 2γ/(γ + 1) where γ=1/√(1 - r²/R₀²).
    !
    ! Notes:
    !   1. the function Φ(r) is defined for r≤R₀ only.
    !   2. as r → 0, γ → 1, so Φ ~ B₀πr².
    !   3. as r → 1⁻, Φ → 2B₀πr² but dΦ/dr → -∞.
    !   4. |B_R|, |B_z| → +-∞.
    !
    if (r_g > self%R0) then
      if (present(psi_n))    psi_n    = -1
      if (present(dpsidr))   dpsidr   = 0
      if (present(dpsidz))   dpsidz   = 0
      if (present(ddpsidrr)) ddpsidrr = 0
      if (present(ddpsidzz)) ddpsidzz = 0
      if (present(ddpsidrz)) ddpsidrz = 0
      return
    end if

    gamma = 1 / sqrt(1 - (r_g/self%R0)**2)
    phi_n = self%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / self%phi_a
    rho_t = sqrt(phi_n)
    rho_p = self%tor2pol(rho_t)

    ! For ∇Φ_n and ∇∇Φ_n we also need:
    !
    !   ∂Φ∂(r²)   = B₀π γ(r)
    !   ∂²Φ∂(r²)² = B₀π γ³(r) / (2 R₀²)
    !
    dphidr2     = self%B0 * pi * gamma / self%phi_a
    ddphidr2dr2 = self%B0 * pi * gamma**3/(2 * self%R0**2) / self%phi_a

    ! ∇Φ_n = ∂Φ_n/∂(r²) ∇(r²)
    ! where ∇(r²) = 2[(R-R₀), (z-z₀)]
    dphidr = dphidr2 * 2*(R - self%R0)
    dphidz = dphidr2 * 2*(z - self%z0)

    ! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²)
    !       = ∂²Φ_n/∂(r²)² ∇(r²)∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²)
    ! where ∇∇(r²) = 2I
    ddphidrdr = ddphidr2dr2 * 4*(R - self%R0)*(R - self%R0) + dphidr2*2
    ddphidzdz = ddphidr2dr2 * 4*(z - self%z0)*(z - self%z0) + dphidr2*2
    ddphidrdz = ddphidr2dr2 * 4*(R - self%R0)*(z - self%z0)

    ! ψ_n = ρ_p(ρ_t)²
    if (present(psi_n)) psi_n = rho_p**2

    ! Using the definitions in `frhotor`:
    !
    !   ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n
    !
    !   ∂ψ_n/∂Φ_n = Φ_a/ψ_a ∂ψ/∂Φ
    !             = Φ_a/ψ_a 1/2πq
    !
    ! Using ψ_a = 1/2π Φ_a / (q₀ + Δq), then:
    !
    !   ∂ψ_n/∂Φ_n = (q₀ + Δq)/q
    !
    q  = self%q0 + (self%q1 - self%q0) * rho_p**self%alpha
    dq = (self%q1 - self%q0) / (self%alpha/2 + 1)
    dpsidphi = (self%q0 + dq) / q

    ! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n
    if (present(dpsidr)) dpsidr = dpsidphi * dphidr
    if (present(dpsidz)) dpsidz = dpsidphi * dphidz

    ! For the second derivatives:
    !
    !   ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n
    !
    !   ∇(∂ψ_n/∂Φ_n) = - (∂ψ_n/∂Φ_n) ∇q/q
    !
    ! From q(ψ) = q₀ + (q₁-q₀) ψ_n^α/2, we have:
    !
    !   ∇q = α/2 (q-q₀) ∇ψ_n/ψ_n
    !      = α/2 (q-q₀)/ψ_n (∂ψ_n/∂Φ_n) ∇Φ_n.
    !
    dqdr = self%alpha/2 * (self%q1 - self%q0)*rho_p**(self%alpha-2) * dpsidphi * dphidr
    dqdz = self%alpha/2 * (self%q1 - self%q0)*rho_p**(self%alpha-2) * dpsidphi * dphidz
    ddpsidphidr = - dpsidphi * dqdr/q
    ddpsidphidz = - dpsidphi * dqdz/q

    ! Combining all of the above:
    !
    !   ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n
    !
    if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * ddphidrdr
    if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * ddphidzdz
    if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz
  end subroutine analytic_pol_flux


  pure subroutine analytic_pol_curr(self, psi_n, fpol, dfpol)
    ! subroutine arguments
    class(analytic_equil), intent(in)            :: self
    real(wp_),             intent(in)            :: psi_n
    real(wp_),             intent(out)           :: fpol   ! poloidal current
    real(wp_),             intent(out), optional :: dfpol  ! its derivative

    unused(psi_n)

    ! The poloidal current function F(ψ) is just a constant:
    !
    !  B_φ = B₀R₀/R φ^ ≡ F(ψ)∇φ  ⇒  F(ψ)=B₀R₀
    !
    fpol = self%B0 * self%R0
    if (present(dfpol)) dfpol = 0
  end subroutine analytic_pol_curr


  pure function analytic_safety(self, psi_n) result(q)
    ! function arguments
    class(analytic_equil), intent(in) :: self
    real(wp_),             intent(in) :: psi_n
    real(wp_)                         :: q

    ! local variables
    real(wp_) :: rho_p

    ! The safety factor is a power law in ρ_p:
    !
    !   q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α
    !
    rho_p = sqrt(psi_n)
    q = abs(self%q0 + (self%q1 - self%q0) * rho_p**self%alpha)
  end function analytic_safety


  pure function analytic_pol2tor(self, rho_in) result(rho_out)
    ! function arguments
    class(analytic_equil), intent(in) :: self
    real(wp_),             intent(in) :: rho_in
    real(wp_)                         :: rho_out

    ! local variables
    real(wp_) :: dq

    ! The change of variable is obtained by integrating
    !
    !   q(ψ) = 1/2π ∂Φ/∂ψ
    !
    ! and defining ψ = ψ_a ρ_p²,  Φ = Φ_a ρ_t².
    ! The result is:
    !
    !   - ψ_a = 1/2π Φ_a / [q₀ + Δq]
    !
    !   - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)]
    !
    ! where Δq = (q₁ - q₀)/(α/2 + 1)
    dq = (self%q1 - self%q0) / (self%alpha/2 + 1)
    rho_out = rho_in * sqrt((self%q0 + dq*rho_in**self%alpha) / (self%q0 + dq))
  end function analytic_pol2tor


  pure function analytic_tor2pol(self, rho_in) result(rho_out)
    ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius
    use const_and_precisions, only : comp_eps
    use minpack,              only : hybrj1

    ! function arguments
    class(analytic_equil), intent(in) :: self
    real(wp_),             intent(in) :: rho_in
    real(wp_)                         :: rho_out

    ! local variables
    real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7)
    integer   :: info

    ! In general there is no closed form for ρ_p(ρ_t) in the
    ! analytical model, we thus solve numerically the equation
    ! ρ_t(ρ_p) = ρ_t₀ for ρ_p.
    rho_p = [rho_in]  ! first guess, ρ_p ≈ ρ_t
    call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, &
                ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7)
    rho_out = rho_p(1)

    contains

      pure subroutine equation(n, x, f, df, ldf, flag)
        ! The equation to solve: f(x) = ρ_t(x) - ρ_t₀ = 0

        ! optimal step size
        real(wp_), parameter :: e = comp_eps**(1/3.0_wp_)

        ! subroutine arguments
        integer,   intent(in)    :: n, ldf, flag
        real(wp_), intent(in)    :: x(n)
        real(wp_), intent(inout) :: f(n), df(ldf,n)

        if (flag == 1) then
          ! return f(x)
          f(1) = self%pol2tor(x(1)) - rho_in
        else
          ! return f'(x), computed numerically
          if (x(1) - e > 0) then
            df(1,1) = (self%pol2tor(x(1) + e) - self%pol2tor(x(1) - e)) / (2*e)
          else
            ! single-sided when close to ρ=0
            df(1,1) = (self%pol2tor(x(1) + e) - self%pol2tor(x(1))) / e
          end if
        end if
      end subroutine

  end function analytic_tor2pol


  pure subroutine analytic_flux_contour(self, psi0, R_min, R, z, &
                                  R_hi, z_hi, R_lo, z_lo)
    use const_and_precisions, only : pi

    ! subroutine arguments
    class(analytic_equil), intent(in)    :: self
    real(wp_),             intent(in)    :: psi0
    real(wp_),             intent(in)    :: R_min
    real(wp_),             intent(out)   :: R(:), z(:)
    real(wp_),             intent(inout) :: R_hi, z_hi, R_lo, z_lo

    ! local variables
    integer   :: n, i
    real(wp_) :: r_g   ! geometric minor radius
    real(wp_) :: theta ! geometric poloidal angle

    unused(R_min)
    unused(R_hi); unused(z_hi); unused(R_lo); unused(z_lo)

    n     = size(R)
    r_g   = sqrt(psi0) * self%a
    theta = 2*pi / (n - 1)

    do concurrent (i=1:n)
      R(i) = self%R0 + r_g * cos(theta*(i-1))
      z(i) = self%z0 + r_g * sin(theta*(i-1))
    end do

  end subroutine analytic_flux_contour


  !
  ! Numerical equilibrium
  !

  pure subroutine numeric_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
                          ddpsidrr, ddpsidzz, ddpsidrz)
    ! function arguments
    class(numeric_equil), intent(in)            :: self
    real(wp_),            intent(in)            :: R, z
    real(wp_),            intent(out), optional :: &
      psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz

    if (self%psi_domain%contains(R, z)) then
      ! Within the safety domain
      if (present(psi_n))    psi_n    = self%psi_spline%eval(R, z)
      if (present(dpsidr))   dpsidr   = self%psi_spline%deriv(R, z, 1, 0)
      if (present(dpsidz))   dpsidz   = self%psi_spline%deriv(R, z, 0, 1)
      if (present(ddpsidrr)) ddpsidrr = self%psi_spline%deriv(R, z, 2, 0)
      if (present(ddpsidzz)) ddpsidzz = self%psi_spline%deriv(R, z, 0, 2)
      if (present(ddpsidrz)) ddpsidrz = self%psi_spline%deriv(R, z, 1, 1)
    else
      ! Values for outside the domain
      if (present(psi_n))    psi_n    = -1
      if (present(dpsidr))   dpsidr   = 0
      if (present(dpsidz))   dpsidz   = 0
      if (present(ddpsidrr)) ddpsidrr = 0
      if (present(ddpsidzz)) ddpsidzz = 0
      if (present(ddpsidrz)) ddpsidrz = 0
    end if

  end subroutine numeric_pol_flux


  pure subroutine numeric_pol_curr(self, psi_n, fpol, dfpol)
    ! subroutine arguments
    class(numeric_equil), intent(in)            :: self
    real(wp_),            intent(in)            :: psi_n
    real(wp_),            intent(out)           :: fpol   ! poloidal current
    real(wp_),            intent(out), optional :: dfpol  ! its derivative

    if (psi_n <= 1 .and. psi_n >= 0) then
      ! Inside plasma
      fpol = self%fpol_spline%eval(psi_n)
      if (present(dfpol)) dfpol = self%fpol_spline%deriv(psi_n)
    else
      ! Outside plasma
      fpol = self%fpol_a
      if (present(dfpol)) dfpol = 0
    end if
  end subroutine numeric_pol_curr
  

  pure function numeric_safety(self, psi_n) result(q)
    ! function arguments
    class(numeric_equil), intent(in) :: self
    real(wp_),            intent(in) :: psi_n
    real(wp_)                        :: q

    if (psi_n < 1) then
      ! Inside plasma
      q = self%q_spline%eval(psi_n)
    else
      ! Outside plasma, q is undefined
      q = 0
    end if
  end function numeric_safety


  pure function numeric_pol2tor(self, rho_in) result(rho_out)
    ! function arguments
    class(numeric_equil), intent(in) :: self
    real(wp_),            intent(in) :: rho_in
    real(wp_)                        :: rho_out

    rho_out = self%rhot_spline%eval(rho_in)
  end function numeric_pol2tor


  pure function numeric_tor2pol(self, rho_in) result(rho_out)
    ! function arguments
    class(numeric_equil), intent(in) :: self
    real(wp_),            intent(in) :: rho_in
    real(wp_)                        :: rho_out

    rho_out = self%rhop_spline%eval(rho_in)
  end function numeric_tor2pol


  subroutine numeric_flux_contour(self, psi0, R_min, R, z, &
                                 R_hi, z_hi, R_lo, z_lo)
    use const_and_precisions, only : pi
    use logger,               only : log_warning
    use dierckx,              only : profil, sproota

    ! subroutine arguments
    class(numeric_equil), intent(in)    :: self
    real(wp_),            intent(in)    :: psi0
    real(wp_),            intent(in)    :: R_min
    real(wp_),            intent(out)   :: R(:), z(:)
    real(wp_),            intent(inout) :: R_hi, z_hi, R_lo, z_lo

    ! local variables
    integer   :: n, np, i
    integer   :: ier, iopt, m
    real(wp_) :: theta
    real(wp_) :: R_hi1, z_hi1, R_lo1, z_lo1, zc
    real(wp_) :: czc(self%psi_spline%nknots_x), zeroc(4)
    character(256) :: msg

    n = size(R)
    np = (n - 1)/2
    theta = pi / np

    call self%find_htg_point(R_hi, z_hi, R_hi1, z_hi1, psi0)
    call self%find_htg_point(R_lo, z_lo, R_lo1, z_lo1, psi0)

    R(1) = R_lo1
    z(1) = z_lo1
    R(n) = R_lo1
    z(n) = z_lo1
    R(np+1) = R_hi1
    z(np+1) = z_hi1

    do i = 2, np
      zc = z_lo1 + (z_hi1 - z_lo1) * (1 - cos(theta*(i-1))) / 2
      iopt = 1
      associate (s => self%psi_spline)
        call profil(iopt, s%knots_x, s%nknots_x, s%knots_y, s%nknots_y, &
                    s%coeffs, 3, 3, zc, s%nknots_x, czc, ier)
        if (ier > 0) then
          write(msg, '(a, a, g0)') &
            'when computing ψ(R,z) contour `profil` returned ier=', ier
          call log_warning(msg, mod='gray_equil', proc='numeric_flux_contour')
        end if
        call sproota(psi0, s%knots_x, s%nknots_x, czc, zeroc, 4, m, ier)
      end associate

      if (zeroc(1) > R_min) then
        R(i) = zeroc(1)
        R(n+1-i) = zeroc(2)
      else
         R(i) = zeroc(2)
         R(n+1-i) = zeroc(3)
      end if
      z(i) = zc
      z(n+1-i) = zc
    end do

    ! Replace the initial guess with the exact values
    R_hi = R_hi1
    z_hi = z_hi1
    R_lo = R_lo1
    z_lo = z_lo1
  end subroutine numeric_flux_contour


  subroutine numeric_init(self, params, data, err)
    ! Initialises a numeric equilibrium
    use const_and_precisions, only : zero, one
    use gray_params,          only : equilibrium_parameters
    use gray_params,          only : EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
    use gray_params,          only : X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING
    use utils,                only : vmaxmin
    use logger,               only : log_debug, log_info

    ! subroutine arguments
    class(numeric_equil),         intent(out) :: self
    type(equilibrium_parameters), intent(in)  :: params
    type(eqdsk_data),             intent(in)  :: data
    integer,                      intent(out) :: err

    ! local variables
    integer :: nr, nz, nrz, npsi, nbnd, ibinf, ibsup
    real(wp_) :: psinoptmp, psinxptmp
    real(wp_) :: rbinf, rbsup, zbinf, zbsup, R1, z1
    real(wp_) :: psiant, psinop
    real(wp_),  dimension(:),  allocatable :: rv1d, zv1d, fvpsi
    integer :: i, j, ij
    character(256) :: msg ! for log messages formatting

    ! compute array sizes
    nr    = size(data%grid_r)
    nz    = size(data%grid_z)
    nrz   = nr*nz
    npsi  = size(data%psi)

    ! length in m !!!
    self%r_range = [data%grid_r(1), data%grid_r(nr)]
    self%z_range = [data%grid_z(1), data%grid_z(nz)]

    call log_debug('computing splines...', mod='gray_equil', proc='numeric_init')

    ! Spline interpolation of ψ(R, z)

    select case (params%iequil)

      case (EQ_EQDSK_PARTIAL)
        ! Data valid only inside boundary (data%psi=0 outside),
        ! presence of boundary anticipated here to filter invalid data
        nbnd = size(data%boundary%R)

        ! allocate knots and spline coefficients arrays
        allocate(self%psi_spline%knots_x(nr + 4), self%psi_spline%knots_y(nz + 4))
        allocate(self%psi_spline%coeffs(nrz))

        ! determine number of valid grid points
        nrz=0
        do j=1,nz
          do i=1,nr
            if (nbnd > 0) then
              if (.not. data%boundary%contains(data%grid_r(i), data%grid_z(j))) cycle
            else
              if (data%psi_map(i,j) <= 0) cycle
            end if
            nrz=nrz+1
          end do
        end do

        ! store valid data
        allocate(rv1d(nrz), zv1d(nrz), fvpsi(nrz))
        ij=0
        do j=1,nz
          do i=1,nr
            if (nbnd > 0) then
              if (.not. data%boundary%contains(data%grid_r(i), data%grid_z(j))) cycle
            else
              if (data%psi_map(i,j) <= 0) cycle
            end if
            ij = ij + 1
            rv1d(ij)  = data%grid_r(i)
            zv1d(ij)  = data%grid_z(j)
            fvpsi(ij) = data%psi_map(i,j)
          end do
        end do

        ! Fit as a scattered set of points
        ! use reduced number of knots to limit memory comsumption ?
        self%psi_spline%nknots_x=nr/4+4
        self%psi_spline%nknots_y=nz/4+4
        call self%psi_spline%init_nonreg( &
          rv1d, zv1d, fvpsi, tension=params%ssplps, &
          range=[self%r_range, self%z_range], err=err)

        ! if failed, re-fit with an interpolating spline (zero tension)
        if (err == -1) then
          err = 0
          self%psi_spline%nknots_x=nr/4+4
          self%psi_spline%nknots_y=nz/4+4
          call self%psi_spline%init_nonreg( &
            rv1d, zv1d, fvpsi, tension=zero, &
            range=[self%r_range, self%z_range], err=err)
        end if
        deallocate(rv1d, zv1d, fvpsi)
        ! reset nrz to the total number of grid points for next allocations
        nrz = nr*nz

    case (EQ_EQDSK_FULL)
      ! Data are valid on the full R,z grid

      ! reshape 2D ψ array to 1D (transposed)
      allocate(fvpsi(nrz))
      fvpsi = reshape(transpose(data%psi_map), [nrz])

      ! compute spline coefficients
      call self%psi_spline%init(data%grid_r, data%grid_z, fvpsi, nr, nz, &
                                range=[self%r_range, self%z_range],  &
                                tension=params%ssplps, err=err)

      ! if failed, re-fit with an interpolating spline (zero tension)
      if (err == -1) then
        call self%psi_spline%init(data%grid_r, data%grid_z, fvpsi, nr, nz, &
                                  range=[self%r_range, self%z_range], tension=zero)
        err = 0
      end if
      deallocate(fvpsi)
    end select

    if (err /= 0) then
      err = 2
      return
    end if

    ! compute spline coefficients for ψ(R,z) partial derivatives
    call self%psi_spline%init_deriv(nr, nz, 1, 0)  ! ∂ψ/∂R
    call self%psi_spline%init_deriv(nr, nz, 0, 1)  ! ∂ψ/∂z
    call self%psi_spline%init_deriv(nr, nz, 1, 1)  ! ∂²ψ/∂R∂z
    call self%psi_spline%init_deriv(nr, nz, 2, 0)  ! ∂²ψ/∂R²
    call self%psi_spline%init_deriv(nr, nz, 0, 2)  ! ∂²ψ/∂z²

    ! set the initial ψ(R,z) domain to the grid boundary
    !
    ! Note: this is required to bootstrap `flux_pol` calls
    ! within this very subroutine.
    self%psi_domain = contour(Rmin=self%r_range(1), Rmax=self%r_range(2), &
                              zmin=self%z_range(1), zmax=self%z_range(2))

    ! Spline interpolation of F(ψ)

    ! give a small weight to the last point
    block
      real(wp_) :: w(npsi)
      w(1:npsi-1) = 1
      w(npsi)     = 1.0e2_wp_
      call self%fpol_spline%init(data%psi, data%fpol, npsi, weights=w, &
                                 range=[zero, one], tension=params%ssplf)
    end block

    ! set vacuum value used outside 0 ≤ ψ ≤ 1 range
    self%fpol_a   = self%fpol_spline%eval(data%psi(npsi))
    self%sgn_bphi = sign(one, self%fpol_a)

    ! Re-normalize ψ_n after the spline computation
    ! Note: this ensures 0 ≤ ψ_n'(R,z) < 1 inside the plasma

    ! Start with un-corrected ψ_n
    self%psi_a = data%psi_a
    psinop     = 0  ! ψ_n(O point)
    psiant     = 1  ! ψ_n(X point) - ψ_n(O point)

    ! Use provided boundary to set an initial guess
    ! for the search of O/X points
    nbnd = size(data%boundary%R)
    if (nbnd > 0) then
      call vmaxmin(data%boundary%z, zbinf, zbsup, ibinf, ibsup)
      rbinf = data%boundary%R(ibinf)
      rbsup = data%boundary%R(ibsup)
    else
      zbinf = data%grid_z(2)
      zbsup = data%grid_z(nz-1)
      rbinf = data%grid_r((nr+1)/2)
      rbsup = rbinf
    end if

    ! Search for exact location of the magnetic axis
    call self%find_ox_point(R0=data%axis(1), z0=data%axis(2), &
                            R1=self%axis(1), z1=self%axis(2), psi1=psinoptmp)

    write (msg, '("O-point found:", 3(x,a,"=",g0.3))') &
      'r', self%axis(1), 'z', self%axis(2), 'ψ', psinoptmp
    call log_info(msg, mod='gray_equil', proc='numeric_init')

    ! Search for X-point
    select case (params%ixp)
      case (X_AT_BOTTOM)
        call self%find_ox_point(R0=rbinf, z0=zbinf, R1=R1, z1=z1, psi1=psinxptmp)
        rbinf = z1
        if (psinxptmp /= -1) then
          write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
            'R', R1, 'z', z1, 'ψ', psinxptmp
          call log_info(msg, mod='gray_equil', proc='numeric_init')
          psinop = psinoptmp
          psiant = psinxptmp-psinop
          call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbsup)/2, &
                                   R1=r1, z1=zbsup, psi0=one)
        end if

      case (X_AT_TOP)
        call self%find_ox_point(R0=rbsup, z0=zbsup, R1=R1, z1=z1, psi1=psinxptmp)
        zbsup = z1
        if (psinxptmp /= -1) then
          write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
            'R', r1, 'z', zbsup, 'ψ', psinxptmp
          call log_info(msg, mod='gray_equil', proc='numeric_init')
          psinop = psinoptmp
          psiant = psinxptmp - psinop
          call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbinf)/2, &
                                   R1=r1, z1=zbinf, psi0=one)
        end if

      case (X_IS_MISSING)
        psinop = psinoptmp
        psiant = 1 - psinop
        ! Find upper horizontal tangent point
        call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbsup)/2, &
                                 R1=rbsup, z1=zbsup, psi0=one)
        ! Find lower horizontal tangent point
        call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbinf)/2, &
                                 R1=rbinf, z1=zbinf, psi0=one)
        write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') &
          'R', rbinf, rbsup, 'z', zbinf, zbsup
        call log_info(msg, mod='gray_equil', proc='numeric_init')
    end select

    ! Correct the spline coefficients: ψ_n → (ψ_n - psinop)/psiant
    call self%psi_spline%transform(1/psiant, -psinop/psiant)

    ! Do the opposite scaling to preserve un-normalised values
    ! Note: this is only used for the poloidal magnetic field
    self%psi_a = self%psi_a * psiant

    ! Compute other constants
    call self%pol_curr(zero, self%b_axis)
    self%b_axis     = self%b_axis / self%axis(1)
    self%b_centre   = self%fpol_a / data%r_ref
    self%r_centre   = data%r_ref
    self%z_boundary = [zbinf, zbsup]

    write (msg, '(2(a,g0.3))') 'B_centre=', self%b_centre, ' B_axis=', self%b_axis
    call log_info(msg, mod='gray_equil', proc='numeric_init')

    ! Compute ρ_p/ρ_t mapping based on the input q profile
    block
      use const_and_precisions, only : pi
      real(wp_), dimension(npsi) :: phi, rho_p, rho_t
      real(wp_) :: dx
      integer   :: k

      call self%q_spline%init(data%psi, abs(data%q), npsi)

      ! Toroidal flux as Φ(ψ) = 2π ∫q(ψ)dψ
      phi(1) = 0
      do k = 1, npsi-1
        dx = self%q_spline%data(k+1) - self%q_spline%data(k)
        phi(k+1) = phi(k) + dx*(self%q_spline%coeffs(k,1)   + dx*(self%q_spline%coeffs(k,2)/2 +  &
                            dx*(self%q_spline%coeffs(k,3)/3 + dx* self%q_spline%coeffs(k,4)/4) ) )
      end do

      self%phi_a = phi(npsi)
      rho_p      = sqrt(data%psi)
      rho_t      = sqrt(phi/self%phi_a)
      self%phi_a = 2*pi * abs(self%psi_a) * self%phi_a

      call self%rhop_spline%init(rho_t, rho_p, size(rho_p))
      call self%rhot_spline%init(rho_p, rho_t, size(rho_t))
    end block

    call log_debug('splines computed', mod='gray_equil', proc='numeric_init')

    ! Compute the domain of the ψ mapping
    self%psi_domain = data%boundary
    call rescale_boundary(self%psi_domain, self%psi_spline, O=self%axis, t0=1.5_wp_)

  end subroutine numeric_init


  pure subroutine rescale_boundary(cont, psi, O, t0)
    ! Given the plasma boundary contour `cont`, the position of the
    ! magnetic axis `O`, and a scaling factor `t0`; this subroutine
    ! rescales the contour by `t0` about `O` while ensuring the
    ! psi_spline stays monotonic within the new boundary.

    ! subroutine arguments
    type(contour),   intent(inout) :: cont ! (R,z) contour
    type(spline_2d), intent(inout) :: psi  ! ψ(R,z) spline
    real(wp_),       intent(in)    :: O(2) ! center point
    real(wp_),       intent(in)    :: t0   ! scaling factor

    ! subroutine variables
    integer :: i
    real(wp_) :: t
    real(wp_), parameter :: dt = 0.05
    real(wp_) :: P(2), N(2)

    do i = 1, size(cont%R)
      ! For each point on the contour compute:
      P = [cont%R(i), cont%z(i)]  ! point on the contour
      N = P - O                   ! direction of the line from O to P

      ! Find the max t: s(t) = ψ(O + tN) is monotonic in [1, t]
      t = 1
      do while (t < t0)
        if (s(t + dt) < s(t)) exit
        t = t + dt
      end do

      ! The final point is thus O + tN
      P = O + t * N
      cont%R(i) = P(1)
      cont%z(i) = P(2)
    end do

  contains

    pure function s(t)
      ! Rescriction of ψ(R, z) on the line Q(t) = O + tN

      real(wp_), intent(in) :: t
      real(wp_)             :: s, Q(2)

      Q = O + t * N
      s = psi%eval(Q(1), Q(2))
    end function

  end subroutine rescale_boundary


  subroutine find_ox_point(self, R0, z0, R1, z1, psi1)
    ! Given the point (R₀,z₀) as an initial guess, finds
    ! the exact location (R₁,z₁) where ∇ψ(R₁,z₁) = 0.
    ! It also returns ψ₁=ψ(R₁,z₁).
    !
    ! Note: this is used to find the magnetic X and O point,
    ! because both are stationary points for ψ(R,z).
    use const_and_precisions, only : comp_eps
    use minpack,              only : hybrj1
    use logger,               only : log_error, log_debug

    ! subroutine arguments
    class(numeric_equil)   :: self
    real(wp_), intent(in)  :: R0, z0
    real(wp_), intent(out) :: R1, z1, psi1

    ! local variables
    integer :: info
    real(wp_) :: sol(2), f(2), df(2,2), wa(15)
    character(256) :: msg

    sol = [R0, z0]  ! first guess

    call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, &
                tol=sqrt(comp_eps), info=info, wa=wa, lwa=15)
    if (info /= 1) then
      write (msg, '("guess: ", g0.3, ", ", g0.3)') R0, z0
      call log_debug(msg, mod='equilibrium', proc='find_ox_point')
      write (msg, '("solution: ", g0.3, ", ", g0.3)') sol
      call log_debug(msg, mod='equilibrium', proc='find_ox_point')
      write (msg, '("hybrj1 failed with error ",g0)') info
      call log_error(msg, mod='equilibrium', proc='find_ox_point')
    end if

    R1 = sol(1)  ! solution
    z1 = sol(2)
    call self%pol_flux(R1, z1, psi1)

  contains

    pure subroutine equation(n, x, f, df, ldf, flag)
      ! The equation to solve: f(R,z) = ∇ψ(R,z) = 0

      ! subroutine arguments
      integer,   intent(in)    :: n, flag, ldf
      real(wp_), intent(in)    :: x(n)
      real(wp_), intent(inout) :: f(n), df(ldf,n)

      if (flag == 1) then
        ! return f(R,z) = ∇ψ(R,z)
        call self%pol_flux(R=x(1), z=x(2), dpsidr=f(1), dpsidz=f(2))
      else
        ! return ∇f(R,z) = ∇∇ψ(R,z)
        call self%pol_flux(R=x(1), z=x(2), ddpsidrr=df(1,1), &
                           ddpsidzz=df(2,2), ddpsidrz=df(1,2))
        df(2,1) = df(1,2)
      end if
    end subroutine equation

  end subroutine find_ox_point


  subroutine find_htg_point(self, R0, z0, R1, z1, psi0)
    ! Given the point (R₀,z₀) as an initial guess, finds
    ! the exact location (R₁,z₁) where:
    !   {     ψ(R₁,z₁) = ψ₀
    !   { ∂ψ/∂R(R₁,z₁) = 0  .
    !
    ! Note: this is used to find the horizontal tangent
    ! point of the contour ψ(R,z)=ψ₀.
    use const_and_precisions, only : comp_eps
    use minpack,              only : hybrj1
    use logger,               only : log_error, log_debug

    ! subroutine arguments
    class(numeric_equil)   :: self
    real(wp_), intent(in)  :: R0, z0, psi0
    real(wp_), intent(out) :: R1, z1

    ! local variables
    integer :: info
    real(wp_) :: sol(2), f(2), df(2,2), wa(15)
    character(256) :: msg

    sol = [R0, z0]  ! first guess

    call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, &
                tol=sqrt(comp_eps), info=info, wa=wa, lwa=15)
    if (info /= 1) then
      write (msg, '("guess: ", g0.3, ", ", g0.3)') R0, z0
      call log_debug(msg, mod='equilibrium', proc='find_ox_point')
      write (msg, '("solution: ", g0.3, ", ", g0.3)') sol
      call log_debug(msg, mod='equilibrium', proc='find_htg_point')
      write (msg, '("hybrj1 failed with error ",g0)') info
      call log_error(msg, mod='equilibrium', proc='find_htg_point')
    end if

    R1 = sol(1)  ! solution
    z1 = sol(2)
  contains

    pure subroutine equation(n, x, f, df, ldf, flag)
      ! The equation to solve: f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R] = 0

      ! subroutine arguments
      integer,   intent(in)    :: n, ldf, flag
      real(wp_), intent(in)    :: x(n)
      real(wp_), intent(inout) :: f(n), df(ldf, n)

      if (flag == 1) then
        ! return f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R]
        call self%pol_flux(R=x(1), z=x(2), psi_n=f(1), dpsidr=f(2))
        f(1) = f(1) - psi0
      else
        ! return ∇f(R,z) = [[∂ψ/∂R, ∂ψ/∂z], [∂²ψ/∂R², ∂²ψ/∂R∂z]]
        call self%pol_flux(R=x(1), z=x(2), dpsidr=df(1,1), dpsidz=df(1,2), &
                           ddpsidrr=df(2,1), ddpsidrz=df(2,2))
      end if
    end subroutine equation

  end subroutine find_htg_point


  !
  ! Vacuum
  !

  pure subroutine vacuum_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
                          ddpsidrr, ddpsidzz, ddpsidrz)
    ! function arguments
    class(vacuum), intent(in)            :: self
    real(wp_),     intent(in)            :: R, z
    real(wp_),     intent(out), optional :: &
      psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz

    unused(self); unused(R); unused(z)

    ! ψ(R,z) is undefined everywhere, return -1
    if (present(psi_n))    psi_n    = -1
    if (present(dpsidr))   dpsidr   = 0
    if (present(dpsidz))   dpsidz   = 0
    if (present(ddpsidrr)) ddpsidrr = 0
    if (present(ddpsidzz)) ddpsidzz = 0
    if (present(ddpsidrz)) ddpsidrz = 0
  end subroutine vacuum_pol_flux


  pure subroutine vacuum_pol_curr(self, psi_n, fpol, dfpol)
    ! subroutine arguments
    class(vacuum), intent(in)            :: self
    real(wp_),     intent(in)            :: psi_n
    real(wp_),     intent(out)           :: fpol   ! poloidal current
    real(wp_),     intent(out), optional :: dfpol  ! its derivative

    unused(self); unused(psi_n)

    ! There is no current, F(ψ)=0
    fpol = 0
    if (present(dfpol)) dfpol = 0
  end subroutine vacuum_pol_curr


  pure function vacuum_safety(self, psi_n) result(q)
    ! function arguments
    class(vacuum), intent(in) :: self
    real(wp_),     intent(in) :: psi_n
    real(wp_)                 :: q

    unused(self); unused(psi_n)

    ! q(ψ) is undefined, return 0
    q = 0
  end function vacuum_safety


  pure function vacuum_conv(self, rho_in) result(rho_out)
    ! function arguments
    class(vacuum), intent(in) :: self
    real(wp_),     intent(in) :: rho_in
    real(wp_)                 :: rho_out

    unused(self)

    ! Neither ρ_p nor ρ_t are defined, do nothing
    rho_out = rho_in
  end function vacuum_conv


  pure subroutine vacuum_flux_contour(self, psi0, R_min, R, z, &
                                 R_hi, z_hi, R_lo, z_lo)
    ! subroutine arguments
    class(vacuum), intent(in)    :: self
    real(wp_),     intent(in)    :: psi0
    real(wp_),     intent(in)    :: R_min
    real(wp_),     intent(out)   :: R(:), z(:)
    real(wp_),     intent(inout) :: R_hi, z_hi, R_lo, z_lo

    unused(self); unused(psi0); unused(R_min)
    unused(R_hi); unused(z_hi); unused(R_lo); unused(z_lo)

    ! flux surfaces are undefined
    R = 0
    z = 0
  end subroutine vacuum_flux_contour


  !
  ! Helpers
  !

  subroutine load_equil(params, equil, limiter, err)
    ! Loads a generic MHD equilibrium and limiter
    ! contour from file (params%filenm)
    use gray_params, only : equilibrium_parameters, EQ_VACUUM
    use gray_params, only : EQ_ANALYTICAL, EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
    use types,       only : contour
    use logger,      only : log_warning, log_debug

    ! subroutine arguments
    type(equilibrium_parameters),       intent(inout) :: params
    class(abstract_equil), allocatable, intent(out)   :: equil
    type(contour),                      intent(out)   :: limiter
    integer,                            intent(out)   :: err

    ! local variables
    type(numeric_equil)  :: ne
    type(analytic_equil) :: ae
    type(eqdsk_data)     :: data

    select case (params%iequil)
      case (EQ_VACUUM)
        call log_debug('vacuum, doing nothing', mod='gray_equil', proc='load_equil')
        allocate(equil, source=vacuum())
        return

      case (EQ_ANALYTICAL)
        call log_debug('loading analytical file', mod='gray_equil', proc='load_equil')
        call load_analytical(params, ae, limiter, err)
        if (err /= 0) return
        allocate(equil, source=ae)

      case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
        call log_debug('loading G-EQDSK file', mod='gray_equil', proc='load_equil')
        call load_eqdsk(params, data, err)
        if (err /= 0) return

        call change_cocos(data, params%icocos, 3)
        call scale_eqdsk(params, data)

        allocate(equil, mold=ne)
        select type (equil)
          type is (numeric_equil)
            call equil%init(params, data, err)
            if (err /= 0) return
        end select

        limiter = data%limiter

    end select

  end subroutine load_equil


  subroutine load_analytical(params, equil, limiter, err)
    ! Loads the analytical equilibrium and limiter contour from file
    use const_and_precisions, only : pi, one
    use gray_params,          only : equilibrium_parameters
    use logger,               only : log_error

    ! subroutine arguments
    type(equilibrium_parameters), intent(in)  :: params
    type(analytic_equil),         intent(out) :: equil
    type(contour),                intent(out) :: limiter
    integer,                      intent(out) :: err

    ! local variables
    integer :: i, u, nlim
    real(wp_), allocatable :: R(:), z(:)

    open(file=params%filenm, status='old', action='read', newunit=u, iostat=err)
    if (err /= 0) then
      call log_error('opening equilibrium ('//trim(params%filenm)//') failed!', &
                     mod='gray_equil', proc='load_analytical')
      err = 1
      return
    end if

    ! The analytical file format is:
    !
    !  1 R0 z0 a
    !  2 B0
    !  3 q0 q1 alpha
    !  4 nlim
    !  5 R z
    !    ...

    ! Read model parameters
    read (u, *) equil%R0, equil%z0, equil%a
    read (u, *) equil%B0
    read (u, *) equil%q0, equil%q1, equil%alpha
    read (u, *) nlim

    ! Read limiter points
    if (nlim > 0) then
      allocate(R(nlim), z(nlim))
      read (u, *) (R(i), z(i), i = 1, nlim)
      limiter = contour(R, z)
    end if
    close(u)

    ! Notes on cocos=3
    ! 1. With both I_p,B_φ > 0 we have ∂ψ/∂r<0 and ∂Φ/∂r>0.
    ! 2. ψ_a = ψ_edge - ψ_axis < 0.
    ! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0.
    ! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ).

    ! Apply forced signs
    if (params%sgni /= 0) then
      equil%q0 = sign(equil%q0, -params%sgni*one)
      equil%q1 = sign(equil%q1, -params%sgni*one)
    end if
    if (params%sgnb /= 0) then
      equil%B0 = sign(equil%B0, +params%sgnb*one)
    end if

    ! Apply rescaling factors
    equil%B0 = equil%B0 * params%factb

    ! Toroidal flux at r=a:
    !
    !   Φ(a) = B₀πa² 2γ/(γ + 1)
    !
    !   where  γ=1/√(1-ε²),
    !          ε=a/R₀ is the tokamak aspect ratio
    block
      real(wp_) :: gamma
      gamma = 1/sqrt(1 - (equil%a/equil%R0)**2)
      equil%phi_a = equil%B0 * pi * equil%a**2 * 2*gamma/(gamma + 1)
    end block

    ! In cocos=3 the safety factor is
    !
    !   q(ψ) = 1/2π ∂Φ/∂ψ.
    !
    ! Given the power law of the model
    !
    !   q(ψ) = q₀ + (q₁-q₀) (ψ/ψa)^(α/2),
    !
    ! we can find ψ_a = ψ(r=a) by integrating:
    !
    !  ∫ q(ψ)dψ = 1/2π ∫ dΦ
    !  ∫₀^ψ_a q(ψ)dψ = 1/2π Φ_a
    !  ψa [q₀ + (q₁-q₀)/(α/2+1)] = Φa/2π
    !
    !  ⇒ ψ_a = Φ_a 1/2π 1/(q₀ + Δq)
    !
    ! where Δq = (q₁ - q₀)/(α/2 + 1)
    block
      real(wp_) :: dq
      dq = (equil%q1 - equil%q0) / (equil%alpha/2 + 1)
      equil%psi_a = 1/(2*pi) * equil%phi_a / (equil%q0 + dq)
    end block

    ! Compute other constants (see abstract_equil)
    equil%b_axis     = equil%B0
    equil%b_centre   = equil%B0
    equil%r_centre   = equil%R0
    equil%sgn_bphi   = sign(one, equil%B0)
    equil%axis       = [equil%R0, equil%z0]
    equil%r_range    = equil%R0 + [-equil%a, equil%a]
    equil%z_range    = equil%z0 + [-equil%a, equil%a]
    equil%z_boundary = equil%z0 + [-equil%a, equil%a]
  end subroutine load_analytical


  subroutine load_eqdsk(params, data, err)
    ! Loads the equilibrium `data` from a G-EQDSK file (params%filenm).
    ! For a description of the G-EQDSK format, see the GRAY user manual.
    use const_and_precisions, only : one
    use gray_params,          only : equilibrium_parameters
    use logger,               only : log_error

    ! subroutine arguments
    type(equilibrium_parameters), intent(in)  :: params
    type(eqdsk_data),             intent(out) :: data
    integer,                      intent(out) :: err

    ! local variables
    integer :: u, i, j, nr, nz, nlim, nbnd
    character(len=48) :: string
    real(wp_) :: r_dim, z_dim, r_left, z_mid, psi_edge, psi_axis
    real(wp_) :: skip_r ! dummy variables, used to discard data
    integer   :: skip_i !

    ! Open the G-EQDSK file
    open(file=params%filenm, status='old', action='read', newunit=u, iostat=err)
    if (err /= 0) then
      call log_error('opening eqdsk file ('//trim(params%filenm)//') failed!', &
                     mod='gray_equil', proc='load_eqdsk')
      err = 1
      return
    end if

    ! Get size of main arrays and allocate them
    if (params%idesc) then
      read (u,'(a48,3i4)') string, skip_i, nr, nz
    else
      read (u,*) nr, nz
    end if
    allocate(data%grid_r(nr), data%grid_z(nz), data%psi_map(nr, nz))
    allocate(data%psi(nr), data%fpol(nr), data%q(nr))

    ! Store 0D data and main arrays
    if (params%ifreefmt) then
      read (u, *) r_dim, z_dim, data%r_ref, r_left, z_mid
      read (u, *) data%axis, psi_axis, psi_edge, skip_r
      read (u, *) skip_r, skip_r, skip_r, skip_r, skip_r
      read (u, *) skip_r, skip_r, skip_r, skip_r, skip_r
      read (u, *) (data%fpol(i), i=1,nr)
      read (u, *) (skip_r,i=1, nr)
      read (u, *) (skip_r,i=1, nr)
      read (u, *) (skip_r,i=1, nr)
      read (u, *) ((data%psi_map(i,j), i=1,nr), j=1,nz)
      read (u, *) (data%q(i), i=1,nr)
    else
      read (u, '(5e16.9)') r_dim, z_dim, data%r_ref, r_left, z_mid
      read (u, '(5e16.9)') data%axis, psi_axis, psi_edge, skip_r
      read (u, '(5e16.9)') skip_r, skip_r, skip_r, skip_r, skip_r
      read (u, '(5e16.9)') skip_r, skip_r, skip_r, skip_r, skip_r
      read (u, '(5e16.9)') (data%fpol(i), i=1,nr)
      read (u, '(5e16.9)') (skip_r, i=1,nr)
      read (u, '(5e16.9)') (skip_r, i=1,nr)
      read (u, '(5e16.9)') (skip_r, i=1,nr)
      read (u, '(5e16.9)') ((data%psi_map(i,j), i=1,nr), j=1,nz)
      read (u, '(5e16.9)') (data%q(i), i=1,nr)
    end if

    ! Get size of boundary and limiter arrays and allocate them
    read (u, *) nbnd, nlim

    ! Load plasma boundary data
      if (nbnd > 0) then
        block
          real(wp_) :: R(nbnd), z(nbnd)
          if (params%ifreefmt) then
            read(u, *) (R(i), z(i), i=1,nbnd)
          else
            read(u, '(5e16.9)') (R(i), z(i), i=1,nbnd)
          end if
          data%boundary = contour(R, z)
        end block
      end if

    ! Load limiter data
    if (nlim > 0) then
      block
        real(wp_) :: R(nlim), z(nlim)
        if (params%ifreefmt) then
          read(u, *) (R(i), z(i), i=1,nlim)
        else
          read(u, '(5e16.9)') (R(i), z(i), i=1,nlim)
        end if
        data%limiter = contour(R, z)
      end block
    end if

    ! End of G-EQDSK file
    close(u)

    ! Build the grid arrays
    block
      real(wp_) :: dpsi, dr, dz
      dr   = r_dim/(nr - 1)
      dz   = z_dim/(nz - 1)
      dpsi = one/(nr - 1)
      do i = 1, nr
        data%psi(i)    = (i-1)*dpsi
        data%grid_r(i) = r_left + (i-1)*dr
      end do
      do i = 1, nz
        data%grid_z(i) = z_mid - z_dim/2 + (i-1)*dz
      end do
    end block

    ! Normalize the poloidal flux
    data%psi_a = psi_edge - psi_axis
    if (.not. params%ipsinorm) data%psi_map = (data%psi_map - psi_axis)/data%psi_a

  end subroutine load_eqdsk


  subroutine scale_eqdsk(params, data)
    ! Rescale the magnetic field (B) and the plasma current (I_p)
    ! and/or force their signs.
    !
    ! Notes:
    ! 1. signi and signb are ignored on input if equal to 0.
    !    They are used to assign the direction of B_φ and I_p BEFORE scaling.
    ! 2. cocos=3 assumed: positive toroidal direction is CCW from above
    ! 3. B_φ and I_p scaled by the same factor factb to keep q unchanged
    ! 4. factb<0 reverses the directions of Bphi and Ipla
    use const_and_precisions, only : one
    use gray_params,          only : equilibrium_parameters

    ! subroutine arguments
    type(equilibrium_parameters), intent(inout) :: params
    type(eqdsk_data),             intent(inout) :: data

    ! Notes on cocos=3
    ! 1. With both I_p,B_φ > 0 we have ∂ψ/∂r<0 and ∂Φ/∂r>0.
    ! 2. ψ_a = ψ_edge - ψ_axis < 0.
    ! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0.
    ! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ).

    ! Apply signs
    if (params%sgni /= 0) data%psi_a = sign(data%psi_a, -params%sgni*one)
    if (params%sgnb /= 0) data%fpol = sign(data%fpol, +params%sgnb*one)

    ! Rescale
    data%psi_a = data%psi_a * params%factb
    data%fpol  =  data%fpol * params%factb

    ! Compute the signs to be shown in the outputs header
    ! Note: this is needed if sgni, sgnb = 0 in gray.ini
    params%sgni = int(sign(one, -data%psi_a))
    params%sgnb = int(sign(one, +data%fpol(size(data%fpol))))
  end subroutine scale_eqdsk


  subroutine change_cocos(data, cocosin, cocosout, err)
    ! Convert the MHD equilibrium data from one coordinate convention
    ! (COCOS) to another. These are specified by `cocosin` and
    ! `cocosout`, respectively.
    !
    ! For more information, see: https://doi.org/10.1016/j.cpc.2012.09.010
    use const_and_precisions, only : zero, one, pi
    use logger,               only : log_warning

    ! subroutine arguments
    type(eqdsk_data), intent(inout)         :: data
    integer,          intent(in)            :: cocosin, cocosout
    integer,          intent(out), optional :: err

    ! local variables
    real(wp_) :: isign, bsign
    integer   :: exp2pi, exp2piout
    logical   :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout

    call decode_cocos(cocosin, exp2pi, phiccw, psiincr, qpos)
    call decode_cocos(cocosout, exp2piout, phiccwout, psiincrout, qposout)

    ! Check sign consistency
    isign = sign(one, data%psi_a)
    if (.not. psiincr) isign = -isign
    bsign = sign(one, data%fpol(size(data%fpol)))
    if (qpos .neqv. isign * bsign * data%q(size(data%q)) > zero) then
      ! Warning: sign inconsistency found among q, I_p and B_ref
      data%q = -data%q
      call log_warning('data not consistent with cocosin', &
                       mod='gray_equil', proc='change_cocos')
      if (present(err)) err = 1
    else
      if (present(err)) err = 0
    end if

    ! Convert cocosin to cocosout

    ! Opposite direction of toroidal angle phi in cocosin and cocosout
    if (phiccw .neqv. phiccwout) data%fpol = -data%fpol

    ! q has opposite sign for given sign of B_phi⋅I_p
    if (qpos .neqv. qposout) data%q = -data%q

    ! psi and Ip signs don't change accordingly
    if ((phiccw .eqv. phiccwout) .neqv. (psiincr .eqv. psiincrout)) &
      data%psi_a = -data%psi_a

    ! Convert Wb to Wb/rad or viceversa
    data%psi_a = data%psi_a * (2*pi)**(exp2piout - exp2pi)
  end subroutine change_cocos


  subroutine decode_cocos(cocos, exp2pi, phiccw, psiincr, qpos)
    ! Extracts the sign and units conventions from a COCOS index

    ! subroutine arguments
    integer, intent(in)  :: cocos
    integer, intent(out) :: exp2pi
    logical, intent(out) :: phiccw, psiincr, qpos

    ! local variables
    integer :: cmod10, cmod4

    cmod10 = mod(cocos, 10)
    cmod4  = mod(cmod10, 4)

    ! cocos>10 ψ in Wb, cocos<10 ψ in Wb/rad
    exp2pi = (cocos - cmod10)/10

    ! cocos mod 10 = 1,3,5,7: toroidal angle φ increasing CCW
    phiccw = (mod(cmod10, 2)== 1)

    ! cocos mod 10 = 1,2,5,6: ψ increasing with positive Ip
    psiincr = (cmod4==1 .or. cmod4==2)

    ! cocos mod 10 = 1,2,7,8: q positive for positive Bφ*Ip
    qpos = (cmod10<3 .or. cmod10>6)
  end subroutine decode_cocos

end module