rework the analytical model

This change modifies the analytical equilibrium in order to simplify the
computation of the poloidal flux normalization and the derivatives.

In the power law parametrisation of the safety factor, ρ_t is replaced
with ρ_p and, similarly, the normalised poloidal radius is now
identified with ρ_p, instead of ρ_t.
With the same parameters (q₀,q₁,α...), this choice slightly changes the
plasma current distribution, but enables us to obtain a closed form for
ψ_a = ψ(r=a) and the relation ρ_t(ρ_p). In fact, both expressions are
now obtained by integrating the q(ρ_p), instead of 1/q(ρ_t), which has
no elementary antiderivative.

As the normalisation is now computed exactly, the values of the
normalised flux ψ_n = ψ/ψ_a and the gradient ∇ψ (entering the raytracing
equations in X and ∇X, respectively) are computed to the same precision.
Previously, ψ_n was computed to a lower precision due to the use of a
simple trapezoid integration of 1/q(ρ_p) for ψ_a, while ∇ψ was computed
up to machine precision using an exact formula.

This error effectively caused a very slight decoupling between X=ω_p²/ω²
and ∇X that introduced a systematic error in the numerical solution of
the raytracing equations.
The error manifests itself as a bias with a weak dependency on X in the
values taken by the dispersion function Λ(r̅, n̅) on the phase-space
points generated by the integrator. More specifically,

  lim h→0 Λ(r̅_i, n̅_i) = -kX(r̅_i)

where h is the integrator step size;
      r̅_i is the position at the i-th step;
      k ≈ -3.258⋅10⁻⁵ and depends only on the number of points used to
        perform the trapedoid integral for ψ_a (as ~ 1/n²).

After this change Λ behaves consistently with being a conserved quantity
(zero) up to the cumulative integration error of the 4° order
Runge-Kutta method. In fact we now have that:

  Λ(r̅_i, n̅_i) ∝ - h⁴ ‖∂⁴X(r̅_i)/∂r̅⁴‖

It must be said that within this model the relation ρ_p(ρ_t) can't be
computed analytically (inverting ρ_t(ρ_p) produces a trascendental
equation of the form b = x + c x^α). However, this relation is not
necessary for raytracing and is easily solved, up to machine
precision, using minpack.

In addition, this change also makes the model consistetly use the
cocos=3 and fully implements the ability to force the signs of I_p, B_φ
(via equilibrium.sgni,sgnb) and rescaling the field (via
equilibrium.factb).

input/equil_a.txt

@@ -1,4 +1,8 @@
- : rr0m,zr0m,rpam  ! rhot[m] = min(sqrt((r-rr0m)**2+(z-zr0m)**2), rpam); tor flux phi = pi*b0*rhot**2
+2.96 0.0 1.25  ! R₀,z₀,a,  where ρ_p(R, z) = √[(R - R₀)² + (z - z₀)²]/a  (m)
- : b0              ! Bphi[T] @ rr0m[m]
+6              ! B₀,       where B_φ(R) = B₀ R₀/R  (T)
- : q0,qa,alq       ! q = q0 + (qa-q0)*rhot**alq
+3.5 10 2       ! q₀,q₁,α,  where q(ρ_p) = q₀ + (q₁-q₀)ρ_p^α
- : nlim            ! number of points in first wall (rlim,zlim) polygon
+0              ! number of points in the limiter contourn, (R,z) pairs
+
+! Notes:
+! 1. use B₀>0 for clockwise B_φ
+! 2. use q₀,q₁<0 for clockwise I_p,B_φ

input/profil_a.txt

@@ -1,3 +1,3 @@
-: dens0,aln1,aln2      ! dens[1e19 m-3] = dens0*(1 - psin**aln1)**aln2
+10 0.3 3  ! n₀,a,b, where n(ψ) = n₀(1 - ψ^a)^b (10¹⁹ m⁻³)
-: te0,dte0,alt1,alt2   ! temp[keV] = (te0-dte0)*(1 - psin**alt1)**alt2 + dte0
+14 0 2 8  ! T₀,T₁,a,b, where T(ψ) = (T₀ - T₁)(1 - ψ^a)^b + T₁  (keV)
-: zeffan
+1.0       ! Z_eff

src/equilibrium.f90

@@ -15,7 +15,11 @@ module equilibrium
   type analytic_model
     real(wp_) :: q0    ! Safety factor at the magnetic axis
     real(wp_) :: q1    ! Safety factor at the edge
-    real(wp_) :: lq ! Exponent for the q(ρ) power law
+    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)
   end type
 
   ! Order of the splines
@@ -33,13 +37,18 @@ module equilibrium
   type(analytic_model), save :: model
 
   ! More parameters
-  real(wp_), save :: fpolas
+  real(wp_), save :: fpolas          ! Poloidal current at the edge (r=a), F_a
-  real(wp_), save :: psia, psiant, psinop
+  real(wp_), save :: psia            ! Poloidal flux at the edge (r=a), ψ_a
-  real(wp_), save :: phitedge, aminor
+  real(wp_), save :: psiant          ! Normalised poloidal flux at the true edge
-  real(wp_), save :: btaxis,rmaxis,zmaxis,sgnbphi
+  real(wp_), save :: psinop          ! Normalised poloidal flux at the true O point
+  real(wp_), save :: phitedge        ! Toroidal flux at the edge
+  real(wp_), save :: btaxis          ! B₀: B_φ = B₀R₀/R
+  real(wp_), save :: rmaxis, zmaxis  ! R,z of the magnetic axis (O point)
+  real(wp_), save :: sgnbphi         ! Sign of B_φ (>0 means CCW)
   real(wp_), save :: btrcen          ! used only for jcd_astra def.
-  real(wp_), save :: rcen                          ! computed as fpol(a)/btrcen
+  real(wp_), save :: rcen            ! Major radius R₀, computed as fpolas/btrcen
-  real(wp_), save :: rmnm, rmxm, zmnm, zmxm
+  real(wp_), save :: rmnm, rmxm      ! R interval of the equilibrium definition
+  real(wp_), save :: zmnm, zmxm      ! z interval of the equilibrium definition
   real(wp_), save :: zbinf, zbsup
 
   private
@@ -48,16 +57,16 @@ module equilibrium
   public equian                              ! Analytical model
   public equinum_psi, equinum_fpol           ! Interpolated data
   public fq, bfield, tor_curr, tor_curr_psi  ! Accessing local quantities
-  public frhopol, frhopolv, frhotor          ! Converting between poloidal/toroidal flux
+  public frhotor, frhopol                    ! Converting between poloidal/toroidal flux
   public set_equil_spline, set_equil_an      ! Initialising internal state
   public unset_equil_spline                  ! Deinitialising internal state
 
   ! Members exposed for magsurf_data
-  public kspl, psi_spline, q_spline, points_tgo
+  public kspl, psi_spline, q_spline, points_tgo, model
 
   ! Members exposed to gray_core and more
   public psia, psiant, psinop
-  public phitedge, aminor
+  public phitedge
   public btaxis, rmaxis, zmaxis, sgnbphi
   public btrcen, rcen
   public rmnm, rmxm, zmnm, zmxm
@@ -213,7 +222,6 @@ contains
 
     ! local variables
     integer :: i, u, nlim
-    real(wp_) :: rr0m, zr0m, rpam, b0
 
     u = get_free_unit(unit)
 
@@ -225,20 +233,9 @@ contains
       return
     end if
 
-    read(u, *) rr0m, zr0m, rpam
+    read(u, *) model%R0, model%z0, model%a
-    read(u, *) b0
+    read(u, *) model%B0
-    read(u, *) model%q0, model%q1, model%lq
+    read(u, *) model%q0, model%q1, model%alpha
-
-    if(allocated(data%rv))   deallocate(data%rv)
-    if(allocated(data%zv))   deallocate(data%zv)
-    if(allocated(data%fpol)) deallocate(data%fpol)
-    if(allocated(data%qpsi)) deallocate(data%qpsi)
-    allocate(data%rv(2), data%zv(1), data%fpol(1), data%qpsi(3))
-
-    data%rv   = [rr0m, rpam]
-    data%zv   = [zr0m]
-    data%fpol = [b0*rr0m]
-    data%qpsi = [model%q0, model%q1, model%lq]
 
     if(ipass >= 2) then
       ! get size of boundary and limiter arrays and allocate them
@@ -252,7 +249,6 @@ contains
         read(u,*) (data%rlim(i), data%zlim(i), i = 1, nlim)
       end if
     end if
-
     close(u)
   end subroutine read_equil_an
 
@@ -340,18 +336,19 @@ contains
 
 
   subroutine scale_equil(params, data)
-    ! Rescale the magnetic field (B) and the plasma current (I)
+    ! 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 Bphi and Ipla BEFORE scaling.
+    !    They are used to assign the direction of B_φ and I_p BEFORE scaling.
-    ! 2. cocos=3 assumed: CCW direction is >0
+    ! 2. cocos=3 assumed: positive toroidal direction is CCW from above
-    ! 3. Bphi and Ipla scaled by the same factor factb to keep q unchanged
+    ! 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, equilibrium_data
+    use gray_params,          only : iequil
 
     implicit none
 
@@ -359,20 +356,42 @@ contains
     type(equilibrium_parameters), intent(inout) :: params
     type(equilibrium_data),       intent(inout) :: data
 
-    ! local variables
+    ! Notes on cocos=3
-    real(wp_) :: last_fpol
+    ! 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_φ).
 
-    last_fpol = data%fpol(size(data%fpol))
+    if (iequil < 2) then
+      ! Analytical model
 
+      ! Apply signs
+      if (params%sgni /= 0) then
+        model%q0 = sign(model%q0, -params%sgni*one)
+        model%q1 = sign(model%q1, -params%sgni*one)
+      end if
+      if (params%sgnb /= 0) then
+        model%B0 = sign(model%B0, +params%sgnb*one)
+      end if
+      ! Rescale
+      model%B0 = model%B0 * params%factb
+    else
+      ! Numeric data
+
+      ! Apply signs
       if (params%sgni /= 0) &
-      data%psia = sign(data%psia, real(-params%sgni, wp_))
+        data%psia = sign(data%psia, -params%sgni*one)
-    if (params%sgnb /=0 .and. params%sgnb * last_fpol < 0) &
+      if (params%sgnb /= 0) &
-      data%fpol = -data%fpol
+        data%fpol = sign(data%fpol, +params%sgnb*one)
-
+      ! Rescale
       data%psia = data%psia * params%factb
       data%fpol = data%fpol * params%factb
+
+      ! Compute the signs to be shown in the outputs header when cocos≠0,10.
+      ! Note: In these cases the values sgni,sgnb from gray.ini are unused.
       params%sgni = int(sign(one, -data%psia))
-    params%sgnb = int(sign(one, last_fpol))
+      params%sgnb = int(sign(one, +data%fpol(size(data%fpol))))
+    end if
   end subroutine scale_equil
 
 
@@ -713,7 +732,7 @@ contains
 
     call q_spline%init(psinq, q, size(q))
 
-    ! Toroidal flux φ = 2π ∫q(ψ)dψ
+    ! Toroidal flux as Φ(ψ) = 2π ∫q(ψ)dψ
     phit(1)=0
     do k=1,q_spline%ndata-1
       dx=q_spline%data(k+1)-q_spline%data(k)
@@ -726,83 +745,32 @@ contains
   end subroutine setqphi_num
 
 
-  subroutine set_equil_an(data, n)
+  subroutine set_equil_an
     ! Computes the analytical equilibrium data and stores them
     ! in their respective global variables, see the top of this file.
-    use const_and_precisions, only : pi, zero, one
+    use const_and_precisions, only : pi, one
-    use gray_params,          only : equilibrium_data
 
     implicit none
 
-    ! subroutine arguments
+    real(wp_) :: dq
-    type(equilibrium_data), intent(in) :: data
-    integer, optional,      intent(in) :: n
 
-    ! local variables
+    btaxis  = model%B0
-    integer, parameter :: nqdef=101
+    rmaxis  = model%R0
-    integer :: i
+    zmaxis  = model%z0
-    real(wp_) :: dr, fq0, fq1, qq, res, rn
+    btrcen  = model%B0
-    real(wp_) :: rax, zax, a, bax, qax, q1, qexp
+    rcen    = model%R0
-    real(wp_), dimension(:), allocatable :: rhotn,rhopn
+    zbinf   = model%z0 - model%a
+    zbsup   = model%z0 + model%a
+    sgnbphi = sign(one, model%B0)
 
-    rax  = data%rv(1)
+    rmxm = model%R0 + model%a
-    zax  = data%zv(1)
+    rmnm = model%R0 - model%a
-    a    = data%rv(2)
-    bax  = data%fpol(1) / rax
-    qax  = data%qpsi(1)
-    q1   = data%qpsi(2)
-    qexp = data%qpsi(3)
-
-    btaxis = bax
-    rmaxis = rax
-    zmaxis = zax
-    btrcen = bax
-    rcen   = rax
-    aminor = a
-    zbinf = zmaxis-a
-    zbsup = zmaxis+a
-    model%q0 = qax
-    model%q1 = q1
-    model%lq = qexp
-    sgnbphi = sign(one, bax)
-
-    rmxm = rmaxis+aminor
-    rmnm = rmaxis-aminor
     zmxm = zbsup
     zmnm = zbinf
 
-    if (present(n)) then
+    phitedge = model%B0 * pi * model%a**2
-      q_spline%ndata = n
+    dq   = (model%q1 - model%q0) / (model%alpha/2 + 1)
-    else
+    psia = 1/(2*pi) * phitedge / (model%q0 + dq)
-      q_spline%ndata = nqdef
-    end if
-
-    if (allocated(q_spline%data)) deallocate(q_spline%data)
-    allocate(q_spline%data(q_spline%ndata))
-    allocate(rhotn(q_spline%ndata), rhopn(q_spline%ndata))
-
-    dr = one/(q_spline%ndata - 1)
-    rhotn(1) = zero
-    q_spline%data(1) = zero
-    res = zero
-    fq0 = zero
-    do i = 2, q_spline%ndata
-      rn = (i - 1)*dr
-      qq  = model%q0 + (model%q1 - model%q0) * rn**model%lq
-      fq1 = rn / qq
-      res = res + (fq1 + fq0)/2 * dr
-      fq0 = fq1
-      rhotn(i) = rn
-      q_spline%data(i) = res
-    end do
-    
-    phitedge = btaxis*aminor**2 ! temporary
-    psia     = res*phitedge
-    phitedge = pi*phitedge    ! final
-    q_spline%data = q_spline%data/res
-    rhopn = sqrt(q_spline%data)
-
-    call set_rho_spline(rhopn, rhotn)
   end subroutine set_equil_an
 
 
@@ -874,120 +842,142 @@ contains
   end subroutine equinum_fpol
 
 
-  subroutine equian(R, z, psi, fpolv, dfpv, dpsidr, dpsidz,  &
+  subroutine equian(R, z, psin, fpolv, dfpv, dpsidr, dpsidz,  &
                     ddpsidrr, ddpsidzz, ddpsidrz)
     ! Computes all MHD equilibrium parameters from a simple analytical model
     implicit none
 
     ! subroutine arguments
     real(wp_), intent(in) :: R, z
-    real(wp_),  intent(out), optional :: psi, fpolv, dfpv, dpsidr, dpsidz,  &
+    real(wp_), intent(out), optional :: psin, fpolv, dfpv, dpsidr, dpsidz,  &
                                         ddpsidrr, ddpsidzz, ddpsidrz
 
-    ! local variables
+    ! variables
-    real(wp_) :: cst, dpsidrp, d2psidrp, dqq, qq, &
+    real(wp_) :: r_p   ! poloidal radius
-                 rn, rpm, snt, rhop, rhot, btaxqq
+    real(wp_) :: rho_p ! poloidal radius normalised
 
-    ! simple model for equilibrium: large aspect ratio
+    ! ρ_p is just the geometrical poloidal radius divided by a
-    ! outside plasma: analytical continuation, not solution Maxwell eqs
+    r_p = hypot(R - model%R0, z - model%z0)
+    rho_p = r_p / model%a
 
-    ! Note: rpm is ρ_t in meters, rn is ρ_t normalised
+    ! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a²
-    rpm = hypot(R - rmaxis, z - zmaxis)
+    if (present(psin)) psin = rho_p**2
-    rn  = rpm/aminor
 
-    snt = 0
+    ! ∂ψ/∂R, ∂ψ/∂z, note that ψ = ψ_n⋅ψ_a
-    cst = 1
+    if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2 * psia
-    if (rpm > 0) then
+    if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**2 * psia
-      snt = (z - zmaxis)/rpm
-      cst = (R - rmaxis)/rpm
-    end if
 
-    if (present(psi)) then
+    ! ∂²ψ/∂R², ∂²ψ/∂z², ∂²ψ/∂R∂z
-      rhot=rn
+    if (present(ddpsidrr)) ddpsidrr = 2/model%a**2 * psia
-      if(rn <= 1) then
+    if (present(ddpsidzz)) ddpsidzz = 2/model%a**2 * psia
-        rhop = frhopol(rhot)
+    if (present(ddpsidrz)) ddpsidrz = 0
-        psi = rhop**2
-      else
-        psi = 1 + btaxis / (2*psia*model%q1) * (rpm**2 - aminor**2)
-        rhop = sqrt(psi)
-      end if
-    end if
 
-    if(rn <= 1) then
+    ! F(ψ) = B₀⋅R₀, a constant
-      qq = model%q0 + (model%q1 - model%q0) * rn**model%lq
+    if (present(fpolv)) fpolv = model%B0 * model%R0
-      btaxqq = btaxis/qq
-      dpsidrp = btaxqq*rpm
-      dqq = model%lq * (model%q1 - model%q0) * rn**(model%lq - 1)
-      d2psidrp=btaxqq*(1 - rn*dqq/qq)
-    else
-      btaxqq = btaxis / model%q1
-      dpsidrp = btaxqq * rpm
-      d2psidrp = btaxqq
-    end if
-
-    if(present(fpolv)) fpolv = btaxis * rmaxis
     if (present(dfpv))   dfpv = 0
-
-    if(present(dpsidr))   dpsidr = dpsidrp*cst
-    if(present(dpsidz))   dpsidz = dpsidrp*snt
-    if(present(ddpsidrr)) ddpsidrr = btaxqq*snt**2 + cst**2*d2psidrp
-    if(present(ddpsidrz)) ddpsidrz = cst * snt * (d2psidrp - btaxqq)
-    if(present(ddpsidzz)) ddpsidzz = btaxqq * cst**2 + snt**2 * d2psidrp
   end subroutine equian
 
 
-  function frhopol(rhot)
+  function frhotor(rho_p)
-    ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius
-
-    implicit none
-
-    ! function arguments
-    real(wp_), intent(in) :: rhot
-    real(wp_)             :: frhopol
-
-    frhopol = rhop_spline%eval(rhot)
-  end function frhopol
-
-
-  function frhopolv(rhot)
-    ! Vector variant of `frhopol`
-    use utils, only : locate
-
-    implicit none
-
-    ! function arguments
-    real(wp_), intent(in) :: rhot(:)
-    real(wp_)             :: frhopolv(size(rhot))
-
-    ! local variables
-    integer :: i, i0, j
-
-    i0 = 1
-    do j = 1, size(rhot)
-      call locate(rhop_spline%xdata(i0:), rhop_spline%ndata-i0+1, rhot(j), i)
-      i = min(max(1,i), rhop_spline%ndata - i0) + i0 - 1
-      frhopolv(j) = rhop_spline%raw_eval(i, rhot(j))
-      i0 = i
-    end do
-  end function frhopolv
-
-
-  function frhotor(rhop)
     ! Converts from poloidal (ρ_p) to toroidal (ρ_t) normalised radius
+    use gray_params, only : iequil
 
     implicit none
 
     ! function arguments
-    real(wp_), intent(in) :: rhop
+    real(wp_), intent(in) :: rho_p
     real(wp_)             :: frhotor
 
-    frhotor = rhot_spline%eval(rhop)
+    if (iequil < 2) then
+      ! Analytical model
+      block
+        ! 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)
+        real(wp_) :: dq
+        dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
+        frhotor = rho_p * sqrt((model%q0 + dq*rho_p**model%alpha) &
+                               / (model%q0 + dq))
+      end block
+    else
+      ! Numerical data
+      frhotor = rhot_spline%eval(rho_p)
+    end if
   end function frhotor
 
 
+  function frhopol(rho_t)
+    ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius
+    use gray_params,          only : iequil
+    use const_and_precisions, only : comp_eps
+
+    implicit none
+
+    ! function arguments
+    real(wp_), intent(in) :: rho_t
+    real(wp_)             :: frhopol
+
+    if (iequil < 2) then
+      ! Analytical model
+      block
+        ! 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.
+        use minpack, only : hybrj1
+
+        real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7)
+        integer   :: info
+
+        rho_p = [rho_t]  ! 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)
+        frhopol = rho_p(1)
+      end block
+    else
+      ! Numerical data
+      frhopol = rhop_spline%eval(rho_t)
+    end if
+
+    contains
+
+      subroutine equation(n, x, f, df, ldf, flag)
+        ! The equation to solve: f(x) = ρ_t(x) - ρ_t₀ = 0
+        implicit none
+
+        ! 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) = frhotor(x(1)) - rho_t
+        else
+          ! return f'(x), computed numerically
+          df(1,1) = (frhotor(x(1) + e) - frhotor(x(1) - e)) / (2*e)
+        end if
+      end subroutine
+
+  end function frhopol
+
+
   function fq(psin)
-    ! Computes the safety factor q as a function of the
+    ! Computes the safety factor q as a function of the normalised
-    ! normalised poloidal flux ψ
+    ! poloidal flux ψ.
+    !
+    ! Note: this returns the absolute value of q.
     use gray_params, only : iequil
 
     implicit none
@@ -996,14 +986,17 @@ contains
     real(wp_), intent(in) :: psin
     real(wp_) :: fq
 
-    ! local variables
-    real(wp_) :: rn
-
     if (iequil < 2) then
-      ! q = q₀ + (q₁ - q₀)ρ^l
+      ! Analytical model
-      rn = frhotor(sqrt(psin))
+      ! The safety factor is a power law in ρ_p:
-      fq = model%q0 + (model%q1 - model%q0) * rn**model%lq
+      !   q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α
+      block
+        real(wp_) :: rho_p
+        rho_p = sqrt(psin)
+        fq = abs(model%q0 + (model%q1 - model%q0) * rho_p**model%alpha)
+      end block
     else
+      ! Numerical data
       fq = q_spline%eval(psin)
     end if
   end function fq
@@ -1024,9 +1017,11 @@ contains
     real(wp_) :: psin,fpol
 
     if (iequil < 2) then
+      ! Analytical model
       call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br)
       if (present(bphi)) bphi=bphi/rpsim
     else
+      ! Numerical data
       call equinum_psi(rpsim,zpsim,psi=bphi,dpsidr=bz,dpsidz=br)
       if (present(bphi)) then
         psin=bphi
@@ -1034,6 +1029,7 @@ contains
         bphi=fpol/rpsim
       end if
     end if
+
     if (present(br)) br=-br/rpsim
     if (present(bz)) bz= bz/rpsim
   end subroutine bfield
@@ -1055,9 +1051,11 @@ contains
     real(wp_) :: dpsidr,ddpsidrr,ddpsidzz
 
     if(iequil < 2) then
+      ! Analytical model
       call equian(rpsim,zpsim,dpsidr=dpsidr, &
                   ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
     else
+      ! Numerical data
       call equinum_psi(rpsim,zpsim,dpsidr=dpsidr, &
                        ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
     end if
@@ -1105,10 +1103,12 @@ contains
     character(64) :: msg
 
     if (iequil < 2) then
-      val=frhotor(sqrt(psin))
+      ! Analytical model
-      r1=rmaxis-val*aminor
+      val = sqrt(psin)
-      r2=rmaxis+val*aminor
+      r1 = model%R0 - val*model%a
+      r2 = model%R0 + val*model%a
     else
+      ! Numerical data
       iopt=1
       zc=zmaxis
       call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, &

src/gray_core.f90

@@ -2328,7 +2328,7 @@ bb: do
       do k = 1, q_spline%ndata
         zk = z(k)
         if (iequil < 2) then
-          call equian(rj, zk, psi=psin, fpolv=bphi, dpsidr=bz, dpsidz=br)
+          call equian(rj, zk, psin=psin, fpolv=bphi, dpsidr=bz, dpsidz=br)
         else
           call equinum_psi(rj, zk, psi=psin, dpsidr=bz, dpsidz=br)
           call equinum_fpol(psin, fpol=bphi)

src/magsurf_data.f90

@@ -447,7 +447,7 @@ contains
     use gray_params, only : iequil
     use logger,      only : log_warning
     use dierckx,     only : profil, sproota
-    use equilibrium, only : rmaxis, zmaxis, aminor, frhotor, psi_spline, &
+    use equilibrium, only : model, frhotor, psi_spline, &
                             kspl, psiant, psinop, points_tgo
     use limiter,     only : rwallm
 
@@ -463,7 +463,7 @@ contains
 
     ! local variables
     integer :: npoints,np,info,ic,ier,iopt,m
-    real(wp_) :: ra,rb,za,zb,rn,th,zc,val
+    real(wp_) :: ra,rb,za,zb,th,zc,val
     real(wp_), dimension(mest) :: zeroc
     real(wp_), dimension(psi_spline%nknots_x) :: czc
 
@@ -472,13 +472,14 @@ contains
     th=pi/dble(np)
 
     if (iequil<2) then
-
+      block
-      rn=frhotor(sqrt(h))
+        real(wp_) :: r_p  ! poloidal radius
+        r_p = sqrt(h) * model%a
         do ic=1,npoints
-        zcn(ic)=zmaxis+aminor*rn*sin(th*(ic-1))
+          rcn(ic) = model%R0 + r_p * cos(th*(ic-1))
-        rcn(ic)=rmaxis+aminor*rn*cos(th*(ic-1))
+          zcn(ic) = model%z0 + r_p * sin(th*(ic-1))
         end do
-
+      end block
     else
 
       ra=rup

src/main.f90

@@ -1,5 +1,5 @@
 program main
-  use const_and_precisions, only : wp_, one, zero
+  use const_and_precisions, only : wp_, zero
   use logger,               only : INFO, ERROR, set_log_level, log_message
   use units,                only : set_active_units, close_units
   use utils,                only : dirname
@@ -226,11 +226,6 @@ contains
                          params%raytracing%ipass,   &
                          data%equilibrium, err)
       if (err /= 0) return
-
-      ! Set psia sign to give the correct sign to Iphi
-      ! (COCOS=3: psia<0 for Iphi>0)
-      data%equilibrium%psia = sign(one, data%equilibrium%qpsi(2) &
-                                      * data%equilibrium%fpol(1))
     else
       ! Numerical equilibrium
       call log_debug('loading G-EQDK file', &
@@ -245,7 +240,7 @@ contains
 
     ! Set global variables (for splines)
     if (params%equilibrium%iequil < 2) then
-      call set_equil_an(data%equilibrium)
+      call set_equil_an
     else
       call log_debug('computing splines...', mod='main', proc='init_equilibrium')
       call set_equil_spline(params%equilibrium, data%equilibrium, err)
@@ -281,7 +276,7 @@ contains
     ! temperature, density and plasma effective charge) and initialises
     ! the respective GRAY data structure.
     use gray_params,  only : profiles_parameters, profiles_data
-    use equilibrium,  only : frhopolv
+    use equilibrium,  only : frhopol
     use coreprofiles, only : read_profiles_an, read_profiles, &
                              scale_profiles, set_profiles_an, &
                              set_profiles_spline
@@ -296,6 +291,9 @@ contains
     type(profiles_data),       intent(out)   :: data
     integer,                   intent(out)   :: err
 
+    ! local variables
+    integer :: i
+
     if (params%iprof == 0) then
       ! Analytical profiles
       call log_debug('loading analytical file',  &
@@ -312,7 +310,7 @@ contains
       ! Convert psrad to ψ
       select case (params%irho)
         case (0) ! psrad is ρ_t = √Φ (toroidal flux)
-          data%psrad = frhopolv(data%psrad)**2
+          data%psrad = [(frhopol(data%psrad(i))**2, i=1,size(data%psrad))]
         case (1) ! psrad is ρ_p = √ψ (poloidal flux)
           data%psrad = data%psrad**2
         case default ! psrad is already ψ
@@ -420,7 +418,10 @@ contains
         ! Max radius, either due to the plasma extent or equilibrium grid
         if (params%equilibrium%iequil < 2) then
           ! analytic equilibrium, use R₀+a
-          R_max = data%equilibrium%rv(1) + data%equilibrium%rv(2)
+          block
+            use equilibrium, only : model
+            R_max = model%R0 + model%a
+          end block
         else
           ! numeric equilibrium, use max R of the grid
           R_max = data%equilibrium%rv(size(data%equilibrium%rv))