From 358fab45cba4dfdcc0c5d61f93401e5d1f2d0a02 Mon Sep 17 00:00:00 2001
From: Peter Hill <peter.hill@york.ac.uk>
Date: Mon, 18 Aug 2025 10:00:52 +0100
Subject: [PATCH] Replace remove `interp2d` with `RectBivariateSpline`

---
 freegs/equilibrium.py      |  4 +---
 freegs/test_equilibrium.py | 33 +++++++++++++++++++++++++++++++++
 2 files changed, 34 insertions(+), 3 deletions(-)

diff --git a/freegs/equilibrium.py b/freegs/equilibrium.py
index 3159761..6a534b4 100644
--- a/freegs/equilibrium.py
+++ b/freegs/equilibrium.py
@@ -559,7 +559,7 @@ def _updateBoundaryPsi(self, psi=None):
                 Z = np.asarray(self.Z[0, :])
 
                 # psi is transposed due to how FreeGS meshgrids R,Z
-                psi_2d = interpolate.interp2d(x=R, y=Z, z=psi.T)
+                psi_2d = interpolate.RectBivariateSpline(x=R, y=Z, z=psi.T)
 
                 # Get psi at the limit points
                 psi_limit_points = np.zeros(len(Rlimit))
@@ -1414,13 +1414,11 @@ def newDomain(eq, Rmin=None, Rmax=None, Zmin=None, Zmax=None, nx=None, ny=None):
 
     tck = interpolate.bisplrep(ravel(eq.R), ravel(eq.Z), ravel(psi))
     spline = interpolate.RectBivariateSpline(eq.R[:, 0], eq.Z[0, :], psi)
-    f = interpolate.interp2d(eq.R[:, 0], eq.Z[0, :], psi, kind="cubic")
 
     plt.plot(eq.R[:, 10], psi[:, 10], "o")
 
     r = linspace(Rmin, Rmax, 1000)
     z = eq.Z[0, 10]
-    plt.plot(r, f(r, z), label="f")
 
     plt.plot(r, spline(r, z), label="spline")
 
diff --git a/freegs/test_equilibrium.py b/freegs/test_equilibrium.py
index 8b430f7..5a60b36 100644
--- a/freegs/test_equilibrium.py
+++ b/freegs/test_equilibrium.py
@@ -2,6 +2,7 @@
 from . import boundary
 from . import jtor
 from . import picard
+from .machine import TestTokamak
 
 import numpy as np
 
@@ -51,6 +52,38 @@ def test_fixed_boundary_psi():
     assert eq.poloidalBeta() > 0.0
 
 
+def test_fixed_boundary_psi_check_limited():
+    # This is adapted from example 5
+
+    tokamak = TestTokamak()
+
+    eq = equilibrium.Equilibrium(
+        tokamak=tokamak,
+        Rmin=0.1,
+        Rmax=2.0,
+        Zmin=-1.0,
+        Zmax=1.0,
+        nx=65,
+        ny=65,
+        boundary=boundary.fixedBoundary,
+        check_limited=True,
+    )
+
+    profiles = jtor.ConstrainPaxisIp(
+        eq, 1e3, 1e5, 1.0  # Plasma pressure on axis [Pascals]  # Plasma current [Amps]
+    )  # fvac = R*Bt
+
+    # Nonlinear solve
+    picard.solve(eq, profiles)
+
+    psi = eq.psi()
+    assert psi[0, 0] == 0.0  # Boundary is fixed
+    assert psi[32, 32] != 0.0  # Solution is not all zero
+
+    assert eq.psi_bndry == 0.0
+    assert eq.poloidalBeta() > 0.0
+
+
 def test_setSolverVcycle():
     eq = equilibrium.Equilibrium(Rmin=0.1, Rmax=2.0, Zmin=-1.0, Zmax=1.0, nx=65, ny=65)