From acacea1a22a709d4551b8c7dfe069c23ef8e9335 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 22 Mar 2018 09:39:56 +0100 Subject: Add code for double-null coordinates. --- doublenull.py | 101 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) create mode 100644 doublenull.py (limited to 'doublenull.py') diff --git a/doublenull.py b/doublenull.py new file mode 100644 index 0000000..01264cf --- /dev/null +++ b/doublenull.py @@ -0,0 +1,101 @@ +import numpy as np +import scipy.interpolate as interp +from scipy.integrate import odeint + +def calc_null_curves(times, spatial_coords, gXX, gXt, gtt): + """ + Compute null curves along a given axis. + + Shoot a null ray from each point in spatial_coords, in the positive and + negative direction and compute its trajectory. + + :param array_like times: 1D array of coordinate times at which the spacetime + curvature is provided + :param array_like spatial_coords: 1D array of spatial coordinates + :param array_like gXX: 2D array containing the values of the XX component + of the spacetime metric, where X is the spatial + coordinate along which the rays are traced. + gXX[i, j] is the value at spacetime point + (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] + """ + + gXX_interp = interp.RectBivariateSpline(times, spatial_coords, gXX) + gXt_interp = interp.RectBivariateSpline(times, spatial_coords, gXt) + gtt_interp = interp.RectBivariateSpline(times, spatial_coords, gtt) + + def _dXdt(coord, t, sign, + coord_min = spatial_coords[0], coord_max = spatial_coords[-1], + gXX = gXX_interp, gtt = gtt_interp, gXt = gXt_interp): + if coord <= coord_min or coord >= coord_max: + return float('nan') + gXt_val = gXt(t, coord) + gXX_val = gXX(t, coord) + gtt_val = gtt(t, coord) + return ((-gXt_val + sign * np.sqrt((gXt_val ** 2) - gtt_val * gXX_val)) / gXX_val).flatten()[0] + + rays_pos = np.empty((spatial_coords.shape[0], times.shape[0])) + 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)): + ray = odeint(_dXdt, X0, times, (sign,))[:, 0] + + # clip the rays to the integration area + k = 0 + while np.isnan(ray[k]): + k += 1 + while k > 0: + ray[k - 1] = ray[k] + k -= 1 + + k = ray.shape[0] - 1 + while np.isnan(ray[k]): + k -= 1 + while k < ray.shape[0] - 1: + ray[k + 1] = ray[k] + k += 1 + + tgt[j] = ray + + return (rays_pos, rays_neg) + +def calc_null_coordinates(times, spatial_coords, u_rays, v_rays, + gXX, gXt, gtt): + """ + Compute double-null coordinates (u, v) as functions of + position and time. + + :param array_like times: 1D array of coordinate times at which the spacetime + curvature is provided + :param array_like spatial_coords: 1D array of spatial coordinates + :param array_like u_rays: 1D array assigning the values of u on the initial + time slice. u_rays[i] is the value of u at + X=spatial_coords[i]. + :param array_like v_rays: same as u_rays, but for v. + :param array_like gXX: 2D array containing the values of the XX component + of the spacetime metric, where X is the spatial + coordinate along which the rays are traced. + gXX[i, j] is the value at spacetime point + (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, values of u and + 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) + + u_of_tx = np.empty((times.shape[0], spatial_coords.shape[0])) + v_of_tx = np.empty_like(u_of_tx) + for i, t in enumerate(times): + Xu = X_of_ut[:, i] + Xv = X_of_vt[:, i] + u_of_tx[i] = interp.interp1d(Xu, u_rays)(spatial_coords) + v_of_tx[i] = interp.interp1d(Xv, v_rays)(spatial_coords) + + return (u_of_tx, v_of_tx) -- cgit v1.2.3