From fd5a437c1f34593031b23e5f42f69d102fc0b889 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Mon, 15 Aug 2022 16:26:57 +0200 Subject: doublenull: add a parameter for integrating backwards in time --- doublenull.py | 29 +++++++++++++++++------------ 1 file 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) -- cgit v1.2.3