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