diff --git a/src/equilibrium.f90 b/src/equilibrium.f90
index 9d6de83..b74bf69 100644
--- a/src/equilibrium.f90
+++ b/src/equilibrium.f90
@@ -782,7 +782,27 @@ contains
     zmxm = zbsup
     zmnm = zbinf
 
+    ! FIXME: it should be that Φ(r) = B₀πr² 1/√(1-ε²)
+    ! where ε=r/R₀ is the tokamak aspect ratio
     phitedge = model%B0 * pi * model%a**2
+
+    ! 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)
     dq   = (model%q1 - model%q0) / (model%alpha/2 + 1)
     psia = 1/(2*pi) * phitedge / (model%q0 + dq)
   end subroutine set_equil_an
@@ -818,26 +838,70 @@ contains
                                         ddpsidrr, ddpsidzz, ddpsidrz
 
     ! local variables
-    real(wp_) :: r_p   ! poloidal radius
-    real(wp_) :: rho_p ! poloidal radius normalised
+    real(wp_) :: rho_t, rho_p              ! √Φ_n, √ψ_n
+    real(wp_) :: dpsidphi                  !  (∂ψ_n/∂Φ_n)
+    real(wp_) :: ddpsidphidr, ddpsidphidz  ! ∇(∂ψ_n/∂Φ_n)
+    real(wp_) :: dphidr, dphidz            ! ∇Φ_n
+    real(wp_) :: q, dq                     ! q(ρ_p), Δq=(q₁-q₀)/(α/2 + 1)
+    real(wp_) :: dqdr, dqdz                ! ∇q
 
     if (iequil < 2) then
       ! Analytical model
-      ! ρ_p is just the geometrical poloidal radius divided by a
-      r_p = hypot(R - model%R0, z - model%z0)
-      rho_p = r_p / model%a
 
-      ! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a²
+      ! ρ_t is mapped to the normalised geometrical radius,
+      ! so ρ_t = √[(R-R₀)² + (z-z₀)²]/a and ρ_p(ρ_t)
+      rho_t = hypot(R - model%R0, z - model%z0)/model%a
+      rho_p = frhopol(rho_t)
+
+      ! ψ_n = ρ_p(ρ_t)²
       if (present(psi_n)) psi_n = rho_p**2
 
-      ! ∂ψ_n/∂R, ∂ψ_n/∂z
-      if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2
-      if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**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  = model%q0 + (model%q1 - model%q0) * rho_p**model%alpha
+      dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
+      dpsidphi = (model%q0 + dq) / q
 
-      ! ∂²ψ_n/∂R², ∂²ψ_n/∂z², ∂²ψ_n/∂R∂z
-      if (present(ddpsidrr)) ddpsidrr = 2/model%a**2
-      if (present(ddpsidzz)) ddpsidzz = 2/model%a**2
-      if (present(ddpsidrz)) ddpsidrz = 0
+      ! 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
+
+      ! 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 = model%alpha/2 * (q - model%q0)/rho_p**2 * dpsidphi * dphidr
+      dqdz = model%alpha/2 * (q - model%q0)/rho_p**2 * dpsidphi * dphidz
+      ddpsidphidr = - dpsidphi * dqdr/q
+      ddpsidphidz = - dpsidphi * dqdz/q
+
+      ! Finally, ∇∇Φ_n = ∇∇ ρ_t²(R,z) = 2/a² I, so:
+      !
+      !   ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) 2/a² I
+      !
+      if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * 2/model%a**2
+      if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * 2/model%a**2
+      if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz
     else
       ! Numerical data
       if (R <= rmxm .and. R >= rmnm .and. &