doublenull: add a parameter for integrating backwards in time

1 files changed, 17 insertions, 12 deletions

diff --git a/doublenull.py b/doublenull.py
index a3f5caf..5aafc06 100644
--- a/doublenull.py
+++ b/doublenull.py
@@ -2,7 +2,7 @@ import numpy as np
 import scipy.interpolate as interp
 from scipy.integrate import solve_ivp
 
-def calc_null_curves(times, spatial_coords, gXX, gXt, gtt):
+def calc_null_curves(times, spatial_coords, gXX, gXt, gtt, reverse = False):
     """
     Compute null curves along a given axis.
 
@@ -19,16 +19,21 @@ def calc_null_curves(times, spatial_coords, gXX, gXt, gtt):
                            (t=times[i], X=spatial_coords[j]).
     :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
-    :return: tuple containing two 2D arrays with, respectively, the positive and
-             negative-direction rays. rays[i, j] contains the X-coordinate of
-             the ray shot from (t=times[0], X=spatial_coords[i]) at time
-             t=times[j]
+    :param bool reverse:   when true, the null curves are traced from the last
+                           time coordinate backwards in time
+    :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      = interp.RectBivariateSpline(times, spatial_coords, gXX)
     gXt_interp      = interp.RectBivariateSpline(times, spatial_coords, gXt)
     gtt_interp      = interp.RectBivariateSpline(times, spatial_coords, gtt)
 
+    if reverse:
+        ray_times = times[::-1]
+    else:
+        ray_times = times
 
     # terminate integration on reaching the outer boundaries
     def event_x_bound_upper(t, x, sign):
@@ -56,22 +61,22 @@ def calc_null_curves(times, spatial_coords, gXX, gXt, gtt):
     rays_neg = np.empty_like(rays_pos)
     for j, X0 in enumerate(spatial_coords):
         for tgt, sign in ((rays_pos, 1.0), (rays_neg, -1.0)):
-            ret = solve_ivp(dXdt, (times[0], times[-1]), (X0,),
-                            method = 'RK45', t_eval = times, args = (sign,),
+            ret = solve_ivp(dXdt, (ray_times[0], ray_times[-1]), (X0,),
+                            method = 'RK45', t_eval = ray_times, args = (sign,),
                             dense_output = True, events = events, rtol = 1e-6, atol = 1e-8)
 
             t, x = ret.t, ret.y[0]
 
-            if len(t) < len(times):
-                x_ext = np.empty_like(times)
+            if len(t) < len(ray_times):
+                x_ext = np.empty_like(ray_times)
                 x_ext[:x.shape[0]] = x
-                x_ext[x.shape[0]:] = x[-1] + sign * (times[x.shape[0]:] - t[-1])
+                x_ext[x.shape[0]:] = x[-1] + sign * (ray_times[x.shape[0]:] - t[-1])
 
                 x = x_ext
 
             tgt[j] = x
 
-    return (rays_pos, rays_neg)
+    return (ray_times, rays_pos, rays_neg)
 
 def calc_null_coordinates(times, spatial_coords, u_rays, v_rays,
                           gXX, gXt, gtt):
@@ -97,7 +102,7 @@ def calc_null_coordinates(times, spatial_coords, u_rays, v_rays,
              v as functions of t and X. u/v[i, j] is the value of u/v at
              t=times[i], X=spatial_coords[j].
     """
-    X_of_ut, X_of_vt = calc_null_curves(times, spatial_coords, gXX, gXt, gtt)
+    _, X_of_ut, X_of_vt = calc_null_curves(times, spatial_coords, gXX, gXt, gtt)
 
     u_of_tx = np.empty((times.shape[0], spatial_coords.shape[0]))
     v_of_tx = np.empty_like(u_of_tx)