doublenull: use interp.py for 2D interpolation

1 files changed, 11 insertions, 4 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