from enum import Enum, auto 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). gXX dX^2 + 2 gXt dX dt - gtt dt^2 = 0 => dx / dt = (-gXt ± √(gXt^2 - gXX gtt)) / gXX """ 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] # terminate integration on reaching the outer boundaries def _events_bnd(spatial_coords): def event_x_bound_upper(t, x, sign, *args): return x[0] - spatial_coords[-1] event_x_bound_upper.terminal = True event_x_bound_upper.direction = 1.0 def event_x_bound_lower(t, x, sign, *args): return x[0] - spatial_coords[0] event_x_bound_lower.terminal = True event_x_bound_lower.direction = -1.0 return [event_x_bound_upper, event_x_bound_lower] def _kernel_time(times, origin, sign, gXX, gXt, gtt, events): idx = np.where(times == origin)[0][0] t_fwd = times[idx:] t_back = times[:idx + 1][::-1] ray_fwd = [0.0] if len(t_fwd) > 1: sol = solve_ivp(_photon_dXdt, (t_fwd[0], t_fwd[-1]), (0,), method = 'RK45', t_eval = t_fwd, args = (sign, gXX, gtt, gXt), dense_output = True, events = events, rtol = 1e-6, atol = 1e-8) ray_fwd = sol.y[0] ray_fwd = np.concatenate((ray_fwd, [ray_fwd[-1]] * (len(t_fwd) - len(ray_fwd)))) ray_back = [] if len(t_back) > 1: sol = solve_ivp(_photon_dXdt, (t_back[0], t_back[-1]), (0,), method = 'RK45', t_eval = t_back, args = (sign, gXX, gtt, gXt), dense_output = True, events = events, rtol = 1e-6, atol = 1e-8) ray_back = sol.y[0] ray_back = np.concatenate((ray_back, [ray_back[-1]] * (len(t_back) - len(ray_back)))) return np.concatenate((ray_back[::-1][:-1], ray_fwd)) def _kernel_space(times, origin, sign, gXX, gXt, gtt, events): sol = solve_ivp(_photon_dXdt, (times[0], times[-1]), (origin,), method = 'RK45', t_eval = times, args = (sign, gXX, gtt, gXt), dense_output = True, events = events, rtol = 1e-6, atol = 1e-8) t, x = sol.t, sol.y[0] if len(t) < len(times): x_ext = np.empty_like(times) x_ext[:x.shape[0]] = x x_ext[x.shape[0]:] = x[-1] + sign * (times[x.shape[0]:] - t[-1]) x = x_ext return x class Curves(Enum): SPATIAL_FORWARD = auto() "photons are traced forward in time from (t=times[0], X=spatial_coords[i])" SPATIAL_BACK = auto() "photons are traced backward in time from (t=times[-1], X=spatial_coords[i])" TEMPORAL = auto() """ photons are traced both forward and backward in time from (t=times[i], X=0); the forward and backward segment are combined to form the resulting curve """ def null_curves(times, spatial_coords, gXX, gXt, gtt, kind = Curves.SPATIAL_FORWARD, interp_order = 6): """ Compute null curves along a given axis. Shoot a null ray from each point in spatial_coords, in the positive and negative spatial 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 :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]. """ 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) \ if gXt is not None else lambda t, X: 0.0 gtt_interp = interp.Interp2D_C([times[0], spatial_coords[0]], [dt, dX], gtt, interp_order) if kind == Curves.TEMPORAL: kernel = _kernel_time origins = times else: kernel = _kernel_space origins = spatial_coords ray_times = times[::-1] if kind == Curves.SPATIAL_BACK else times events = _events_bnd(spatial_coords) rays_pos = np.empty((origins.shape[0], times.shape[0])) rays_neg = np.empty_like(rays_pos) for tgt, sign in ((rays_pos, 1.0), (rays_neg, -1.0)): for i, origin in enumerate(origins): tgt[i] = kernel(ray_times, origin, sign, gXX_interp, gXt_interp, gtt_interp, events) return (ray_times, rays_pos, rays_neg) def null_coordinates(times, spatial_coords, u_rays, v_rays, gXX, gXt, gtt, kind = Curves.SPATIAL_FORWARD): """ 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 = null_curves(times, spatial_coords, gXX, gXt, gtt, kind = kind) 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] = interp1d(Xu, u_rays, fill_value = 'extrapolate')(spatial_coords) v_of_tx[i] = interp1d(Xv, v_rays, fill_value = 'extrapolate')(spatial_coords) return (u_of_tx, v_of_tx)