src/equilibrium.f90: use the exact toroidal flux

To improve the analytical model correctness this changes the formula
of the toroidal flux to lift the large aspect ratio approximation (a << R₀).
In fact, Φ(r) = B₀πr² is technically inconsistent with the field varying
as B₀R₀/R. The exact expression is:

  Φ(r) = B₀πr² 2/[1 + √(1 - r²/R₀²)],

which is approximately equal to the former for r << R₀.

Note that this change introduces a divergence in the poloidal field at
r=R₀ (since ∂Φ/∂r → +∞), so the domain of the equilibrium has been
restricted to r<R₀, as expected.
This should not be a concern because the field outside the plasma
boundary is never directly, particularly not by the integrator.

src/equilibrium.f90

@@ -766,7 +766,7 @@ contains
 
     implicit none
 
-    real(wp_) :: dq
+    real(wp_) :: dq, gamma
 
     btaxis  = model%B0
     rmaxis  = model%R0
@@ -782,9 +782,14 @@ contains
     zmxm = zbsup
     zmnm = zbinf
 
-    ! FIXME: it should be that Φ(r) = B₀πr² 1/√(1-ε²)
+    ! Toroidal flux at r=a:
-    ! where ε=r/R₀ is the tokamak aspect ratio
+    !
-    phitedge = model%B0 * pi * model%a**2
+    !   Φ(a) = B₀πa² 2γ/(γ + 1)
+    !
+    !   where  γ=1/√(1-ε²),
+    !          ε=a/R₀ is the tokamak aspect ratio
+    gamma = 1/sqrt(1 - (model%a/model%R0)**2)
+    phitedge = model%B0 * pi * model%a**2 * 2*gamma/(gamma + 1)
 
     ! In cocos=3 the safety factor is
     !
@@ -828,7 +833,8 @@ contains
     ! derivatives wrt (R, z) up to the second order.
     !
     ! Note: all output arguments are optional.
-    use gray_params, only : iequil
+    use gray_params,          only : iequil
+    use const_and_precisions, only : one, pi
 
     implicit none
 
@@ -838,21 +844,75 @@ contains
                                         ddpsidrr, ddpsidzz, ddpsidrz
 
     ! local variables
-    real(wp_) :: rho_t, rho_p              ! √Φ_n, √ψ_n
+    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_) :: dphidr, dphidz            ! ∇Φ_n
+    real(wp_) :: phi_n                     !   Φ_n
+    real(wp_) :: dphidr, dphidz            !  ∇Φ_n
+    real(wp_) :: ddphidrdr, ddphidzdz      ! ∇∇Φ_n
+    real(wp_) :: ddphidrdz                 !
     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²)²
 
     if (iequil < 2) then
       ! Analytical model
+      !
+      ! The normalised poloidal flux ψ_n(R, z) is computed as follows:
+      !   1. ψ_n = ρ_p²
+      !   2. ρ_p = ρ_p(ρ_t), using `frhopol`, which in turns uses q(ψ)
+      !   3. ρ_t = √Φ_n
+      !   4. Φ_n = Φ(r)/Φ(a), where Φ(r) is the flux of B_φ=B₀R₀/R
+      !      through a circular surface
+      !   5. r = √[(R-R₀)²+(z-z₀)²] is the geometric minor radius
+      r_g = hypot(R - model%R0, z - model%z0)
 
-      ! ρ_t is mapped to the normalised geometrical radius,
+      ! The exact flux of the toroidal field B_φ = B₀R₀/R is:
-      ! so ρ_t = √[(R-R₀)² + (z-z₀)²]/a and ρ_p(ρ_t)
+      !
-      rho_t = hypot(R - model%R0, z - model%z0)/model%a
+      !   Φ(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 > model%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/model%R0)**2)
+      phi_n = model%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / phitedge
+      rho_t = sqrt(phi_n)
       rho_p = frhopol(rho_t)
 
+      ! For ∇Φ_n and ∇∇Φ_n we also need:
+      !
+      !   ∂Φ∂(r²)   = B₀π γ(r)
+      !   ∂²Φ∂(r²)² = B₀π γ³(r) / (2 R₀²)
+      !
+      dphidr2     = model%B0 * pi * gamma / phitedge
+      ddphidr2dr2 = model%B0 * pi * gamma**3/(2 * model%R0**2) / phitedge
+
+      ! ∇Φ_n = ∂Φ_n/∂(r²) ∇(r²)
+      ! where ∇(r²) = 2[(R-R₀), (z-z₀)]
+      dphidr = dphidr2 * 2*(R - model%R0)
+      dphidz = dphidr2 * 2*(z - model%z0)
+
+      ! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²)
+      !       = ∂²Φ_n/∂(r²)² ∇(r²)∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²)
+      ! where ∇∇(r²) = 2I
+      ddphidrdr = ddphidr2dr2 * 4*(R - model%R0)*(R - model%R0) + dphidr2*2
+      ddphidzdz = ddphidr2dr2 * 4*(z - model%z0)*(z - model%z0) + dphidr2*2
+      ddphidrdz = ddphidr2dr2 * 4*(R - model%R0)*(z - model%z0)
+
       ! ψ_n = ρ_p(ρ_t)²
       if (present(psi_n)) psi_n = rho_p**2
 
@@ -871,10 +931,6 @@ contains
       dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
       dpsidphi = (model%q0 + dq) / q
 
-      ! Since Φ_n = ρ_t²(R,z): ∇Φ_n = 2/a² [R-R₀, z-z₀]
-      dphidr = 2*(R - model%R0)/model%a**2
-      dphidz = 2*(z - model%z0)/model%a**2
-
       ! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n
       if (present(dpsidr)) dpsidr = dpsidphi * dphidr
       if (present(dpsidz)) dpsidz = dpsidphi * dphidz
@@ -890,18 +946,18 @@ contains
       !   ∇q = α/2 (q-q₀) ∇ψ_n/ψ_n
       !      = α/2 (q-q₀)/ψ_n (∂ψ_n/∂Φ_n) ∇Φ_n.
       !
-      dqdr = model%alpha/2 * (q - model%q0)/rho_p**2 * dpsidphi * dphidr
+      dqdr = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidr
-      dqdz = model%alpha/2 * (q - model%q0)/rho_p**2 * dpsidphi * dphidz
+      dqdz = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidz
       ddpsidphidr = - dpsidphi * dqdr/q
       ddpsidphidz = - dpsidphi * dqdz/q
 
-      ! Finally, ∇∇Φ_n = ∇∇ ρ_t²(R,z) = 2/a² I, so:
+      ! Combining all of the above:
       !
-      !   ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) 2/a² I
+      !   ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n
       !
-      if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * 2/model%a**2
+      if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * ddphidrdr
-      if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * 2/model%a**2
+      if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * ddphidzdz
-      if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz
+      if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz
     else
       ! Numerical data
       if (R <= rmxm .and. R >= rmnm .and. &