Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <cstddef>
       7              : 
       8              : #include "DataStructures/Tensor/TypeAliases.hpp"
       9              : #include "Utilities/Gsl.hpp"
      10              : 
      11              : namespace gr {
      12              : 
      13              : /*!
      14              :  * \brief First-order formulation of the geodesic equation for null geodesics
      15              :  * that is suitable for ray tracing.
      16              :  *
      17              :  * This is the formulation of the geodesic equation that is originally used for
      18              :  * ray tracing in \cite Bohn:2014xxa (Eq. (4)) and \cite Bohn:2016afc (Eq. (6)).
      19              :  *
      20              :  * For null rays with position $x^\mu$, proper time (or affine geodesic
      21              :  * parameter) $\tau$, and four-momentum $p^\mu = dx^\mu / d\tau$, we first
      22              :  * define the momentum variable
      23              :  * \begin{equation}
      24              :  *   \Pi_i = \frac{p_i}{\alpha p^0} = \frac{p_i}{\sqrt{\gamma^{jk} p_j p_k}}
      25              :  *   \text{.}
      26              :  * \end{equation}
      27              :  * Here we work in a 3+1 decomposition of spacetime with time coordinate $t$,
      28              :  * lapse $\alpha$, shift $\beta^i$, spatial metric $\gamma_{ij}$, and
      29              :  * extrinsic curvature $K_{ij}$.
      30              :  * Then, the geodesic equation is given by
      31              :  * \begin{align}
      32              :  *   \frac{d \Pi_i}{d t} &= -\partial_i \alpha + (\Pi^j \partial_j \alpha
      33              :  *     -\alpha K_{jk}\Pi^j\Pi^k) \Pi_i) + \Pi_k \partial_i \beta^k
      34              :  *     - \frac{1}{2} \alpha \Pi_j\Pi_k \partial_i \gamma^{jk} \\
      35              :  *   \frac{d x^i}{d t} &= \alpha \Pi^i - \beta^i
      36              :  *   \text{.}
      37              :  * \end{align}
      38              :  *
      39              :  * This function also computes the evolution of the additional redshift variable
      40              :  * $\ln(\alpha p^0)$ (Eq. (5) in \cite Bohn:2014xxa) as
      41              :  * \begin{equation}
      42              :  *   \frac{d \ln(\alpha p^0)}{d t} = -\Pi^i \partial_i \alpha
      43              :  *     + \alpha K_{ij} \Pi^i \Pi^j
      44              :  *   \text{.}
      45              :  * \end{equation}
      46              :  */
      47              : template <typename DataType, size_t Dim, typename Frame>
      48            1 : void geodesic_equation(
      49              :     // Output time derivs
      50              :     gsl::not_null<tnsr::I<DataType, Dim, Frame>*> dt_x,
      51              :     gsl::not_null<tnsr::i<DataType, Dim, Frame>*> dt_pi,
      52              :     gsl::not_null<Scalar<DataType>*> dt_lnp0,
      53              :     // Current state
      54              :     const tnsr::I<DataType, Dim, Frame>& x,
      55              :     const tnsr::i<DataType, Dim, Frame>& pi, const Scalar<DataType>& lnp0,
      56              :     // Background spacetime
      57              :     const Scalar<DataType>& lapse,
      58              :     const tnsr::i<DataType, Dim, Frame>& deriv_lapse,
      59              :     const tnsr::I<DataType, Dim, Frame>& shift,
      60              :     const tnsr::iJ<DataType, Dim, Frame>& deriv_shift,
      61              :     const tnsr::II<DataType, Dim, Frame>& inv_spatial_metric,
      62              :     const tnsr::iJJ<DataType, Dim, Frame>& deriv_inv_spatial_metric,
      63              :     const tnsr::ii<DataType, Dim, Frame>& extrinsic_curvature);
      64              : 
      65              : }  // namespace gr