diff --git a/src/equilibrium.f90 b/src/equilibrium.f90
index 97ffdcb..d5cc50e 100644
--- a/src/equilibrium.f90
+++ b/src/equilibrium.f90
@@ -22,6 +22,12 @@ module equilibrium
     real(wp_) :: B0    ! Magnetic field at the magnetic axis (T)
   end type
 
+  ! A 2D contour in the (R,z) plane
+  type contour
+    real(wp_), allocatable :: R(:)
+    real(wp_), allocatable :: z(:)
+  end type
+
   ! Order of the splines
   integer, parameter :: kspl=3, ksplp=kspl + 1
 
@@ -36,6 +42,11 @@ module equilibrium
   ! Analytic model
   type(analytic_model), save :: model
 
+  ! Contour of the region where the ψ(R,z) spline is monotonically
+  ! increasing. The mapping from ψ to other plasma parameters is
+  ! physically meaningful only inside this domain.
+  type(contour), save :: psi_domain
+
   ! More parameters
   real(wp_), save :: fpolas          ! Poloidal current at the edge (r=a), F_a
   real(wp_), save :: psia            ! Poloidal flux at the edge (r=a), ψ_a
@@ -45,8 +56,8 @@ module equilibrium
   real(wp_), save :: sgnbphi         ! Sign of B_φ (>0 means CCW)
   real(wp_), save :: btrcen          ! used only for jcd_astra def.
   real(wp_), save :: rcen            ! Major radius R₀, computed as fpolas/btrcen
-  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 :: rmnm, rmxm      ! R interval of the equilibrium grid
+  real(wp_), save :: zmnm, zmxm      ! z interval of the equilibrium grid
   real(wp_), save :: zbinf, zbsup    ! z interval of the plasma boundary
 
   private
@@ -538,6 +549,13 @@ contains
     call psi_spline%init_deriv(nr, nz, 2, 0)  ! ∂²ψ/∂R²
     call 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.
+    psi_domain%R = [rmnm, rmnm, rmxm, rmxm]
+    psi_domain%z = [zmnm, zmxm, zmxm, zmnm]
+
     ! Spline interpolation of F(ψ)
 
     ! give a small weight to the last point
@@ -672,9 +690,67 @@ contains
     call setqphi_num(data%psinr, abs(data%qpsi), abs(psia), rhotn)
     call set_rho_spline(sqrt(data%psinr), rhotn)
 
+    ! Compute the domain of the ψ mapping
+    psi_domain%R = data%rbnd
+    psi_domain%z = data%zbnd
+    call rescale_boundary(psi_domain, O=[rmaxis, zmaxis], t0=1.5_wp_)
+
   end subroutine set_equil_spline
 
 
+  subroutine rescale_boundary(cont, 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.
+
+    implicit none
+
+    ! subroutine arguments
+    type(contour), intent(inout) :: cont ! (R,z) contour
+    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
+
+    function s(t)
+      ! Rescriction of ψ(R, z) on the line Q(t) = O + tN
+      implicit none
+
+      real(wp_), intent(in) :: t
+      real(wp_)             :: s, Q(2)
+
+      Q = O + t * N
+      s = psi_spline%eval(Q(1), Q(2))
+    end function
+
+  end subroutine rescale_boundary
+
+
   subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
                         tx,nknt_x,ty,nknt_y,coeff,ierr)
     ! Computes the spline interpolation of a surface when
@@ -834,6 +910,7 @@ contains
     !
     ! Note: all output arguments are optional.
     use gray_params,          only : iequil
+    use reflections,          only : inside
     use const_and_precisions, only : one, pi
 
     implicit none
@@ -960,8 +1037,7 @@ contains
       if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz
     else
       ! Numerical data
-      if (R <= rmxm .and. R >= rmnm .and. &
-          z <= zmxm .and. z >= zmnm) then
+      if (inside(psi_domain%R, psi_domain%z, psi_domain%npoints, R, z)) then
         ! Within the interpolation range
         if (present(psi_n))    psi_n    = psi_spline%eval(R, z)
         if (present(dpsidr))   dpsidr   = psi_spline%deriv(R, z, 1, 0)
@@ -1431,6 +1507,9 @@ contains
     call q_spline%deinit
     call rhop_spline%deinit
     call rhot_spline%deinit
+
+    if (allocated(psi_domain%R)) deallocate(psi_domain%R, psi_domain%z)
+
   end subroutine unset_equil_spline
 
 end module equilibrium