doublenull: use interp.py for 2D interpolation

3 files changed, 23 insertions, 7 deletions

diff --git a/doublenull.py b/doublenull.py
index 622f156..576dc32 100644
--- a/doublenull.py
+++ b/doublenull.py
@@ -4,6 +4,8 @@ import numpy as np
 from scipy.interpolate import RectBivariateSpline, interp1d
 from scipy.integrate import solve_ivp
 
+from . import interp
+
 def _photon_dXdt(t, coord, sign, gXX, gtt, gXt):
     """
     Null curve equation RHS (not a geodesic - not affinely parametrized).
@@ -89,7 +91,7 @@ class Curves(Enum):
     """
 
 def null_curves(times, spatial_coords, gXX, gXt, gtt,
-                kind = Curves.SPATIAL_FORWARD):
+                kind = Curves.SPATIAL_FORWARD, interp_order = 6):
     """
     Compute null curves along a given axis.
 
@@ -107,14 +109,19 @@ def null_curves(times, spatial_coords, gXX, gXt, gtt,
     :param array_like gXt: same as gXX, but for the Xt component of the metric
     :param array_like gtt: same as gXX, but for the tt component of the metric
     :param Curves kind:    specifies where to integrate the curves from
+    :param int interp_order: Order of interpolation used for metric quantities.
+
     :return: Tuple of (ray_times, rays_pos, rays_neg). rays_*[i, j] contains the
              X-coordinate of the ray shot from (t=ray_times[0],
              X=spatial_coords[i]) at time t=ray_times[j].
     """
 
-    gXX_interp = RectBivariateSpline(times, spatial_coords, gXX)
-    gXt_interp = RectBivariateSpline(times, spatial_coords, gXt)
-    gtt_interp = RectBivariateSpline(times, spatial_coords, gtt)
+    dt = times[1]          - times[0]
+    dX = spatial_coords[1] - spatial_coords[0]
+
+    gXX_interp = interp.Interp2D_C([times[0], spatial_coords[0]], [dt, dX], gXX, interp_order)
+    gXt_interp = interp.Interp2D_C([times[0], spatial_coords[0]], [dt, dX], gXt, interp_order)
+    gtt_interp = interp.Interp2D_C([times[0], spatial_coords[0]], [dt, dX], gtt, interp_order)
 
     if kind == Curves.TEMPORAL:
         kernel  = _kernel_time
diff --git a/test/test_doublenull.npz b/test/test_doublenull.npz
index 4745593..c668973 100644
--- a/test/test_doublenull.npz
+++ b/test/test_doublenull.npz
diff --git a/test/test_doublenull.py b/test/test_doublenull.py
index a61af33..204870d 100644
--- a/test/test_doublenull.py
+++ b/test/test_doublenull.py
@@ -9,10 +9,15 @@ from nr_analysis_axi import doublenull, invars
 
 class TestDoubleNull(TestCase):
 
-    def setUp(self):
-        datafile = os.path.splitext(__file__)[0] + '.npz'
+    def _data_path(self):
+        return os.path.splitext(__file__)[0] + '.npz'
+
+    def update_refs(self, **kwargs):
+        data = dict(self.data) | kwargs
+        np.savez(self._data_path(), **data)
 
-        self.data = np.load(datafile)
+    def setUp(self):
+        self.data = np.load(self._data_path())
 
     def test_null_curves(self):
         rays = doublenull.null_curves(self.data['t'], self.data['x'],
@@ -23,6 +28,8 @@ class TestDoubleNull(TestCase):
         assert_allclose(rays[1], self.data['rays_pos'],  1e-12)
         assert_allclose(rays[2], self.data['rays_neg'],  1e-12)
 
+        #self.update_refs(ray_times = rays[0], rays_pos = rays[1], rays_neg = rays[2])
+
     def test_null_coordinates(self):
         x = self.data['x']
         idx = x.shape[0] // 2
@@ -34,3 +41,5 @@ class TestDoubleNull(TestCase):
                                              -(self.data['alpha'] ** 2))
         assert_allclose(coords[0], self.data['uxt'], 1e-12)
         assert_allclose(coords[1], self.data['vxt'], 1e-12)
+
+        #self.update_refs(uxt = coords[0], vxt = coords[1])