diff --git a/src/gray_equil.f90 b/src/gray_equil.f90
index 886b6cc..2ad4741 100644
--- a/src/gray_equil.f90
+++ b/src/gray_equil.f90
@@ -81,7 +81,7 @@ module gray_equil
                                 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,z are the contour points ordered clockwise
       !  - 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
@@ -435,7 +435,7 @@ contains
     !   q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α
     !
     rho_p = sqrt(psi_n)
-    q = abs(self%q0 + (self%q1 - self%q0) * rho_p**self%alpha)
+    q = self%q0 + (self%q1 - self%q0) * rho_p**self%alpha
   end function analytic_safety
 
 
@@ -541,8 +541,8 @@ contains
     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))
+      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
@@ -955,9 +955,10 @@ contains
       real(wp_) :: dx
       integer   :: k
 
-      call self%q_spline%init(data%psi, abs(data%q), npsi)
+      call self%q_spline%init(data%psi, data%q, npsi)
 
-      ! Toroidal flux as Φ(ψ) = 2π ∫q(ψ)dψ
+      ! Toroidal flux as Φ(ψ) = 2π ∫ q(ψ)dψ
+      ! Specifically, Φ(ψ_n) = 2π|ψ_a| ∫₀ψ_n q(ψ_n)dψ_n, with ψ_n ∈ [0,1]
       phi(1) = 0
       do k = 1, npsi-1
         dx = self%q_spline%data(k+1) - self%q_spline%data(k)
diff --git a/src/gray_tables.f90 b/src/gray_tables.f90
index d0de7ce..ec77af8 100644
--- a/src/gray_tables.f90
+++ b/src/gray_tables.f90
@@ -366,10 +366,10 @@ contains
 
       if (flag == 1) then
         ! return f(x)
-        f(1) = equil%safety(x(1)) - qvals(i)
+        f(1) = abs(equil%safety(x(1))) - qvals(i)
       else
         ! return f'(x), computed numerically
-        df(1,1) = (equil%safety(x(1) + e) - equil%safety(x(1) - e)) / (2*e)
+        df(1,1) = (abs(equil%safety(x(1) + e)) - abs(equil%safety(x(1) - e))) / (2*e)
       end if
     end subroutine
 
diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90
index 849b02b..2256b31 100644
--- a/src/magsurf_data.f90
+++ b/src/magsurf_data.f90
@@ -108,7 +108,7 @@ contains
     block
       real(wp_) :: a, b
       call equil%pol_flux(equil%axis(1), equil%axis(2), ddpsidrr=a, ddpsidzz=b)
-      q(1) = abs(equil%b_axis) / sqrt(a*b)/equil%psi_a
+      q(1) = equil%b_axis / sqrt(a*b)/equil%psi_a
     end block
 
     ! Compute the uniform λ grid and H(0, λ)
@@ -140,8 +140,8 @@ contains
 
       block
         ! Variables for the previous and current contour point
-        real(wp_) :: B_tot0, B_pol0, B_tot1, B_pol1, J_phi0, J_phi1
-        real(wp_) :: dl  ! line element
+        real(wp_) :: B_tot0, B_tot1, B_pol0, B_pol1, J_phi0, J_phi1
+        real(wp_) :: dl, dir  ! line element, sgn(dl⋅B_p)
 
         ! Initialise all integrals
         area(i)      = 0
@@ -184,6 +184,10 @@ contains
 
             ! The line differential is difference of the two vertices
             dl = norm2(P - Q)
+
+            ! We always assume dl ∥ B_p (numerically they're not),
+            ! but compute the sign of dl⋅B_p for the case I_p < 0.
+            dir = sign(one, dot_product([B_R, B_z], P - Q))
           end block
 
           ! Evaluate integrands at the current point for line integrals ∮ dl
@@ -197,18 +201,18 @@ contains
           ! Do one step of the trapezoid method
           ! N     = ∮ dl/B_p
           ! dA/dψ = ∮ (1/R) dl/B_p
-          ! I_p   = 1/μ₀ ∫ B_p⋅dS_φ = 1/μ₀ ∮ B_p dl
           ! ⟨B⟩   = 1/N ∮ B  dl/B_p
           ! ⟨B²⟩  = 1/N ∮ B² dl/B_p
           ! ⟨R⁻²⟩ = 1/N ∮ R⁻² dl/B_p
           ! ⟨J_φ⟩ = 1/N ∮ J_φ dl/B_p
+          ! I_p   = 1/μ₀ ∮ B_p⋅dl
           norm         = norm         + dl*(1/B_pol1 + 1/B_pol0)/2
           dAdpsi       = dAdpsi       + dl*(1/(B_pol1*R(j+1)) + 1/(B_pol0*R(j)))/2
-          I_p(i)       = I_p(i)       + dl*(B_pol1 + B_pol0)/2 * mu0inv
           B_avg(i)     = B_avg(i)     + dl*(B_tot1/B_pol1 + B_tot0/B_pol0)/2
           B2_avg       = B2_avg       + dl*(B_tot1**2/B_pol1 + B_tot0**2/B_pol0)/2
           R2_inv_avg   = R2_inv_avg   + dl*(1/(B_pol1*R(j+1)**2) + 1/(B_pol0*R(j)**2))/2
           J_phi_avg(i) = J_phi_avg(i) + dl*(J_phi0/(B_pol0*R(j)) + J_phi1/(B_pol1*R(j+1)))/2
+          I_p(i)       = I_p(i)       + dl*(B_pol1 + B_pol0)/2 * dir * mu0inv
 
           ! Update the previous values
           J_phi0 = J_phi1
@@ -255,7 +259,7 @@ contains
       !
       !   q(ρ) = F/2π ∮ 1/R² dl/B_p = F/2π N ⟨R⁻²⟩
       !
-      q(i) = abs(fpol/(2*pi) * norm * R2_inv_avg)
+      q(i) = fpol/(2*pi) * norm * R2_inv_avg
 
       ! Using Φ=Φ_a⋅ρ_t² and q=1/2π⋅∂Φ/∂ψ (see gray_equil.f90), we have
       !