# Configuration definitions for thorn BaikalX

# Interface definition for thorn BaikalX
implements: BaikalX
inherits:

# Parameter definitions for thorn BaikalX

# Schedule definitions for thorn BaikalX

# Main make.code.defn file for thorn BaikalX
# Source files in this directory
SRCS = 
# Subdirectories containing source files
SUBDIRS = 

% *======================================================================*
%  Cactus Thorn template for ThornGuide documentation
%  Author: Ian Kelley
%  Date: Sun Jun 02, 2002
%  $Header$
%
%  Thorn documentation in the latex file doc/documentation.tex
%  will be included in ThornGuides built with the Cactus make system.
%  The scripts employed by the make system automatically include
%  pages about variables, parameters and scheduling parsed from the
%  relevant thorn CCL files.
%
%  This template contains guidelines which help to assure that your
%  documentation will be correctly added to ThornGuides. More
%  information is available in the Cactus UsersGuide.
%
%  Guidelines:
%   - Do not change anything before the line
%       % START CACTUS THORNGUIDE",
%     except for filling in the title, author, date, etc. fields.
%        - Each of these fields should only be on ONE line.
%        - Author names should be separated with a \\ or a comma.
%   - You can define your own macros, but they must appear after
%     the START CACTUS THORNGUIDE line, and must not redefine standard
%     latex commands.
%   - To avoid name clashes with other thorns, 'labels', 'citations',
%     'references', and 'image' names should conform to the following
%     convention:
%       ARRANGEMENT_THORN_LABEL
%     For example, an image wave.eps in the arrangement CactusWave and
%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
%   - Graphics should only be included using the graphicx package.
%     More specifically, with the "\includegraphics" command.  Do
%     not specify any graphic file extensions in your .tex file. This
%     will allow us to create a PDF version of the ThornGuide
%     via pdflatex.
%   - References should be included with the latex "\bibitem" command.
%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
%   - Do not use \appendix, instead include any appendices you need as
%     standard sections.
%   - For the benefit of our Perl scripts, and for future extensions,
%     please use simple latex.
%
% *======================================================================*
%
% Example of including a graphic image:
%    \begin{figure}[ht]
% 	\begin{center}
%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
% 	\end{center}
% 	\caption{Illustration of this and that}
% 	\label{MyArrangement_MyThorn_MyLabel}
%    \end{figure}
%
% Example of using a label:
%   \label{MyArrangement_MyThorn_MyLabel}
%
% Example of a citation:
%    \cite{MyArrangement_MyThorn_Author99}
%
% Example of including a reference
%   \bibitem{MyArrangement_MyThorn_Author99}
%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
%   1--16. {\tt http://www.nowhere.com/}}
%
% *======================================================================*
% If you are using CVS use this line to give version information
% $Header$
\documentclass{article}
% Use the Cactus ThornGuide style file
% (Automatically used from Cactus distribution, if you have a
%  thorn without the Cactus Flesh download this from the Cactus
%  homepage at www.cactuscode.org)
\usepackage{../../../../doc/latex/cactus}
\begin{document}
% The author of the documentation
\author{Roland Haas \textless rhaas@illinois.edu\textgreater}
% The title of the document (not necessarily the name of the Thorn)
\title{BaikalX}
% the date your document was last changed, if your document is in CVS,
% please use:
%    \date{$ $Date$ $}
\date{December 16 2019}
\maketitle
% Do not delete next line
% START CACTUS THORNGUIDE
% Add all definitions used in this documentation here
%   \def\mydef etc
% Add an abstract for this thorn's documentation
\begin{abstract}
\end{abstract}
% The following sections are suggestive only.
% Remove them or add your own.
\section{Introduction}
\section{Physical System}
\section{Numerical Implementation}
\section{Using This Thorn}
\subsection{Obtaining This Thorn}
\subsection{Basic Usage}
\subsection{Special Behaviour}
\subsection{Interaction With Other Thorns}
\subsection{Examples}
\subsection{Support and Feedback}
\section{History}
\subsection{Thorn Source Code}
\subsection{Thorn Documentation}
\subsection{Acknowledgements}
\begin{thebibliography}{9}
\end{thebibliography}
% Do not delete next line
% END CACTUS THORNGUIDE
\end{document}

# As documented in the NRPy+ tutorial module
#   Tutorial-ADMBSSN_tofrom_4metric.ipynb,
#   this module will construct expressions for
#   ADM or BSSN quantities in terms of the
#   4-metric g4DD, and g4DD/g4UU in terms of
#   ADM/BSSN quantities.
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import sys                        # Standard Python modules for multiplatform OS-level functions
def setup_ADM_quantities(inputvars):
    if inputvars == "ADM":
        gammaDD = ixp.declarerank2("gammaDD", "sym01")
        betaU = ixp.declarerank1("betaU")
        alpha = sp.symbols("alpha", real=True)
    elif inputvars == "BSSN":
        import BSSN.ADM_in_terms_of_BSSN as AitoB
        # Construct gamma_{ij} in terms of cf & gammabar_{ij}
        AitoB.ADM_in_terms_of_BSSN()
        gammaDD = AitoB.gammaDD
        # Next construct beta^i in terms of vet^i and reference metric quantities
        import BSSN.BSSN_quantities as Bq
        Bq.BSSN_basic_tensors()
        betaU = Bq.betaU
        alpha = sp.symbols("alpha", real=True)
    else:
        print("inputvars = " + str(inputvars) + " not supported. Please choose ADM or BSSN.")
        sys.exit(1)
    return gammaDD,betaU,alpha
# g_{mu nu} in terms of BSSN (if inputvars=="BSSN") or ADM (if inputvars=="ADM") variables.
def g4DD_ito_BSSN_or_ADM(inputvars,gammaDD=None,betaU=None,alpha=None):
    # Step 0: Declare g4DD as globals, to make interfacing with other modules/functions easier
    global g4DD
    # Step 1: Set gammaDD, betaU, and alpha if not already input.
    if gammaDD==None and betaU==None and alpha==None:
        gammaDD,betaU,alpha = setup_ADM_quantities(inputvars)
    # Step 2: Compute g4DD = g_{mu nu}:
    # To get \gamma_{\mu \nu} = gamma4DD[mu][nu], we'll need to construct the 4-metric, using Eq. 2.122 in B&S:
    g4DD = ixp.zerorank2(DIM=4)
    # Step 2.a: Compute beta_i via Eq. 2.121 in B&S
    betaD = ixp.zerorank1()
    for i in range(3):
        for j in range(3):
            betaD[i] += gammaDD[i][j] * betaU[j]
    # Step 2.b: Compute beta_i beta^i, the beta contraction.
    beta2 = sp.sympify(0)
    for i in range(3):
        beta2 += betaU[i] * betaD[i]
    # Step 2.c: Construct g4DD via Eq. 2.122 in B&S
    g4DD[0][0] = -alpha ** 2 + beta2
    for mu in range(1, 4):
        g4DD[mu][0] = g4DD[0][mu] = betaD[mu - 1]
    for mu in range(1, 4):
        for nu in range(1, 4):
            g4DD[mu][nu] = gammaDD[mu - 1][nu - 1]
# g^{mu nu} in terms of BSSN (if inputvars=="BSSN") or ADM (if inputvars=="ADM") variables.
def g4UU_ito_BSSN_or_ADM(inputvars,gammaDD=None,betaU=None,alpha=None, gammaUU=None):
    # Step 0: Declare g4UU as globals, to make interfacing with other modules/functions easier
    global g4UU
    # Step 1: Set gammaDD, betaU, and alpha if not already input.
    if gammaDD==None and betaU==None and alpha==None:
        gammaDD,betaU,alpha = setup_ADM_quantities(inputvars)
    # Step 2: Compute g4UU = g_{mu nu}:
    # To get \gamma^{\mu \nu} = gamma4UU[mu][nu], we'll need to use Eq. 2.119 in B&S.
    g4UU = ixp.zerorank2(DIM=4)
    # Step 3: Construct g4UU = g^{mu nu}
    # Step 3.a: Compute gammaUU based on provided gammaDD:
    if gammaUU==None:
        gammaUU, gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)
    # Then evaluate g4UU:
    g4UU = ixp.zerorank2(DIM=4)
    g4UU[0][0] = -1 / alpha ** 2
    for mu in range(1, 4):
        g4UU[0][mu] = g4UU[mu][0] = betaU[mu - 1] / alpha ** 2
    for mu in range(1, 4):
        for nu in range(1, 4):
            g4UU[mu][nu] = gammaUU[mu - 1][nu - 1] - betaU[mu - 1] * betaU[nu - 1] / alpha ** 2
# BSSN (if inputvars=="BSSN") or ADM (if inputvars=="ADM") metric variables in terms of g_{mu nu}
def BSSN_or_ADM_ito_g4DD(inputvars,g4DD=None):
    # Step 0: Declare output variables as globals, to make interfacing with other modules/functions easier
    if inputvars == "ADM":
        global gammaDD, betaU, alpha
    elif inputvars == "BSSN":
        global hDD, cf, vetU, alpha
    else:
        print("inputvars = " + str(inputvars) + " not supported. Please choose ADM or BSSN.")
        sys.exit(1)
    # Step 1: declare g4DD as symmetric rank-4 tensor:
    g4DD_is_input_into_this_function = True
    if g4DD == None:
        g4DD = ixp.declarerank2("g4DD", "sym01", DIM=4)
        g4DD_is_input_into_this_function = False
    # Step 2: Compute gammaDD & betaD
    betaD = ixp.zerorank1()
    gammaDD = ixp.zerorank2()
    for i in range(3):
        betaD[i] = g4DD[0][i]
        for j in range(3):
            gammaDD[i][j] = g4DD[i + 1][j + 1]
    # Step 3: Compute betaU
    # Step 3.a: Compute gammaUU based on provided gammaDD
    gammaUU, gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)
    # Step 3.b: Use gammaUU to raise betaU
    betaU = ixp.zerorank1()
    for i in range(3):
        for j in range(3):
            betaU[i] += gammaUU[i][j] * betaD[j]
    # Step 4:   Compute alpha = sqrt(beta^2 - g_{00}):
    # Step 4.a: Compute beta^2 = beta^k beta_k:
    beta_squared = sp.sympify(0)
    for k in range(3):
        beta_squared += betaU[k] * betaD[k]
    # Step 4.b: alpha = sqrt(beta^2 - g_{00}):
    if g4DD_is_input_into_this_function == False:
        alpha = sp.sqrt(sp.simplify(beta_squared) - g4DD[0][0])
    else:
        alpha = sp.sqrt(beta_squared - g4DD[0][0])
    # Step 5: If inputvars == "ADM", we are finished. Return.
    if inputvars == "ADM":
        return
    # Step 6: If inputvars == "BSSN", convert ADM to BSSN
    import BSSN.BSSN_in_terms_of_ADM as BitoA
    dummyBU = ixp.zerorank1()
    BitoA.gammabarDD_hDD(gammaDD)
    BitoA.cf_from_gammaDD(gammaDD)
    BitoA.betU_vetU(betaU, dummyBU)
    hDD  = BitoA.hDD
    cf   = BitoA.cf
    vetU = BitoA.vetU

# This module performs the conversion between ADM
# spacetime variables in Spherical or Cartesian coordinates
# given as *exact* SymPy expressions (i.e., direct functions
# of r,th,ph or x,y,z), to rescaled BSSN-in-curvilinear
# coordinate quantities, as defined in BSSN_RHSs.py
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step P0: Import needed Python/NRPy+ modules
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
import sys                        # Standard Python module for multiplatform OS-level functions
def Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear(CoordType_in, Sph_r_th_ph_or_Cart_xyz,
        gammaDD_inSphorCart, KDD_inSphorCart, alpha_inSphorCart, betaU_inSphorCart, BU_inSphorCart):
    # This routine converts the ADM variables
    #    $$\left\{\gamma_{ij}, K_{ij}, \alpha, \beta^i\right\}$$
    #    in Spherical or Cartesian basis+coordinates, first to the BSSN variables
    #    in the chosen reference_metric::CoordSystem coordinate system+basis:
    # $$\left\{\bar{\gamma}_{i j},\bar{A}_{i j},\phi, K, \bar{\Lambda}^{i}, \alpha, \beta^i, B^i\right\},$$
    #    and then to the rescaled variables:
    # $$\left\{h_{i j},a_{i j},\phi, K, \lambda^{i}, \alpha, \mathcal{V}^i, \mathcal{B}^i\right\}.$$
    # The ADM & BSSN formalisms only work in 3D; they are 3+1 decompositions of Einstein's equations.
    #    To implement axisymmetry or spherical symmetry, simply set all spatial derivatives in
    #    the relevant angular directions to zero; DO NOT SET DIM TO ANYTHING BUT 3.
    # Step P1: Set spatial dimension (must be 3 for BSSN)
    DIM = 3
    
    # Step P2: Copy gammaSphDD_in to gammaSphDD, KSphDD_in to KSphDD, etc.
    #    This ensures that the input arrays are not modified below;
    #    modifying them would result in unexpected effects outside
    #    this function.
    alphaSphorCart = alpha_inSphorCart
    betaSphorCartU   = ixp.zerorank1()
    BSphorCartU      = ixp.zerorank1()
    gammaSphorCartDD = ixp.zerorank2()
    KSphorCartDD     = ixp.zerorank2()
    for i in range(DIM):
        betaSphorCartU[i] = betaU_inSphorCart[i]
        BSphorCartU[i]    = BU_inSphorCart[i]
        for j in range(DIM):
            gammaSphorCartDD[i][j] = gammaDD_inSphorCart[i][j]
            KSphorCartDD[i][j]     = KDD_inSphorCart[i][j]
    # Make sure that rfm.reference_metric() has been called.
    #    We'll need the variables it defines throughout this module.
    if rfm.have_already_called_reference_metric_function == False:
        print("Error. Called Convert_Spherical_ADM_to_BSSN_curvilinear() without")
        print("       first setting up reference metric, by calling rfm.reference_metric().")
        sys.exit(1)
    
    # Step 1: All input quantitiefs are in terms of r,th,ph or x,y,z. We want them in terms
    #         of xx0,xx1,xx2, so here we call sympify_integers__replace_rthph() to replace
    #         r,th,ph or x,y,z, respectively, with the appropriate functions of xx0,xx1,xx2
    #         as defined for this particular reference metric in reference_metric.py's 
    #         xxSph[] or xxCart[], respectively:
    # Note that substitution only works when the variable is not an integer. Hence the 
    #         if isinstance(...,...) stuff:
    def sympify_integers__replace_rthph_or_Cartxyz(obj, rthph_or_xyz, rthph_or_xyz_of_xx):
        if isinstance(obj, int):
            return sp.sympify(obj)
        else:
            return obj.subs(rthph_or_xyz[0], rthph_or_xyz_of_xx[0]).\
                subs(rthph_or_xyz[1], rthph_or_xyz_of_xx[1]).\
                subs(rthph_or_xyz[2], rthph_or_xyz_of_xx[2])
    r_th_ph_or_Cart_xyz_of_xx = []
    if CoordType_in == "Spherical":
        r_th_ph_or_Cart_xyz_of_xx = rfm.xxSph
    elif CoordType_in == "Cartesian":
        r_th_ph_or_Cart_xyz_of_xx = rfm.xxCart
    else:
        print("Error: Can only convert ADM Cartesian or Spherical initial data to BSSN Curvilinear coords.")
        sys.exit(1)
    alphaSphorCart = sympify_integers__replace_rthph_or_Cartxyz(
        alphaSphorCart, Sph_r_th_ph_or_Cart_xyz, r_th_ph_or_Cart_xyz_of_xx)
    for i in range(DIM):
        betaSphorCartU[i] = sympify_integers__replace_rthph_or_Cartxyz(
            betaSphorCartU[i], Sph_r_th_ph_or_Cart_xyz, r_th_ph_or_Cart_xyz_of_xx)
        BSphorCartU[i]    = sympify_integers__replace_rthph_or_Cartxyz(
            BSphorCartU[i],    Sph_r_th_ph_or_Cart_xyz, r_th_ph_or_Cart_xyz_of_xx)
        for j in range(DIM):
            gammaSphorCartDD[i][j] = sympify_integers__replace_rthph_or_Cartxyz(
                gammaSphorCartDD[i][j], Sph_r_th_ph_or_Cart_xyz, r_th_ph_or_Cart_xyz_of_xx)
            KSphorCartDD[i][j]     = sympify_integers__replace_rthph_or_Cartxyz(
                KSphorCartDD[i][j],     Sph_r_th_ph_or_Cart_xyz, r_th_ph_or_Cart_xyz_of_xx)
    # Step 2: All ADM initial data quantities are now functions of xx0,xx1,xx2, but
    #         they are still in the Spherical or Cartesian basis. We can now directly apply
    #         Jacobian transformations to get them in the correct xx0,xx1,xx2 basis:
    # alpha is a scalar, so no Jacobian transformation is necessary.
    alpha = alphaSphorCart
    Jac_dUSphorCart_dDrfmUD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            Jac_dUSphorCart_dDrfmUD[i][j] = sp.diff(r_th_ph_or_Cart_xyz_of_xx[i],rfm.xx[j])
    Jac_dUrfm_dDSphorCartUD, dummyDET = ixp.generic_matrix_inverter3x3(Jac_dUSphorCart_dDrfmUD)
    betaU   = ixp.zerorank1()
    BU      = ixp.zerorank1()
    gammaDD = ixp.zerorank2()
    KDD     = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            betaU[i] += Jac_dUrfm_dDSphorCartUD[i][j] * betaSphorCartU[j]
            BU[i]    += Jac_dUrfm_dDSphorCartUD[i][j] * BSphorCartU[j]
            for k in range(DIM):
                for l in range(DIM):
                    gammaDD[i][j] += Jac_dUSphorCart_dDrfmUD[k][i]*Jac_dUSphorCart_dDrfmUD[l][j] * gammaSphorCartDD[k][l]
                    KDD[i][j]     += Jac_dUSphorCart_dDrfmUD[k][i]*Jac_dUSphorCart_dDrfmUD[l][j] *     KSphorCartDD[k][l]
    # Step 3: All ADM quantities were input into this function in the Spherical or Cartesian
    #         basis, as functions of r,th,ph or x,y,z, respectively. In Steps 1 and 2 above,
    #         we converted them to the xx0,xx1,xx2 basis, and as functions of xx0,xx1,xx2.
    #         Here we convert ADM quantities in the "rfm" basis to their BSSN Curvilinear
    #         counterparts:
    import BSSN.BSSN_in_terms_of_ADM as BitoA
    BitoA.gammabarDD_hDD(gammaDD)
    BitoA.trK_AbarDD_aDD(gammaDD, KDD)
    BitoA.LambdabarU_lambdaU__exact_gammaDD(gammaDD)
    BitoA.cf_from_gammaDD(gammaDD)
    BitoA.betU_vetU(betaU, BU)
    # Step 4: Return the BSSN Curvilinear variables in the desired xx0,xx1,xx2
    #         basis, and as functions of the consistent xx0,xx1,xx2 coordinates.
    return BitoA.cf, BitoA.hDD, BitoA.lambdaU, BitoA.aDD, BitoA.trK, alpha, BitoA.vetU, BitoA.betU

# This module performs the conversion between ADM
# spacetime variables in Spherical or Cartesian coordinates
# given as *numerical* expressions (i.e., ADM quantities
# are only given to floating-point precision; e.g., in the
# case of an initial data solver), to rescaled BSSN-in-curvilinear
# coordinate quantities, as defined in BSSN_RHSs.py
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step P1: Initialize core Python/NRPy+ modules
from outputC import *             # NRPy+: Core C code output module
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import finite_difference as fin   # NRPy+: Finite difference C code generation module
import grid as gri                # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
import BSSN.BSSN_quantities as Bq # NRPy+: Computes useful BSSN quantities; e.g., gammabarUU & GammabarUDD needed below
import os, sys                    # Standard Python modules for multiplatform OS-level functions
def Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear(CoordType_in, ADM_input_function_name, 
                                                          Ccodesdir = "BSSN", pointer_to_ID_inputs=False,loopopts=",oldloops"):
    # The ADM & BSSN formalisms only work in 3D; they are 3+1 decompositions of Einstein's equations.
    #    To implement axisymmetry or spherical symmetry, simply set all spatial derivatives in
    #    the relevant angular directions to zero; DO NOT SET DIM TO ANYTHING BUT 3.
    # Step 0: Set spatial dimension (must be 3 for BSSN)
    DIM = 3
    # Step 1: All ADM initial data quantities are now functions of xx0,xx1,xx2, but
    #         they are still in the Spherical or Cartesian basis. We can now directly apply
    #         Jacobian transformations to get them in the correct xx0,xx1,xx2 basis:
    #         All input quantities are in terms of r,th,ph or x,y,z. We want them in terms
    #         of xx0,xx1,xx2, so here we call sympify_integers__replace_rthph() to replace
    #         r,th,ph or x,y,z, respectively, with the appropriate functions of xx0,xx1,xx2
    #         as defined for this particular reference metric in reference_metric.py's
    #         xxSph[] or xxCart[], respectively:
    # Define the input variables:
    gammaSphorCartDD = ixp.declarerank2("gammaSphorCartDD", "sym01")
    KSphorCartDD = ixp.declarerank2("KSphorCartDD", "sym01")
    alphaSphorCart = sp.symbols("alphaSphorCart")
    betaSphorCartU = ixp.declarerank1("betaSphorCartU")
    BSphorCartU = ixp.declarerank1("BSphorCartU")
    # Make sure that rfm.reference_metric() has been called.
    #    We'll need the variables it defines throughout this module.
    if rfm.have_already_called_reference_metric_function == False:
        print("Error. Called Convert_Spherical_ADM_to_BSSN_curvilinear() without")
        print("       first setting up reference metric, by calling rfm.reference_metric().")
        sys.exit(1)
    r_th_ph_or_Cart_xyz_oID_xx = []
    if CoordType_in == "Spherical":
        r_th_ph_or_Cart_xyz_oID_xx = rfm.xxSph
    elif CoordType_in == "Cartesian":
        r_th_ph_or_Cart_xyz_oID_xx = rfm.xxCart
    else:
        print("Error: Can only convert ADM Cartesian or Spherical initial data to BSSN Curvilinear coords.")
        sys.exit(1)
    # Step 2: All ADM initial data quantities are now functions of xx0,xx1,xx2, but
    #         they are still in the Spherical or Cartesian basis. We can now directly apply
    #         Jacobian transformations to get them in the correct xx0,xx1,xx2 basis:
    # alpha is a scalar, so no Jacobian transformation is necessary.
    alpha = alphaSphorCart
    Jac_dUSphorCart_dDrfmUD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            Jac_dUSphorCart_dDrfmUD[i][j] = sp.diff(r_th_ph_or_Cart_xyz_oID_xx[i], rfm.xx[j])
    Jac_dUrfm_dDSphorCartUD, dummyDET = ixp.generic_matrix_inverter3x3(Jac_dUSphorCart_dDrfmUD)
    
    betaU = ixp.zerorank1()
    BU = ixp.zerorank1()
    gammaDD = ixp.zerorank2()
    KDD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            betaU[i] += Jac_dUrfm_dDSphorCartUD[i][j] * betaSphorCartU[j]
            BU[i] += Jac_dUrfm_dDSphorCartUD[i][j] * BSphorCartU[j]
            for k in range(DIM):
                for l in range(DIM):
                    gammaDD[i][j] += Jac_dUSphorCart_dDrfmUD[k][i] * Jac_dUSphorCart_dDrfmUD[l][j] * \
                                     gammaSphorCartDD[k][l]
                    KDD[i][j] += Jac_dUSphorCart_dDrfmUD[k][i] * Jac_dUSphorCart_dDrfmUD[l][j] * KSphorCartDD[k][l]
    # Step 3: All ADM quantities were input into this function in the Spherical or Cartesian
    #         basis, as functions of r,th,ph or x,y,z, respectively. In Steps 1 and 2 above,
    #         we converted them to the xx0,xx1,xx2 basis, and as functions of xx0,xx1,xx2.
    #         Here we convert ADM quantities in the "rfm" basis to their BSSN Curvilinear
    #         counterparts, for all BSSN quantities *except* lambda^i:
    import BSSN.BSSN_in_terms_of_ADM as BitoA
    BitoA.gammabarDD_hDD(gammaDD)
    BitoA.trK_AbarDD_aDD(gammaDD, KDD)
    BitoA.cf_from_gammaDD(gammaDD)
    BitoA.betU_vetU(betaU, BU)
    hDD = BitoA.hDD
    trK = BitoA.trK
    aDD = BitoA.aDD
    cf = BitoA.cf
    vetU = BitoA.vetU
    betU = BitoA.betU
    # Step 4: Compute $\bar{\Lambda}^i$ (Eqs. 4 and 5 of
    #         [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf)),
    #         from finite-difference derivatives of rescaled metric
    #         quantities $h_{ij}$:
    # \bar{\Lambda}^i = \bar{\gamma}^{jk}\left(\bar{\Gamma}^i_{jk} - \hat{\Gamma}^i_{jk}\right).
    # The reference_metric.py module provides us with analytic expressions for
    #         $\hat{\Gamma}^i_{jk}$, so here we need only compute
    #         finite-difference expressions for $\bar{\Gamma}^i_{jk}$, based on
    #         the values for $h_{ij}$ provided in the initial data. Once
    #         $\bar{\Lambda}^i$ has been computed, we apply the usual rescaling
    #         procedure:
    # \lambda^i = \bar{\Lambda}^i/\text{ReU[i]},
    # and then output the result to a C file using the NRPy+
    #         finite-difference C output routine.
    # We will need all BSSN gridfunctions to be defined, as well as
    #     expressions for gammabarDD_dD in terms of exact derivatives of
    #     the rescaling matrix and finite-difference derivatives of
    #     hDD's. This functionality is provided by BSSN.BSSN_unrescaled_and_barred_vars,
    #     which we call here to overwrite above definitions of gammabarDD,gammabarUU, etc.
    Bq.gammabar__inverse_and_derivs() # Provides gammabarUU and GammabarUDD
    gammabarUU    = Bq.gammabarUU
    GammabarUDD   = Bq.GammabarUDD
    # Next evaluate \bar{\Lambda}^i, based on GammabarUDD above and GammahatUDD
    #       (from the reference metric):
    LambdabarU = ixp.zerorank1()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                LambdabarU[i] += gammabarUU[j][k] * (GammabarUDD[i][j][k] - rfm.GammahatUDD[i][j][k])
    # Finally apply rescaling:
    # lambda^i = Lambdabar^i/\text{ReU[i]}
    lambdaU = ixp.zerorank1()
    for i in range(DIM):
        lambdaU[i] = LambdabarU[i] / rfm.ReU[i]
    if ADM_input_function_name == "DoNotOutputADMInputFunction":
        return hDD,aDD,trK,vetU,betU,alpha,cf,lambdaU
    # Step 5.A: Output files containing finite-differenced lambdas.
    outCparams = "preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False"
    lambdaU_expressions = [lhrh(lhs=gri.gfaccess("in_gfs", "lambdaU0"), rhs=lambdaU[0]),
                           lhrh(lhs=gri.gfaccess("in_gfs", "lambdaU1"), rhs=lambdaU[1]),
                           lhrh(lhs=gri.gfaccess("in_gfs", "lambdaU2"), rhs=lambdaU[2])]
    desc = "Output lambdaU[i] for BSSN, built using finite-difference derivatives."
    name = "ID_BSSN_lambdas"
    params = "const paramstruct *restrict params,REAL *restrict xx[3],REAL *restrict in_gfs"
    preloop = ""
    opts = ""
    idx4replace = "IDX4S"
    if "oldloops" in loopopts:
        params = "const int Nxx[3],const int Nxx_plus_2NGHOSTS[3],REAL *xx[3],const REAL dxx[3],REAL *in_gfs"
        opts = "DisableCparameters"
        idx4replace = "IDX4"
        preloop = """
const REAL invdx0 = 1.0/dxx[0];
const REAL invdx1 = 1.0/dxx[1];
const REAL invdx2 = 1.0/dxx[2];
"""
    outCfunction(
        outfile=os.path.join(Ccodesdir, name + ".h"), desc=desc, name=name, params=params,
        preloop=preloop,
        body=fin.FD_outputC("returnstring", lambdaU_expressions, outCparams).replace("IDX4",idx4replace),
        loopopts="InteriorPoints,Read_xxs"+loopopts, opts=opts)
    # Step 5: Output all ADM-to-BSSN expressions to a C function. This function
    #         must first call the ID_ADM_SphorCart() defined above. Using these
    #         Spherical or Cartesian data, it sets up all quantities needed for
    #         BSSNCurvilinear initial data, *except* $\lambda^i$, which must be
    #         computed from numerical data using finite-difference derivatives.
    ID_inputs_param = "ID_inputs other_inputs,"
    if pointer_to_ID_inputs == True:
        ID_inputs_param = "ID_inputs *other_inputs,"
    desc = "Write BSSN variables in terms of ADM variables at a given point xx0,xx1,xx2"
    name = "ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs"
    opts = ""
    params = "const paramstruct *restrict params, "
    if "oldloops" in loopopts:
        opts = "DisableCparameters"
        params = ""
    params += "const REAL xx0xx1xx2[3]," + ID_inputs_param + """
                    REAL *hDD00,REAL *hDD01,REAL *hDD02,REAL *hDD11,REAL *hDD12,REAL *hDD22,
                    REAL *aDD00,REAL *aDD01,REAL *aDD02,REAL *aDD11,REAL *aDD12,REAL *aDD22,
                    REAL *trK, 
                    REAL *vetU0,REAL *vetU1,REAL *vetU2,
                    REAL *betU0,REAL *betU1,REAL *betU2,
                    REAL *alpha,  REAL *cf"""
    outCparams = "preindent=1,outCverbose=False,includebraces=False"
    outCfunction(
        outfile=os.path.join(Ccodesdir, name + ".h"), desc=desc, name=name, params=params,
        body="""
      REAL gammaSphorCartDD00,gammaSphorCartDD01,gammaSphorCartDD02,
           gammaSphorCartDD11,gammaSphorCartDD12,gammaSphorCartDD22;
      REAL KSphorCartDD00,KSphorCartDD01,KSphorCartDD02,
           KSphorCartDD11,KSphorCartDD12,KSphorCartDD22;
      REAL alphaSphorCart,betaSphorCartU0,betaSphorCartU1,betaSphorCartU2;
      REAL BSphorCartU0,BSphorCartU1,BSphorCartU2;
      const REAL xx0 = xx0xx1xx2[0];
      const REAL xx1 = xx0xx1xx2[1];
      const REAL xx2 = xx0xx1xx2[2];
      REAL xyz_or_rthph[3];\n""" +
             outputC(r_th_ph_or_Cart_xyz_oID_xx[0:3], ["xyz_or_rthph[0]", "xyz_or_rthph[1]", "xyz_or_rthph[2]"],
                     "returnstring",
                     outCparams + ",CSE_enable=False") + "      " + ADM_input_function_name + """(xyz_or_rthph, other_inputs,
                       &gammaSphorCartDD00,&gammaSphorCartDD01,&gammaSphorCartDD02,
                       &gammaSphorCartDD11,&gammaSphorCartDD12,&gammaSphorCartDD22,
                       &KSphorCartDD00,&KSphorCartDD01,&KSphorCartDD02,
                       &KSphorCartDD11,&KSphorCartDD12,&KSphorCartDD22,
                       &alphaSphorCart,&betaSphorCartU0,&betaSphorCartU1,&betaSphorCartU2,
                       &BSphorCartU0,&BSphorCartU1,&BSphorCartU2);
      // Next compute all rescaled BSSN curvilinear quantities:\n""" +
             outputC([hDD[0][0], hDD[0][1], hDD[0][2], hDD[1][1], hDD[1][2], hDD[2][2],
                      aDD[0][0], aDD[0][1], aDD[0][2], aDD[1][1], aDD[1][2], aDD[2][2],
                      trK, vetU[0], vetU[1], vetU[2], betU[0], betU[1], betU[2],
                      alpha, cf],
                     ["*hDD00", "*hDD01", "*hDD02", "*hDD11", "*hDD12", "*hDD22",
                      "*aDD00", "*aDD01", "*aDD02", "*aDD11", "*aDD12", "*aDD22",
                      "*trK", "*vetU0", "*vetU1", "*vetU2", "*betU0", "*betU1", "*betU2",
                      "*alpha", "*cf"], "returnstring", params=outCparams),
        opts = opts)
    # Step 5.a: Output the driver function for the above
    #           function ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs()
    # Next write the driver function for ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs():
    desc = """Driver function for ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs(), 
which writes BSSN variables in terms of ADM variables at a given point xx0,xx1,xx2"""
    name = "ID_BSSN__ALL_BUT_LAMBDAs"
    params = "const paramstruct *restrict params,REAL *restrict xx[3]," + ID_inputs_param + "REAL *in_gfs"
    opts = ""
    funccallparams = "params, "
    idx3replace   = "IDX3S"
    idx4ptreplace = "IDX4ptS"
    if "oldloops" in loopopts:
        params = "const int Nxx_plus_2NGHOSTS[3],REAL *xx[3]," + ID_inputs_param + "REAL *in_gfs"
        opts = "DisableCparameters"
        funccallparams = ""
        idx3replace   = "IDX3"
        idx4ptreplace = "IDX4pt"
    outCfunction(
        outfile=os.path.join(Ccodesdir, name + ".h"), desc=desc, name=name, params=params,
        body="""
const int idx = IDX3(i0,i1,i2);
const REAL xx0xx1xx2[3] = {xx0,xx1,xx2};
ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs(""".replace("IDX3",idx3replace)+funccallparams+"""xx0xx1xx2,other_inputs,
                    &in_gfs[IDX4pt(HDD00GF,idx)],&in_gfs[IDX4pt(HDD01GF,idx)],&in_gfs[IDX4pt(HDD02GF,idx)],
                    &in_gfs[IDX4pt(HDD11GF,idx)],&in_gfs[IDX4pt(HDD12GF,idx)],&in_gfs[IDX4pt(HDD22GF,idx)],
                    &in_gfs[IDX4pt(ADD00GF,idx)],&in_gfs[IDX4pt(ADD01GF,idx)],&in_gfs[IDX4pt(ADD02GF,idx)],
                    &in_gfs[IDX4pt(ADD11GF,idx)],&in_gfs[IDX4pt(ADD12GF,idx)],&in_gfs[IDX4pt(ADD22GF,idx)],
                    &in_gfs[IDX4pt(TRKGF,idx)],
                    &in_gfs[IDX4pt(VETU0GF,idx)],&in_gfs[IDX4pt(VETU1GF,idx)],&in_gfs[IDX4pt(VETU2GF,idx)],
                    &in_gfs[IDX4pt(BETU0GF,idx)],&in_gfs[IDX4pt(BETU1GF,idx)],&in_gfs[IDX4pt(BETU2GF,idx)],
                    &in_gfs[IDX4pt(ALPHAGF,idx)],&in_gfs[IDX4pt(CFGF,idx)]);
""".replace("IDX4pt",idx4ptreplace),
        loopopts="AllPoints,Read_xxs"+loopopts, opts=opts)

# As documented in the NRPy+ tutorial module
#   Tutorial-ADM_in_terms_of_BSSN.ipynb,
#   this module will construct expressions for ADM
#   quantities in terms of BSSN quantities.
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step 1.a: import all needed modules from NRPy+:
from outputC import *             # NRPy+: Core C code output module
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
import sys                        # Standard Python module for multiplatform OS-level functions
def ADM_in_terms_of_BSSN():
    global gammaDD, gammaDDdD, gammaDDdDD, gammaUU, detgamma, GammaUDD, KDD, KDDdD
    # Step 1.c: Given the chosen coordinate system, set up
    #           corresponding reference metric and needed
    #           reference metric quantities
    # The following function call sets up the reference metric
    #    and related quantities, including rescaling matrices ReDD,
    #    ReU, and hatted quantities.
    rfm.reference_metric()
    # Step 1.d: Set spatial dimension (must be 3 for BSSN, as BSSN is
    #           a 3+1-dimensional decomposition of the general
    #           relativistic field equations)
    DIM = 3
    # Step 1.e: Import all basic (unrescaled) BSSN scalars & tensors
    import BSSN.BSSN_quantities as Bq
    Bq.BSSN_basic_tensors()
    gammabarDD = Bq.gammabarDD
    cf         = Bq.cf
    AbarDD     = Bq.AbarDD
    trK        = Bq.trK
    Bq.gammabar__inverse_and_derivs()
    gammabarDD_dD  = Bq.gammabarDD_dD
    gammabarDD_dDD = Bq.gammabarDD_dDD
    Bq.AbarUU_AbarUD_trAbar_AbarDD_dD()
    AbarDD_dD = Bq.AbarDD_dD
    # Step 2: The ADM three-metric gammaDD and its
    #         derivatives in terms of BSSN quantities.
    gammaDD = ixp.zerorank2()
    exp4phi = sp.sympify(0)
    if par.parval_from_str("EvolvedConformalFactor_cf") == "phi":
        exp4phi = sp.exp(4 * cf)
    elif par.parval_from_str("EvolvedConformalFactor_cf") == "chi":
        exp4phi = (1 / cf)
    elif par.parval_from_str("EvolvedConformalFactor_cf") == "W":
        exp4phi = (1 / cf ** 2)
    else:
        print("Error EvolvedConformalFactor_cf type = \"" + par.parval_from_str("EvolvedConformalFactor_cf") + "\" unknown.")
        sys.exit(1)
    for i in range(DIM):
        for j in range(DIM):
            gammaDD[i][j] = exp4phi * gammabarDD[i][j]
    # Step 2.a: Derivatives of $e^{4\phi}$
    phidD = ixp.zerorank1()
    phidDD = ixp.zerorank2()
    cf_dD  = ixp.declarerank1("cf_dD")
    cf_dDD = ixp.declarerank2("cf_dDD","sym01")
    if par.parval_from_str("EvolvedConformalFactor_cf") == "phi":
        for i in range(DIM):
            phidD[i]  = cf_dD[i]
            for j in range(DIM):
                phidDD[i][j] = cf_dDD[i][j]
    elif par.parval_from_str("EvolvedConformalFactor_cf") == "chi":
        for i in range(DIM):
            phidD[i]  = -sp.Rational(1,4)*exp4phi*cf_dD[i]
            for j in range(DIM):
                phidDD[i][j] = sp.Rational(1,4)*( exp4phi**2*cf_dD[i]*cf_dD[j] - exp4phi*cf_dDD[i][j] )
    elif par.parval_from_str("EvolvedConformalFactor_cf") == "W":
        exp2phi = (1 / cf)
        for i in range(DIM):
            phidD[i]  = -sp.Rational(1,2)*exp2phi*cf_dD[i]
            for j in range(DIM):
                phidDD[i][j] = sp.Rational(1,2)*( exp4phi*cf_dD[i]*cf_dD[j] - exp2phi*cf_dDD[i][j] )
    else:
        print("Error EvolvedConformalFactor_cf type = \""+par.parval_from_str("EvolvedConformalFactor_cf")+"\" unknown.")
        sys.exit(1)
    exp4phidD  = ixp.zerorank1()
    exp4phidDD = ixp.zerorank2()
    for i in range(DIM):
        exp4phidD[i] = 4*exp4phi*phidD[i]
        for j in range(DIM):
            exp4phidDD[i][j] = 16*exp4phi*phidD[i]*phidD[j] + 4*exp4phi*phidDD[i][j]
    # Step 2.b: Derivatives of gammaDD, the ADM three-metric
    gammaDDdD = ixp.zerorank3()
    gammaDDdDD = ixp.zerorank4()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                gammaDDdD[i][j][k] = exp4phidD[k] * gammabarDD[i][j] + exp4phi * gammabarDD_dD[i][j][k]
                for l in range(DIM):
                    gammaDDdDD[i][j][k][l] = exp4phidDD[k][l] * gammabarDD[i][j] + \
                                             exp4phidD[k] * gammabarDD_dD[i][j][l] + \
                                             exp4phidD[l] * gammabarDD_dD[i][j][k] + \
                                             exp4phi * gammabarDD_dDD[i][j][k][l]
    # Step 2.c: 3-Christoffel symbols associated with ADM 3-metric gammaDD
    # Step 2.c.i: First compute the inverse 3-metric gammaUU:
    gammaUU, detgamma = ixp.symm_matrix_inverter3x3(gammaDD)
    GammaUDD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    GammaUDD[i][j][k] += sp.Rational(1,2)*gammaUU[i][l]* \
                                    (gammaDDdD[l][j][k] + gammaDDdD[l][k][j] - gammaDDdD[j][k][l])
                    
    # Step 3: Define ADM extrinsic curvature KDD and
    #         its first spatial derivatives KDDdD
    #         in terms of BSSN quantities
    KDD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            KDD[i][j] = exp4phi * AbarDD[i][j] + sp.Rational(1, 3) * gammaDD[i][j] * trK
    KDDdD = ixp.zerorank3()
    trK_dD = ixp.declarerank1("trK_dD")
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                KDDdD[i][j][k] = exp4phidD[k] * AbarDD[i][j] + exp4phi * AbarDD_dD[i][j][k] + \
                                 sp.Rational(1, 3) * (gammaDDdD[i][j][k] * trK + gammaDD[i][j] * trK_dD[k])

# This module sets up a standard initial data function used for
#   setting up SENR initial data at all gridpoints.
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
from outputC import *
def BSSN_ID_function_string(cf,hDD,lambdaU,aDD,trK,alpha,vetU,betU):
    rhss = [trK,alpha,cf]
    lhss = ["in_gfs[IDX4S(TRKGF,i0,i1,i2)]","in_gfs[IDX4S(ALPHAGF,i0,i1,i2)]","in_gfs[IDX4S(CFGF,i0,i1,i2)]"]
    for i in range(3):
        rhss.append(lambdaU[i])
        lhss.append("in_gfs[IDX4S(LAMBDAU"+str(i)+"GF,i0,i1,i2)]")
        rhss.append(vetU[i])
        lhss.append("in_gfs[IDX4S(VETU"+str(i)+"GF,i0,i1,i2)]")
        rhss.append(betU[i])
        lhss.append("in_gfs[IDX4S(BETU"+str(i)+"GF,i0,i1,i2)]")
        for j in range(i,3):
            rhss.append(hDD[i][j])
            lhss.append("in_gfs[IDX4S(HDD" + str(i) + str(j) + "GF,i0,i1,i2)]")
            rhss.append(aDD[i][j])
            lhss.append("in_gfs[IDX4S(ADD" + str(i) + str(j) + "GF,i0,i1,i2)]")
    # Sort the lhss list alphabetically, and rhss to match:
    lhss,rhss = [list(x) for x in zip(*sorted(zip(lhss, rhss), key=lambda pair: pair[0]))]
    body = outputC(rhss, lhss, filename="returnstring",
                   params="preindent=1,CSE_enable=True,outCverbose=False",  # outCverbose=False to prevent
                                                                            # enormous output files.
                   prestring="", poststring="")
    desc = "Set up the initial data at all points on the numerical grid."
    add_to_Cfunction_dict(
        desc    =desc,
        name    ="initial_data",
        params  ="const paramstruct *restrict params,REAL *restrict xx[3], REAL *restrict in_gfs",
        body    =body,
        loopopts="AllPoints,Read_xxs")

# As documented in the NRPy+ tutorial module
#   Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb,
#   this module will construct the right-hand sides (RHSs)
#   expressions of the BSSN time evolution equations.
#
# Time-evolution equations for the BSSN gauge conditions are
#   specified in the BSSN_gauge_RHSs module and documented in
#   the Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb
#   NRPy+ tutorial module.
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step 1.a: import all needed modules from NRPy+:
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par    # NRPy+: Parameter interface
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
import sys                        # Standard Python modules for multiplatform OS-level functions
import BSSN.BSSN_quantities as Bq # NRPy+: This module depends on the parameter EvolvedConformalFactor_cf,
                                  #        which is defined in BSSN.BSSN_quantities
have_already_called_BSSN_RHSs_function = False
# Step 1.b: Set the coordinate system for the numerical grid:
#  DO NOT SET IN STANDALONE PYTHON MODULE
# par.set_parval_from_str("reference_metric::CoordSystem","Spherical")
def BSSN_RHSs():
    # Step 1.c: Given the chosen coordinate system, set up
    #           corresponding reference metric and needed
    #           reference metric quantities
    # The following function call sets up the reference metric
    #    and related quantities, including rescaling matrices ReDD,
    #    ReU, and hatted quantities.
    rfm.reference_metric()
    global have_already_called_BSSN_RHSs_function # setting to global enables other modules to see updated value.
    have_already_called_BSSN_RHSs_function = True
    # Step 1.d: Set spatial dimension (must be 3 for BSSN, as BSSN is
    #           a 3+1-dimensional decomposition of the general
    #           relativistic field equations)
    DIM = 3
    # Step 1.e: Import all basic (unrescaled) BSSN scalars & tensors
    import BSSN.BSSN_quantities as Bq
    Bq.BSSN_basic_tensors()
    gammabarDD = Bq.gammabarDD
    AbarDD     = Bq.AbarDD
    LambdabarU = Bq.LambdabarU
    trK   = Bq.trK
    alpha = Bq.alpha
    betaU = Bq.betaU
    # Step 1.f: Import all neeeded rescaled BSSN tensors:
    aDD = Bq.aDD
    cf  = Bq.cf
    lambdaU = Bq.lambdaU
    # Step 2.a.i: Import derivative expressions for betaU defined in the BSSN.BSSN_quantities module:
    Bq.betaU_derivs()
    betaU_dD = Bq.betaU_dD
    betaU_dDD = Bq.betaU_dDD
    # Step 2.a.ii: Import derivative expression for gammabarDD
    Bq.gammabar__inverse_and_derivs()
    gammabarDD_dupD = Bq.gammabarDD_dupD
    # Step 2.a.iii: First term of \partial_t \bar{\gamma}_{i j} right-hand side:
    # \beta^k \bar{\gamma}_{ij,k} + \beta^k_{,i} \bar{\gamma}_{kj} + \beta^k_{,j} \bar{\gamma}_{ik}
    gammabar_rhsDD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                gammabar_rhsDD[i][j] += betaU[k] * gammabarDD_dupD[i][j][k] + betaU_dD[k][i] * gammabarDD[k][j] \
                                        + betaU_dD[k][j] * gammabarDD[i][k]
    # Step 2.b.i: First import \bar{A}_{ij} = AbarDD[i][j], and its contraction trAbar = \bar{A}^k_k
    #           from BSSN.BSSN_quantities
    Bq.AbarUU_AbarUD_trAbar_AbarDD_dD()
    trAbar = Bq.trAbar
    # Step 2.b.ii: Import detgammabar quantities from BSSN.BSSN_quantities:
    Bq.detgammabar_and_derivs()
    detgammabar = Bq.detgammabar
    detgammabar_dD = Bq.detgammabar_dD
    # Step 2.b.ii: Compute the contraction \bar{D}_k \beta^k = \beta^k_{,k} + \frac{\beta^k \bar{\gamma}_{,k}}{2 \bar{\gamma}}
    Dbarbetacontraction = sp.sympify(0)
    for k in range(DIM):
        Dbarbetacontraction += betaU_dD[k][k] + betaU[k] * detgammabar_dD[k] / (2 * detgammabar)
    # Step 2.b.iii: Second term of \partial_t \bar{\gamma}_{i j} right-hand side:
    # \frac{2}{3} \bar{\gamma}_{i j} \left (\alpha \bar{A}_{k}^{k} - \bar{D}_{k} \beta^{k}\right )
    for i in range(DIM):
        for j in range(DIM):
            gammabar_rhsDD[i][j] += sp.Rational(2, 3) * gammabarDD[i][j] * (alpha * trAbar - Dbarbetacontraction)
    # Step 2.c: Third term of \partial_t \bar{\gamma}_{i j} right-hand side:
    # -2 \alpha \bar{A}_{ij}
    for i in range(DIM):
        for j in range(DIM):
            gammabar_rhsDD[i][j] += -2 * alpha * AbarDD[i][j]
    # Step 3.a: First term of \partial_t \bar{A}_{i j}:
    # \beta^k \partial_k \bar{A}_{ij} + \partial_i \beta^k \bar{A}_{kj} + \partial_j \beta^k \bar{A}_{ik}
    # First define AbarDD_dupD:
    AbarDD_dupD = Bq.AbarDD_dupD # From Bq.AbarUU_AbarUD_trAbar_AbarDD_dD()
    Abar_rhsDD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                Abar_rhsDD[i][j] += betaU[k] * AbarDD_dupD[i][j][k] + betaU_dD[k][i] * AbarDD[k][j] \
                                    + betaU_dD[k][j] * AbarDD[i][k]
    # Step 3.b: Second term of \partial_t \bar{A}_{i j}:
    # - (2/3) \bar{A}_{i j} \bar{D}_{k} \beta^{k} - 2 \alpha \bar{A}_{i k} {\bar{A}^{k}}_{j} + \alpha \bar{A}_{i j} K
    gammabarUU = Bq.gammabarUU  # From Bq.gammabar__inverse_and_derivs()
    AbarUD = Bq.AbarUD  # From Bq.AbarUU_AbarUD_trAbar()
    for i in range(DIM):
        for j in range(DIM):
            Abar_rhsDD[i][j] += -sp.Rational(2, 3) * AbarDD[i][j] * Dbarbetacontraction + alpha * AbarDD[i][j] * trK
            for k in range(DIM):
                Abar_rhsDD[i][j] += -2 * alpha * AbarDD[i][k] * AbarUD[k][j]
    # Step 3.c.i: Define partial derivatives of \phi in terms of evolved quantity "cf":
    Bq.phi_and_derivs()
    phi_dD = Bq.phi_dD
    phi_dupD = Bq.phi_dupD
    phi_dDD = Bq.phi_dDD
    exp_m4phi = Bq.exp_m4phi
    phi_dBarD = Bq.phi_dBarD  # phi_dBarD = Dbar_i phi = phi_dD (since phi is a scalar)
    phi_dBarDD = Bq.phi_dBarDD  # phi_dBarDD = Dbar_i Dbar_j phi (covariant derivative)
    # Step 3.c.ii: Define RbarDD
    Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()
    RbarDD = Bq.RbarDD
    # Step 3.c.iii: Define first and second derivatives of \alpha, as well as
    #         \bar{D}_i \bar{D}_j \alpha, which is defined just like phi
    alpha_dD = ixp.declarerank1("alpha_dD")
    alpha_dDD = ixp.declarerank2("alpha_dDD", "sym01")
    alpha_dBarD = alpha_dD
    alpha_dBarDD = ixp.zerorank2()
    GammabarUDD = Bq.GammabarUDD  # Defined in Bq.gammabar__inverse_and_derivs()
    for i in range(DIM):
        for j in range(DIM):
            alpha_dBarDD[i][j] = alpha_dDD[i][j]
            for k in range(DIM):
                alpha_dBarDD[i][j] += - GammabarUDD[k][i][j] * alpha_dD[k]
    # Step 3.c.iv: Define the terms in curly braces:
    curlybrackettermsDD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            curlybrackettermsDD[i][j] = -2 * alpha * phi_dBarDD[i][j] + 4 * alpha * phi_dBarD[i] * phi_dBarD[j] \
                                        + 2 * alpha_dBarD[i] * phi_dBarD[j] \
                                        + 2 * alpha_dBarD[j] * phi_dBarD[i] \
                                        - alpha_dBarDD[i][j] + alpha * RbarDD[i][j]
    # Step 3.c.v: Compute the trace:
    curlybracketterms_trace = sp.sympify(0)
    for i in range(DIM):
        for j in range(DIM):
            curlybracketterms_trace += gammabarUU[i][j] * curlybrackettermsDD[i][j]
    # Step 3.c.vi: Third and final term of Abar_rhsDD[i][j]:
    for i in range(DIM):
        for j in range(DIM):
            Abar_rhsDD[i][j] += exp_m4phi * (curlybrackettermsDD[i][j] -
                                             sp.Rational(1, 3) * gammabarDD[i][j] * curlybracketterms_trace)
    # Step 4: Right-hand side of conformal factor variable "cf". Supported
    #          options include: cf=phi, cf=W=e^(-2*phi) (default), and cf=chi=e^(-4*phi)
    # \partial_t phi = \left[\beta^k \partial_k \phi \right] <- TERM 1
    #                  + \frac{1}{6} \left (\bar{D}_{k} \beta^{k} - \alpha K \right ) <- TERM 2
    global cf_rhs
    cf_rhs = sp.Rational(1, 6) * (Dbarbetacontraction - alpha * trK)  # Term 2
    for k in range(DIM):
        cf_rhs += betaU[k] * phi_dupD[k]  # Term 1
    # Next multiply to convert phi_rhs to cf_rhs.
    if par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") == "phi":
        pass  # do nothing; cf_rhs = phi_rhs
    elif par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") == "W":
        cf_rhs *= -2 * cf  # cf_rhs = -2*cf*phi_rhs
    elif par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") == "chi":
        cf_rhs *= -4 * cf  # cf_rhs = -4*cf*phi_rhs
    else:
        print("Error: EvolvedConformalFactor_cf == " +
              par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") + " unsupported!")
        sys.exit(1)
    # Step 5: right-hand side of trK (trace of extrinsic curvature):
    # \partial_t K = \beta^k \partial_k K <- TERM 1
    #           + \frac{1}{3} \alpha K^{2} <- TERM 2
    #           + \alpha \bar{A}_{i j} \bar{A}^{i j} <- TERM 3
    #           - - e^{-4 \phi} (\bar{D}_{i} \bar{D}^{i} \alpha + 2 \bar{D}^{i} \alpha \bar{D}_{i} \phi ) <- TERM 4
    global trK_rhs
    # TERM 2:
    trK_rhs = sp.Rational(1, 3) * alpha * trK * trK
    trK_dupD = ixp.declarerank1("trK_dupD")
    for i in range(DIM):
        # TERM 1:
        trK_rhs += betaU[i] * trK_dupD[i]
    for i in range(DIM):
        for j in range(DIM):
            # TERM 4:
            trK_rhs += -exp_m4phi * gammabarUU[i][j] * (alpha_dBarDD[i][j] + 2 * alpha_dBarD[j] * phi_dBarD[i])
    AbarUU = Bq.AbarUU  # From Bq.AbarUU_AbarUD_trAbar()
    for i in range(DIM):
        for j in range(DIM):
            # TERM 3:
            trK_rhs += alpha * AbarDD[i][j] * AbarUU[i][j]
    # Step 6: right-hand side of \partial_t \bar{\Lambda}^i:
    # \partial_t \bar{\Lambda}^i = \beta^k \partial_k \bar{\Lambda}^i - \partial_k \beta^i \bar{\Lambda}^k <- TERM 1
    #                            + \bar{\gamma}^{j k} \hat{D}_{j} \hat{D}_{k} \beta^{i} <- TERM 2
    #                            + \frac{2}{3} \Delta^{i} \bar{D}_{j} \beta^{j} <- TERM 3
    #                            + \frac{1}{3} \bar{D}^{i} \bar{D}_{j} \beta^{j} <- TERM 4
    #                            - 2 \bar{A}^{i j} (\partial_{j} \alpha - 6 \partial_{j} \phi) <- TERM 5
    #                            + 2 \alpha \bar{A}^{j k} \Delta_{j k}^{i} <- TERM 6
    #                            - \frac{4}{3} \alpha \bar{\gamma}^{i j} \partial_{j} K <- TERM 7
    # Step 6.a: Term 1 of \partial_t \bar{\Lambda}^i: \beta^k \partial_k \bar{\Lambda}^i - \partial_k \beta^i \bar{\Lambda}^k
    # First we declare \bar{\Lambda}^i and \bar{\Lambda}^i_{,j} in terms of \lambda^i and \lambda^i_{,j}
    global LambdabarU_dupD # Used on the RHS of the Gamma-driving shift conditions
    LambdabarU_dupD = ixp.zerorank2()
    lambdaU_dupD = ixp.declarerank2("lambdaU_dupD", "nosym")
    for i in range(DIM):
        for j in range(DIM):
            LambdabarU_dupD[i][j] = lambdaU_dupD[i][j] * rfm.ReU[i] + lambdaU[i] * rfm.ReUdD[i][j]
    global Lambdabar_rhsU # Used on the RHS of the Gamma-driving shift conditions
    Lambdabar_rhsU = ixp.zerorank1()
    for i in range(DIM):
        for k in range(DIM):
            Lambdabar_rhsU[i] += betaU[k] * LambdabarU_dupD[i][k] - betaU_dD[i][k] * LambdabarU[k]  # Term 1
    # Step 6.b: Term 2 of \partial_t \bar{\Lambda}^i = \bar{\gamma}^{jk} (Term 2a + Term 2b + Term 2c)
    # Term 2a: \bar{\gamma}^{jk} \beta^i_{,kj}
    Term2aUDD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                Term2aUDD[i][j][k] += betaU_dDD[i][k][j]
    # Term 2b: \hat{\Gamma}^i_{mk,j} \beta^m + \hat{\Gamma}^i_{mk} \beta^m_{,j}
    #          + \hat{\Gamma}^i_{dj}\beta^d_{,k} - \hat{\Gamma}^d_{kj} \beta^i_{,d}
    Term2bUDD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for m in range(DIM):
                    Term2bUDD[i][j][k] += rfm.GammahatUDDdD[i][m][k][j] * betaU[m] \
                                          + rfm.GammahatUDD[i][m][k] * betaU_dD[m][j] \
                                          + rfm.GammahatUDD[i][m][j] * betaU_dD[m][k] \
                                          - rfm.GammahatUDD[m][k][j] * betaU_dD[i][m]
    # Term 2c: \hat{\Gamma}^i_{dj}\hat{\Gamma}^d_{mk} \beta^m - \hat{\Gamma}^d_{kj} \hat{\Gamma}^i_{md} \beta^m
    Term2cUDD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for m in range(DIM):
                    for d in range(DIM):
                        Term2cUDD[i][j][k] += (rfm.GammahatUDD[i][d][j] * rfm.GammahatUDD[d][m][k] \
                                               - rfm.GammahatUDD[d][k][j] * rfm.GammahatUDD[i][m][d]) * betaU[m]
    Lambdabar_rhsUpieceU = ixp.zerorank1()
    # Put it all together to get Term 2:
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                Lambdabar_rhsU[i] += gammabarUU[j][k] * (Term2aUDD[i][j][k] + Term2bUDD[i][j][k] + Term2cUDD[i][j][k])
                Lambdabar_rhsUpieceU[i] += gammabarUU[j][k] * (
                            Term2aUDD[i][j][k] + Term2bUDD[i][j][k] + Term2cUDD[i][j][k])
    # Step 6.c: Term 3 of \partial_t \bar{\Lambda}^i:
    #    \frac{2}{3} \Delta^{i} \bar{D}_{j} \beta^{j}
    DGammaU = Bq.DGammaU  # From Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()
    for i in range(DIM):
        Lambdabar_rhsU[i] += sp.Rational(2, 3) * DGammaU[i] * Dbarbetacontraction  # Term 3
    # Step 6.d: Term 4 of \partial_t \bar{\Lambda}^i:
    #           \frac{1}{3} \bar{D}^{i} \bar{D}_{j} \beta^{j}
    detgammabar_dDD = Bq.detgammabar_dDD  # From Bq.detgammabar_and_derivs()
    Dbarbetacontraction_dBarD = ixp.zerorank1()
    for k in range(DIM):
        for m in range(DIM):
            Dbarbetacontraction_dBarD[m] += betaU_dDD[k][k][m] + \
                                            (betaU_dD[k][m] * detgammabar_dD[k] +
                                             betaU[k] * detgammabar_dDD[k][m]) / (2 * detgammabar) \
                                            - betaU[k] * detgammabar_dD[k] * detgammabar_dD[m] / (
                                                        2 * detgammabar * detgammabar)
    for i in range(DIM):
        for m in range(DIM):
            Lambdabar_rhsU[i] += sp.Rational(1, 3) * gammabarUU[i][m] * Dbarbetacontraction_dBarD[m]
    # Step 6.e: Term 5 of \partial_t \bar{\Lambda}^i:
    #           - 2 \bar{A}^{i j} (\partial_{j} \alpha - 6 \alpha \partial_{j} \phi)
    for i in range(DIM):
        for j in range(DIM):
            Lambdabar_rhsU[i] += -2 * AbarUU[i][j] * (alpha_dD[j] - 6 * alpha * phi_dD[j])
    # Step 6.f: Term 6 of \partial_t \bar{\Lambda}^i:
    #           2 \alpha \bar{A}^{j k} \Delta^{i}_{j k}
    DGammaUDD = Bq.DGammaUDD  # From RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                Lambdabar_rhsU[i] += 2 * alpha * AbarUU[j][k] * DGammaUDD[i][j][k]
    # Step 6.g: Term 7 of \partial_t \bar{\Lambda}^i:
    #           -\frac{4}{3} \alpha \bar{\gamma}^{i j} \partial_{j} K
    trK_dD = ixp.declarerank1("trK_dD")
    for i in range(DIM):
        for j in range(DIM):
            Lambdabar_rhsU[i] += -sp.Rational(4, 3) * alpha * gammabarUU[i][j] * trK_dD[j]
    # Step 7: Rescale the RHS quantities so that the evolved
    #         variables are smooth across coord singularities
    global h_rhsDD,a_rhsDD,lambda_rhsU
    h_rhsDD = ixp.zerorank2()
    a_rhsDD = ixp.zerorank2()
    lambda_rhsU = ixp.zerorank1()
    for i in range(DIM):
        lambda_rhsU[i] = Lambdabar_rhsU[i] / rfm.ReU[i]
        for j in range(DIM):
            h_rhsDD[i][j] = gammabar_rhsDD[i][j] / rfm.ReDD[i][j]
            a_rhsDD[i][j] = Abar_rhsDD[i][j] / rfm.ReDD[i][j]
    # print(str(Abar_rhsDD[2][2]).replace("**","^").replace("_","").replace("xx","x").replace("sin(x2)","Sin[x2]").replace("sin(2*x2)","Sin[2*x2]").replace("cos(x2)","Cos[x2]").replace("detgbaroverdetghat","detg"))
    # print(str(Dbarbetacontraction).replace("**","^").replace("_","").replace("xx","x").replace("sin(x2)","Sin[x2]").replace("detgbaroverdetghat","detg"))
    # print(betaU_dD)
    # print(str(trK_rhs).replace("xx2","xx3").replace("xx1","xx2").replace("xx0","xx1").replace("**","^").replace("_","").replace("sin(xx2)","Sinx2").replace("xx","x").replace("sin(2*x2)","Sin2x2").replace("cos(x2)","Cosx2").replace("detgbaroverdetghat","detg"))
    # print(str(bet_rhsU[0]).replace("xx2","xx3").replace("xx1","xx2").replace("xx0","xx1").replace("**","^").replace("_","").replace("sin(xx2)","Sinx2").replace("xx","x").replace("sin(2*x2)","Sin2x2").replace("cos(x2)","Cosx2").replace("detgbaroverdetghat","detg"))

# This module registers rescaled BSSN T^{mu nu} source term variables
#    as AUX (i.e., not evolved) gridfunctions
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step P1: import all needed modules from NRPy+:
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import grid as gri                # NRPy+: Functions having to do with numerical grids
def define_BSSN_T4UUmunu_rescaled_source_terms():
    # Step 0: First check to see if this function has already been called.
    #   If so, do not register the gridfunctions again!
    for i in range(len(gri.glb_gridfcs_list)):
        if "sDD00" in gri.glb_gridfcs_list[i].name:
            return
    # Step 1: Declare as globals all quantities declared in this function.
    global rho,S,sD,sDD
    # Step 2: Register all needed *evolved* gridfunctions.
    # Step 2a: Register indexed quantities, using ixp.register_... functions
    sDD = ixp.register_gridfunctions_for_single_rank2("AUX", "sDD", "sym01")
    sD  = ixp.register_gridfunctions_for_single_rank1("AUX", "sD")
    # Step 2b: Register scalar quantities, using gri.register_gridfunctions()
    rho, S = gri.register_gridfunctions("AUX",["rho","S"])

# As documented in the NRPy+ tutorial module
#   Tutorial-BSSN_constraints.ipynb,
#   this module will construct expressions for the
#   BSSN Hamiltonian and momentum constraint equations.
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step 1: Initialize needed Python/NRPy+ modules
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par    # NRPy+: Parameter interface
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
import BSSN.BSSN_quantities as Bq # NRPy+: Computes useful BSSN quantities
import BSSN.BSSN_T4UUmunu_vars as BTmunu
thismodule = __name__
def BSSN_constraints(add_T4UUmunu_source_terms=False):
    # Step 1.a: Set spatial dimension (must be 3 for BSSN, as BSSN is
    #           a 3+1-dimensional decomposition of the general
    #           relativistic field equations)
    DIM = 3
    # Step 1.b: Given the chosen coordinate system, set up
    #           corresponding reference metric and needed
    #           reference metric quantities
    # The following function call sets up the reference metric
    #    and related quantities, including rescaling matrices ReDD,
    #    ReU, and hatted quantities.
    rfm.reference_metric()
    # Step 2: Hamiltonian constraint.
    # First declare all needed variables
    Bq.declare_BSSN_gridfunctions_if_not_declared_already()  # Sets trK
    Bq.BSSN_basic_tensors()  # Sets AbarDD
    Bq.gammabar__inverse_and_derivs()  # Sets gammabarUU
    Bq.AbarUU_AbarUD_trAbar_AbarDD_dD()  # Sets AbarUU and AbarDD_dD
    Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()  # Sets RbarDD
    Bq.phi_and_derivs()  # Sets phi_dBarD & phi_dBarDD
    ###############################
    ###############################
    #  HAMILTONIAN CONSTRAINT
    ###############################
    ###############################
    # Term 1: 2/3 K^2
    global H
    H = sp.Rational(2, 3) * Bq.trK ** 2
    # Term 2: -A_{ij} A^{ij}
    for i in range(DIM):
        for j in range(DIM):
            H += -Bq.AbarDD[i][j] * Bq.AbarUU[i][j]
    # Term 3a: trace(Rbar)
    Rbartrace = sp.sympify(0)
    for i in range(DIM):
        for j in range(DIM):
            Rbartrace += Bq.gammabarUU[i][j] * Bq.RbarDD[i][j]
    # Term 3b: -8 \bar{\gamma}^{ij} \bar{D}_i \phi \bar{D}_j \phi = -8*phi_dBar_times_phi_dBar
    # Term 3c: -8 \bar{\gamma}^{ij} \bar{D}_i \bar{D}_j \phi      = -8*phi_dBarDD_contraction
    phi_dBar_times_phi_dBar = sp.sympify(0)  # Term 3b
    phi_dBarDD_contraction = sp.sympify(0)  # Term 3c
    for i in range(DIM):
        for j in range(DIM):
            phi_dBar_times_phi_dBar += Bq.gammabarUU[i][j] * Bq.phi_dBarD[i] * Bq.phi_dBarD[j]
            phi_dBarDD_contraction += Bq.gammabarUU[i][j] * Bq.phi_dBarDD[i][j]
    # Add Term 3:
    H += Bq.exp_m4phi * (Rbartrace - 8 * (phi_dBar_times_phi_dBar + phi_dBarDD_contraction))
    if add_T4UUmunu_source_terms:
        M_PI = par.Cparameters("#define", thismodule, "M_PI","")  # M_PI is pi as defined in C
        BTmunu.define_BSSN_T4UUmunu_rescaled_source_terms()
        rho = BTmunu.rho
        H += -16*M_PI*rho
    # FIXME: ADD T4UUmunu SOURCE TERMS TO MOMENTUM CONSTRAINT!
    # Step 3: M^i, the momentum constraint
    ###############################
    ###############################
    #  MOMENTUM CONSTRAINT
    ###############################
    ###############################
    # SEE Tutorial-BSSN_constraints.ipynb for full documentation.
    global MU
    MU = ixp.zerorank1()
    # Term 2: 6 A^{ij} \partial_j \phi:
    for i in range(DIM):
        for j in range(DIM):
            MU[i] += 6*Bq.AbarUU[i][j]*Bq.phi_dD[j]
    # Term 3: -2/3 \bar{\gamma}^{ij} K_{,j}
    trK_dD = ixp.declarerank1("trK_dD") # Not defined in BSSN_RHSs; only trK_dupD is defined there.
    for i in range(DIM):
        for j in range(DIM):
            MU[i] += -sp.Rational(2,3)*Bq.gammabarUU[i][j]*trK_dD[j]
    # Next evaluate the conformal covariant derivative \bar{D}_j \bar{A}_{lm}
    AbarDD_dBarD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                AbarDD_dBarD[i][j][k] = Bq.AbarDD_dD[i][j][k]
                for l in range(DIM):
                    AbarDD_dBarD[i][j][k] += -Bq.GammabarUDD[l][k][i] * Bq.AbarDD[l][j]
                    AbarDD_dBarD[i][j][k] += -Bq.GammabarUDD[l][k][j] * Bq.AbarDD[i][l]
    # Term 1: Contract twice with the metric to make \bar{D}_{j} \bar{A}^{ij}
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    MU[i] += Bq.gammabarUU[i][k] * Bq.gammabarUU[j][l] * AbarDD_dBarD[k][l][j]
    # Finally, we multiply by e^{-4 phi} and rescale the momentum constraint:
    for i in range(DIM):
        MU[i] *= Bq.exp_m4phi / rfm.ReU[i]

# As documented in the NRPy+ tutorial module
#   Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb
#   this module will construct the right-hand sides (RHSs)
#   expressions for the time evolution equations of the
#   BSSN gauge quantities alpha and beta^i (i.e., the
#   lapse and shift)
#
# Non-gauge BSSN time-evolution equations are documented in the
#   NRPy+ tutorial module
#   Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb,
#   and separately constructed in the BSSN.BSSN_RHSs
#   Python module.
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step 1: Import all needed modules from NRPy+:
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par    # NRPy+: Parameter interface
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
import BSSN.BSSN_quantities as Bq # NRPy+: Computes useful BSSN quantities
import BSSN.BSSN_RHSs as Brhs     # NRPy+: Constructs BSSN right-hand-side expressions
import sys                        # Standard Python modules for multiplatform OS-level functions
# Step 1.a: Declare/initialize parameters for this module
thismodule = __name__
par.initialize_param(par.glb_param("char", thismodule, "LapseEvolutionOption", "OnePlusLog"))
par.initialize_param(par.glb_param("char", thismodule, "ShiftEvolutionOption", "GammaDriving2ndOrder_Covariant"))
def BSSN_gauge_RHSs():
    # Step 1.d: Set spatial dimension (must be 3 for BSSN, as BSSN is
    #           a 3+1-dimensional decomposition of the general
    #           relativistic field equations)
    DIM = 3
    # Step 1.e: Given the chosen coordinate system, set up
    #           corresponding reference metric and needed
    #           reference metric quantities
    # The following function call sets up the reference metric
    #    and related quantities, including rescaling matrices ReDD,
    #    ReU, and hatted quantities.
    rfm.reference_metric()
    # Step 1.f: Define needed BSSN quantities:
    # Declare scalars & tensors (in terms of rescaled BSSN quantities)
    Bq.BSSN_basic_tensors()
    Bq.betaU_derivs()
    # Declare BSSN_RHSs (excluding the time evolution equations for the gauge conditions),
    #    if they haven't already been declared.
    if Brhs.have_already_called_BSSN_RHSs_function == False:
        print("BSSN_gauge_RHSs() Error: You must call BSSN_RHSs() before calling BSSN_gauge_RHSs().")
        sys.exit(1)
    # Step 2: Lapse conditions
    LapseEvolOption = par.parval_from_str(thismodule + "::LapseEvolutionOption")
    # Step 2.a: The 1+log lapse condition:
    #   \partial_t \alpha = \beta^i \alpha_{,i} - 2*\alpha*K
    # First import expressions from BSSN_quantities
    cf = Bq.cf
    trK = Bq.trK
    alpha = Bq.alpha
    betaU = Bq.betaU
    # Implement the 1+log lapse condition
    global alpha_rhs
    alpha_rhs = sp.sympify(0)
    if LapseEvolOption == "OnePlusLog":
        alpha_rhs = -2 * alpha * trK
        alpha_dupD = ixp.declarerank1("alpha_dupD")
        for i in range(DIM):
            alpha_rhs += betaU[i] * alpha_dupD[i]
    # Step 2.b: Implement the harmonic slicing lapse condition
    elif LapseEvolOption == "HarmonicSlicing":
        if par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") == "W":
            alpha_rhs = -3 * cf ** (-4) * Brhs.cf_rhs
        elif par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") == "phi":
            alpha_rhs = 6 * sp.exp(6 * cf) * Brhs.cf_rhs
        else:
            print("Error LapseEvolutionOption==HarmonicSlicing unsupported for EvolvedConformalFactor_cf!=(W or phi)")
            sys.exit(1)
    # Step 2.c: Frozen lapse
    #    \partial_t \alpha = 0
    elif LapseEvolOption == "Frozen":
        alpha_rhs = sp.sympify(0)
        
    else:
        print("Error: "+thismodule + "::LapseEvolutionOption == "+LapseEvolOption+" not supported!")
        sys.exit(1)
    # Step 3.a: Set \partial_t \beta^i
    # First check that ShiftEvolutionOption parameter choice is supported.
    ShiftEvolOption = par.parval_from_str(thismodule + "::ShiftEvolutionOption")
    if ShiftEvolOption != "Frozen" and \
        ShiftEvolOption != "GammaDriving2ndOrder_NoCovariant" and \
        ShiftEvolOption != "GammaDriving2ndOrder_Covariant"  and \
        ShiftEvolOption != "GammaDriving2ndOrder_Covariant__Hatted" and \
        ShiftEvolOption != "GammaDriving1stOrder_Covariant" and \
        ShiftEvolOption != "GammaDriving1stOrder_Covariant__Hatted":
        print("Error: ShiftEvolutionOption == " + ShiftEvolOption + " unsupported!")
        sys.exit(1)
    # Next import expressions from BSSN_quantities
    BU = Bq.BU
    betU = Bq.betU
    betaU_dupD = Bq.betaU_dupD
    # Define needed quantities
    beta_rhsU = ixp.zerorank1()
    B_rhsU = ixp.zerorank1()
    # In the case of Frozen shift condition, we
    #    explicitly set the betaU and BU RHS's to zero
    #    instead of relying on the ixp.zerorank1()'s above,
    #    for safety.
    if ShiftEvolOption == "Frozen":
        for i in range(DIM):
            beta_rhsU[i] = sp.sympify(0)
            BU[i]        = sp.sympify(0)
    if ShiftEvolOption == "GammaDriving2ndOrder_NoCovariant":
        # Step 3.a.i: Compute right-hand side of beta^i
        # *  \partial_t \beta^i = \beta^j \beta^i_{,j} + B^i
        for i in range(DIM):
            beta_rhsU[i] += BU[i]
            for j in range(DIM):
                beta_rhsU[i] += betaU[j] * betaU_dupD[i][j]
        # Compute right-hand side of B^i:
        eta = par.Cparameters("REAL", thismodule, ["eta"],2.0)
        # Step 3.a.ii: Compute right-hand side of B^i
        # *  \partial_t B^i     = \beta^j B^i_{,j} + 3/4 * \partial_0 \Lambda^i - eta B^i
        # Step 3.a.iii: Define BU_dupD, in terms of derivative of rescaled variable \bet^i
        BU_dupD = ixp.zerorank2()
        betU_dupD = ixp.declarerank2("betU_dupD", "nosym")
        for i in range(DIM):
            for j in range(DIM):
                BU_dupD[i][j] = betU_dupD[i][j] * rfm.ReU[i] + betU[i] * rfm.ReUdD[i][j]
        # Step 3.a.iv: Compute \partial_0 \bar{\Lambda}^i = (\partial_t - \beta^i \partial_i) \bar{\Lambda}^j
        Lambdabar_partial0 = ixp.zerorank1()
        for i in range(DIM):
            Lambdabar_partial0[i] = Brhs.Lambdabar_rhsU[i]
        for i in range(DIM):
            for j in range(DIM):
                Lambdabar_partial0[j] += -betaU[i] * Brhs.LambdabarU_dupD[j][i]
        # Step 3.a.v: Evaluate RHS of B^i:
        for i in range(DIM):
            B_rhsU[i] += sp.Rational(3, 4) * Lambdabar_partial0[i] - eta * BU[i]
            for j in range(DIM):
                B_rhsU[i] += betaU[j] * BU_dupD[i][j]
    # Step 3.b: The right-hand side of the \partial_t \beta^i equation
    if "GammaDriving2ndOrder_Covariant" in ShiftEvolOption:
        # Step 3.b Option 2: \partial_t \beta^i = \left[\beta^j \bar{D}_j \beta^i\right] + B^{i}
        # First we need GammabarUDD, defined in Bq.gammabar__inverse_and_derivs()
        Bq.gammabar__inverse_and_derivs()
        ConnectionUDD = Bq.GammabarUDD
        # If instead we wish to use the Hatted covariant derivative, we replace
        #    ConnectionUDD with GammahatUDD:
        if ShiftEvolOption == "GammaDriving2ndOrder_Covariant__Hatted":
            ConnectionUDD = rfm.GammahatUDD
        # Then compute right-hand side:
        # Term 1: \beta^j \beta^i_{,j}
        for i in range(DIM):
            for j in range(DIM):
                beta_rhsU[i] += betaU[j] * betaU_dupD[i][j]
        # Term 2: \beta^j \bar{\Gamma}^i_{mj} \beta^m
        for i in range(DIM):
            for j in range(DIM):
                for m in range(DIM):
                    beta_rhsU[i] += betaU[j] * ConnectionUDD[i][m][j] * betaU[m]
        # Term 3: B^i
        for i in range(DIM):
            beta_rhsU[i] += BU[i]
    if "GammaDriving2ndOrder_Covariant" in ShiftEvolOption:
        ConnectionUDD = Bq.GammabarUDD
        # If instead we wish to use the Hatted covariant derivative, we replace
        #    ConnectionUDD with GammahatUDD:
        if ShiftEvolOption == "GammaDriving2ndOrder_Covariant__Hatted":
            ConnectionUDD = rfm.GammahatUDD
        # Step 3.c: Covariant option:
        #  \partial_t B^i = \beta^j \bar{D}_j B^i
        #               + \frac{3}{4} ( \partial_t \bar{\Lambda}^{i} - \beta^j \bar{D}_j \bar{\Lambda}^{i} )
        #               - \eta B^{i}
        #                 = \beta^j B^i_{,j} + \beta^j \bar{\Gamma}^i_{mj} B^m
        #               + \frac{3}{4}[ \partial_t \bar{\Lambda}^{i}
        #                            - \beta^j (\bar{\Lambda}^i_{,j} + \bar{\Gamma}^i_{mj} \bar{\Lambda}^m)]
        #               - \eta B^{i}
        # Term 1, part a: First compute B^i_{,j} using upwinded derivative
        BU_dupD = ixp.zerorank2()
        betU_dupD = ixp.declarerank2("betU_dupD", "nosym")
        for i in range(DIM):
            for j in range(DIM):
                BU_dupD[i][j] = betU_dupD[i][j] * rfm.ReU[i] + betU[i] * rfm.ReUdD[i][j]
        # Term 1: \beta^j B^i_{,j}
        for i in range(DIM):
            for j in range(DIM):
                B_rhsU[i] += betaU[j] * BU_dupD[i][j]
        # Term 2: \beta^j \bar{\Gamma}^i_{mj} B^m
        for i in range(DIM):
            for j in range(DIM):
                for m in range(DIM):
                    B_rhsU[i] += betaU[j] * ConnectionUDD[i][m][j] * BU[m]
        # Term 3: \frac{3}{4}\partial_t \bar{\Lambda}^{i}
        for i in range(DIM):
            B_rhsU[i] += sp.Rational(3, 4) * Brhs.Lambdabar_rhsU[i]
        # Term 4: -\frac{3}{4}\beta^j \bar{\Lambda}^i_{,j}
        for i in range(DIM):
            for j in range(DIM):
                B_rhsU[i] += -sp.Rational(3, 4) * betaU[j] * Brhs.LambdabarU_dupD[i][j]
        # Term 5: -\frac{3}{4}\beta^j \bar{\Gamma}^i_{mj} \bar{\Lambda}^m
        for i in range(DIM):
            for j in range(DIM):
                for m in range(DIM):
                    B_rhsU[i] += -sp.Rational(3, 4) * betaU[j] * ConnectionUDD[i][m][j] * Bq.LambdabarU[m]
        # Term 6: - \eta B^i
        # eta is a free parameter; we declare it here:
        eta = par.Cparameters("REAL", thismodule, ["eta"],2.0)
        for i in range(DIM):
            B_rhsU[i] += -eta * BU[i]
    if "GammaDriving1stOrder_Covariant" in ShiftEvolOption:
        # Step 3.c: \partial_t \beta^i = \left[\beta^j \bar{D}_j \beta^i\right] + 3/4 Lambdabar^i - eta*beta^i
        # First set \partial_t B^i = 0:
        B_rhsU = ixp.zerorank1()  # \partial_t B^i = 0
        # Second, set \partial_t beta^i RHS:
        # Compute covariant advection term:
        #  We need GammabarUDD, defined in Bq.gammabar__inverse_and_derivs()
        Bq.gammabar__inverse_and_derivs()
        ConnectionUDD = Bq.GammabarUDD
        # If instead we wish to use the Hatted covariant derivative, we replace
        #    ConnectionUDD with GammahatUDD:
        if ShiftEvolOption == "GammaDriving1stOrder_Covariant__Hatted":
            ConnectionUDD = rfm.GammahatUDD
        # Term 1: \beta^j \beta^i_{,j}
        for i in range(DIM):
            for j in range(DIM):
                beta_rhsU[i] += betaU[j] * betaU_dupD[i][j]
        # Term 2: \beta^j \bar{\Gamma}^i_{mj} \beta^m
        for i in range(DIM):
            for j in range(DIM):
                for m in range(DIM):
                    beta_rhsU[i] += betaU[j] * ConnectionUDD[i][m][j] * betaU[m]
        # Term 3: 3/4 Lambdabar^i - eta*beta^i
        eta = par.Cparameters("REAL", thismodule, ["eta"], 2.0)
        for i in range(DIM):
            beta_rhsU[i] += sp.Rational(3, 4) * Bq.LambdabarU[i] - eta * betaU[i]
    # Step 4: Rescale the BSSN gauge RHS quantities so that the evolved
    #         variables may remain smooth across coord singularities
    global vet_rhsU,bet_rhsU
    vet_rhsU = ixp.zerorank1()
    bet_rhsU = ixp.zerorank1()
    for i in range(DIM):
        vet_rhsU[i] = beta_rhsU[i] / rfm.ReU[i]
        bet_rhsU[i] = B_rhsU[i] / rfm.ReU[i]
    # print(str(Abar_rhsDD[2][2]).replace("**","^").replace("_","").replace("xx","x").replace("sin(x2)","Sin[x2]").replace("sin(2*x2)","Sin[2*x2]").replace("cos(x2)","Cos[x2]").replace("detgbaroverdetghat","detg"))
    # print(str(Dbarbetacontraction).replace("**","^").replace("_","").replace("xx","x").replace("sin(x2)","Sin[x2]").replace("detgbaroverdetghat","detg"))
    # print(betaU_dD)
    # print(str(trK_rhs).replace("xx2","xx3").replace("xx1","xx2").replace("xx0","xx1").replace("**","^").replace("_","").replace("sin(xx2)","Sinx2").replace("xx","x").replace("sin(2*x2)","Sin2x2").replace("cos(x2)","Cosx2").replace("detgbaroverdetghat","detg"))
    # print(str(bet_rhsU[0]).replace("xx2","xx3").replace("xx1","xx2").replace("xx0","xx1").replace("**","^").replace("_","").replace("sin(xx2)","Sinx2").replace("xx","x").replace("sin(2*x2)","Sin2x2").replace("cos(x2)","Cosx2").replace("detgbaroverdetghat","detg"))

# As documented in the NRPy+ tutorial module
#   Tutorial-BSSN_in_terms_of_ADM.ipynb,
#   this module will construct expressions for BSSN
#   quantities in terms of ADM quantities.
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step 1: Import needed core NRPy+ modules
from outputC import *             # NRPy+: Core C code output module
import NRPy_param_funcs as par    # NRPy+: Parameter interface
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
import sys                        # Standard Python modules for multiplatform OS-level functions
import BSSN.BSSN_quantities as Bq # NRPy+: This module depends on the parameter EvolvedConformalFactor_cf,
                                  #        which is defined in BSSN.BSSN_quantities
# Step 1.a: Set DIM=3, as we're using a 3+1 decomposition of Einstein's equations
DIM=3
# Step 2: All ADM quantities were input into this function in the Spherical or Cartesian
#         basis, as functions of r,th,ph or x,y,z, respectively. In Steps 1 and 2 above,
#         we converted them to the xx0,xx1,xx2 basis, and as functions of xx0,xx1,xx2.
#         Here we convert ADM quantities to their BSSN Curvilinear counterparts:
# Step 2.a: Convert ADM $\gamma_{ij}$ to BSSN $\bar{gamma}_{ij}$:
#           We have (Eqs. 2 and 3 of [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):
def gammabarDD_hDD(gammaDD):
    global gammabarDD,hDD
    if gammaDD == None:
        gammaDD = ixp.declarerank2("gammaDD","sym01")
    if rfm.have_already_called_reference_metric_function == False:
        print("BSSN.BSSN_in_terms_of_ADM.hDD_given_ADM(): Must call reference_metric() first!")
        sys.exit(1)
    # \bar{gamma}_{ij} = (\frac{\bar{gamma}}{gamma})^{1/3}*gamma_{ij}.
    gammaUU, gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)
    gammabarDD = ixp.zerorank2()
    hDD        = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            gammabarDD[i][j] = (rfm.detgammahat/gammaDET)**(sp.Rational(1,3))*gammaDD[i][j]
            hDD[i][j] = (gammabarDD[i][j] - rfm.ghatDD[i][j]) / rfm.ReDD[i][j]
# Step 2.b: Convert the extrinsic curvature K_{ij} to the trace-free extrinsic
#           curvature \bar{A}_{ij}, plus the trace of the extrinsic curvature K,
#           where (Eq. 3 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf)):
def trK_AbarDD_aDD(gammaDD,KDD):
    global trK,AbarDD,aDD
    if gammaDD == None:
        gammaDD = ixp.declarerank2("gammaDD","sym01")
    if KDD == None:
        KDD = ixp.declarerank2("KDD","sym01")
    if rfm.have_already_called_reference_metric_function == False:
        print("BSSN.BSSN_in_terms_of_ADM.trK_AbarDD(): Must call reference_metric() first!")
        sys.exit(1)
    # \bar{gamma}_{ij} = (\frac{\bar{gamma}}{gamma})^{1/3}*gamma_{ij}.
    gammaUU, gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)
    # K = gamma^{ij} K_{ij}, and
    # \bar{A}_{ij} &= (\frac{\bar{gamma}}{gamma})^{1/3}*(K_{ij} - \frac{1}{3}*gamma_{ij}*K)
    trK = sp.sympify(0)
    for i in range(DIM):
        for j in range(DIM):
            trK += gammaUU[i][j]*KDD[i][j]
    AbarDD = ixp.zerorank2()
    aDD    = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            AbarDD[i][j] = (rfm.detgammahat/gammaDET)**(sp.Rational(1,3))*(KDD[i][j] - sp.Rational(1,3)*gammaDD[i][j]*trK)
            aDD[i][j]    = AbarDD[i][j] / rfm.ReDD[i][j]
# Step 2.c: Define \bar{Lambda}^i (Eqs. 4 and 5 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf)):
def LambdabarU_lambdaU__exact_gammaDD(gammaDD):
    global LambdabarU, lambdaU
    if gammaDD == None:
        gammaDD = ixp.declarerank2("gammaDD","sym01")
    # \bar{Lambda}^i = \bar{gamma}^{jk}(\bar{Gamma}^i_{jk} - \hat{Gamma}^i_{jk}).
    gammabarDD_hDD(gammaDD)
    gammabarUU, gammabarDET = ixp.symm_matrix_inverter3x3(gammabarDD)
    # First compute Christoffel symbols \bar{Gamma}^i_{jk}, with respect to barred metric:
    GammabarUDD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    GammabarUDD[i][j][k] += sp.Rational(1, 2) * gammabarUU[i][l] * (
                                sp.diff(gammabarDD[l][j], rfm.xx[k]) +
                                sp.diff(gammabarDD[l][k], rfm.xx[j]) -
                                sp.diff(gammabarDD[j][k], rfm.xx[l]))
    # Next evaluate \bar{Lambda}^i, based on GammabarUDD above and GammahatUDD
    #       (from the reference metric):
    LambdabarU = ixp.zerorank1()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                LambdabarU[i] += gammabarUU[j][k] * (GammabarUDD[i][j][k] - rfm.GammahatUDD[i][j][k])
    lambdaU    = ixp.zerorank1()
    for i in range(DIM):
        lambdaU[i] = LambdabarU[i] / rfm.ReU[i]
# Step 2.d: Set the conformal factor variable cf, which is set
#           by the "BSSN_quantities::EvolvedConformalFactor_cf" parameter. For example if
#           "EvolvedConformalFactor_cf" is set to "phi", we can use Eq. 3 of
#           [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf),
#           which in arbitrary coordinates is written:
def cf_from_gammaDD(gammaDD):
    global cf
    if gammaDD == None:
        gammaDD = ixp.declarerank2("gammaDD","sym01")
    # \bar{Lambda}^i = \bar{gamma}^{jk}(\bar{Gamma}^i_{jk} - \hat{Gamma}^i_{jk}).
    gammabarDD_hDD(gammaDD)
    gammabarUU, gammabarDET = ixp.symm_matrix_inverter3x3(gammabarDD)
    gammaUU, gammaDET       = ixp.symm_matrix_inverter3x3(gammaDD)
    cf = sp.sympify(0)
    if par.parval_from_str("EvolvedConformalFactor_cf") == "phi":
        # phi = \frac{1}{12} log(\frac{gamma}{\bar{gamma}}).
        cf = sp.Rational(1, 12) * sp.log(gammaDET / gammabarDET)
    elif par.parval_from_str("EvolvedConformalFactor_cf") == "chi":
        # chi = exp(-4*phi) = exp(-4*\frac{1}{12}*(\frac{gamma}{\bar{gamma}}))
        #      = exp(-\frac{1}{3}*log(\frac{gamma}{\bar{gamma}})) = (\frac{gamma}{\bar{gamma}})^{-1/3}.
        #
        cf = (gammaDET / gammabarDET) ** (-sp.Rational(1, 3))
    elif par.parval_from_str("EvolvedConformalFactor_cf") == "W":
        # W = exp(-2*phi) = exp(-2*\frac{1}{12}*log(\frac{gamma}{\bar{gamma}}))
        #   = exp(-\frac{1}{6}*log(\frac{gamma}{\bar{gamma}})) = (\frac{gamma}{bar{gamma}})^{-1/6}.
        cf = (gammaDET / gammabarDET) ** (-sp.Rational(1, 6))
    else:
        print("Error EvolvedConformalFactor_cf type = \"" + par.parval_from_str(
            "EvolvedConformalFactor_cf") + "\" unknown.")
        sys.exit(1)
# Step 2.e: Rescale beta^i and B^i according to the prescription described in
#         the [BSSN in curvilinear coordinates tutorial module](Tutorial-BSSNCurvilinear.ipynb)
#         (also [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):
#
# \mathcal{V}^i &= beta^i/(ReU[i])
# \mathcal{B}^i &= B^i/(ReU[i])
def betU_vetU(betaU,BU):
    global vetU,betU
    if betaU == None:
        betaU = ixp.declarerank1("betaU")
    if BU == None:
        BU = ixp.declarerank1("BU")
    if rfm.have_already_called_reference_metric_function == False:
        print("BSSN.BSSN_in_terms_of_ADM.bet_vet(): Must call reference_metric() first!")
        sys.exit(1)
    vetU    = ixp.zerorank1()
    betU    = ixp.zerorank1()
    for i in range(DIM):
        vetU[i]    =      betaU[i] / rfm.ReU[i]
        betU[i]    =         BU[i] / rfm.ReU[i]

# This module provides functions that declare and define useful BSSN quantities
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step 1: Import all needed modules from NRPy+:
import NRPy_param_funcs as par    # NRPy+: Parameter interface
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import grid as gri                # NRPy+: Functions having to do with numerical grids
import reference_metric as rfm    # NRPy+: Reference metric support
import sys                        # Standard Python modules for multiplatform OS-level functions
# Step 1.a: Set the coordinate system for the numerical grid
#  DO NOT SET IN STANDALONE PYTHON MODULE
# par.set_parval_from_str("reference_metric::CoordSystem","Spherical")
# Step 1.b: Given the chosen coordinate system, set up 
#           corresponding reference metric and needed
#           reference metric quantities
# The following function call sets up the reference metric
#    and related quantities, including rescaling matrices ReDD,
#    ReU, and hatted quantities.
#  DO NOT CALL IN STANDALONE PYTHON MODULE
# rfm.reference_metric()
# Step 1.c: Set spatial dimension (must be 3 for BSSN, as BSSN is 
#           a 3+1-dimensional decomposition of the general 
#           relativistic field equations)
#  DO NOT CALL IN STANDALONE PYTHON MODULE
# DIM = 3
# par.set_parval_from_str("grid::DIM",DIM)
# Step 1.d: Declare/initialize parameters for this module
thismodule = __name__
par.initialize_param(par.glb_param("char", thismodule, "EvolvedConformalFactor_cf", "W"))
par.initialize_param(par.glb_param("bool", thismodule, "detgbarOverdetghat_equals_one", "True"))
par.initialize_param(par.glb_param("bool", thismodule, "LeaveRicciSymbolic", "False"))
def declare_BSSN_gridfunctions_if_not_declared_already():
    # Step 2: Register all needed BSSN gridfunctions.
    # Declare as globals all variables that may be
    # used outside this function
    global hDD,aDD,lambdaU,vetU,betU,trK,cf,alpha
    #   Check to see if this function has already been called.
    #   If so, do not register the gridfunctions again!
    for i in range(len(gri.glb_gridfcs_list)):
        if "hDD00" in gri.glb_gridfcs_list[i].name:
            hDD = ixp.declarerank2("hDD", "sym01")
            aDD = ixp.declarerank2("aDD", "sym01")
            lambdaU = ixp.declarerank1("lambdaU")
            vetU = ixp.declarerank1("vetU")
            betU = ixp.declarerank1("betU")
            trK, cf, alpha = sp.symbols('trK cf alpha', real=True)
            return hDD, aDD, lambdaU, vetU, betU, trK, cf, alpha
    # Step 2.a: Register indexed quantities, using ixp.register_... functions
    hDD = ixp.register_gridfunctions_for_single_rank2("EVOL", "hDD", "sym01")
    aDD = ixp.register_gridfunctions_for_single_rank2("EVOL", "aDD", "sym01")
    lambdaU = ixp.register_gridfunctions_for_single_rank1("EVOL", "lambdaU")
    vetU = ixp.register_gridfunctions_for_single_rank1("EVOL", "vetU")
    betU = ixp.register_gridfunctions_for_single_rank1("EVOL", "betU")
    # Step 2.b: Register scalar quantities, using gri.register_gridfunctions()
    trK, cf, alpha = gri.register_gridfunctions("EVOL", ["trK", "cf", "alpha"])
    return hDD, aDD, lambdaU, vetU, betU, trK, cf, alpha
# Step 3: Define all basic conformal BSSN tensors
#        gammabarDD,AbarDD,LambdabarU,betaU,BU
#        in terms of BSSN gridfunctions.
def BSSN_basic_tensors():
    # Step 3.a: Declare as globals all variables that may be used
    #           outside this function, declare BSSN gridfunctions
    #           if not defined already, and set DIM=3.
    global gammabarDD,AbarDD,LambdabarU,betaU,BU
    hDD, aDD, lambdaU, vetU, betU, trK, cf, alpha = declare_BSSN_gridfunctions_if_not_declared_already()
    DIM = 3
    # Step 3.a.i: gammabarDD and AbarDD:
    gammabarDD = ixp.zerorank2()
    AbarDD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            # gammabar_{ij}  = h_{ij}*ReDD[i][j] + gammahat_{ij}
            gammabarDD[i][j] = hDD[i][j] * rfm.ReDD[i][j] + rfm.ghatDD[i][j]
            # Abar_{ij}      = a_{ij}*ReDD[i][j]
            AbarDD[i][j] = aDD[i][j] * rfm.ReDD[i][j]
    # Step 3.a.ii: LambdabarU, betaU, and BU:
    LambdabarU = ixp.zerorank1()
    betaU = ixp.zerorank1()
    BU = ixp.zerorank1()
    for i in range(DIM):
        LambdabarU[i] = lambdaU[i] * rfm.ReU[i]
        betaU[i] = vetU[i] * rfm.ReU[i]
        BU[i] = betU[i] * rfm.ReU[i]
# Step 4: gammabarUU and spatial derivatives of gammabarDD,
#         including GammabarUDD
def gammabar__inverse_and_derivs():
    # Step 4.a: Declare as globals all expressions that may be used
    #           outside this function, declare BSSN gridfunctions
    #           if not defined already, and set DIM=3.
    global gammabarUU, gammabarDD_dD, gammabarDD_dupD, gammabarDD_dDD, GammabarUDD
    hDD, aDD, lambdaU, vetU, betU, trK, cf, alpha = declare_BSSN_gridfunctions_if_not_declared_already()
    DIM = 3
    # This function needs gammabarDD, defined in BSSN_basic_tensors()
    BSSN_basic_tensors()
    # Step 4.a.i: gammabarUU:
    gammabarUU, dummydet = ixp.symm_matrix_inverter3x3(gammabarDD)
    # Step 4.b.i: gammabarDDdD[i][j][k]
    #               = \hat{\gamma}_{ij,k} + h_{ij,k} \text{ReDD[i][j]} + h_{ij} \text{ReDDdD[i][j][k]}.
    gammabarDD_dD = ixp.zerorank3()
    gammabarDD_dupD = ixp.zerorank3()
    hDD_dD = ixp.declarerank3("hDD_dD", "sym01")
    hDD_dupD = ixp.declarerank3("hDD_dupD", "sym01")
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                gammabarDD_dD[i][j][k] = rfm.ghatDDdD[i][j][k] + \
                                         hDD_dD[i][j][k] * rfm.ReDD[i][j] + hDD[i][j] * rfm.ReDDdD[i][j][k]
                # Compute associated upwinded derivative, needed for the \bar{\gamma}_{ij} RHS
                gammabarDD_dupD[i][j][k] = rfm.ghatDDdD[i][j][k] + \
                                           hDD_dupD[i][j][k] * rfm.ReDD[i][j] + hDD[i][j] * rfm.ReDDdD[i][j][k]
    # Step 4.b.ii: Compute gammabarDD_dDD in terms of the rescaled BSSN quantity hDD
    #      and its derivatives, as well as the reference metric and rescaling
    #      matrix, and its derivatives (expression given below):
    hDD_dDD = ixp.declarerank4("hDD_dDD", "sym01_sym23")
    gammabarDD_dDD = ixp.zerorank4()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    # gammabar_{ij,kl} = gammahat_{ij,kl}
                    #                  + h_{ij,kl} ReDD[i][j]
                    #                  + h_{ij,k} ReDDdD[i][j][l] + h_{ij,l} ReDDdD[i][j][k]
                    #                  + h_{ij} ReDDdDD[i][j][k][l]
                    gammabarDD_dDD[i][j][k][l] = rfm.ghatDDdDD[i][j][k][l]
                    gammabarDD_dDD[i][j][k][l] += hDD_dDD[i][j][k][l] * rfm.ReDD[i][j]
                    gammabarDD_dDD[i][j][k][l] += hDD_dD[i][j][k] * rfm.ReDDdD[i][j][l] + \
                                                  hDD_dD[i][j][l] * rfm.ReDDdD[i][j][k]
                    gammabarDD_dDD[i][j][k][l] += hDD[i][j] * rfm.ReDDdDD[i][j][k][l]
    # Step 4.b.iii: Define barred Christoffel symbol \bar{\Gamma}^{i}_{kl} = GammabarUDD[i][k][l] (see expression below)
    GammabarUDD = ixp.zerorank3()
    for i in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                for m in range(DIM):
                    # Gammabar^i_{kl} = 1/2 * gammabar^{im} ( gammabar_{mk,l} + gammabar_{ml,k} - gammabar_{kl,m}):
                    GammabarUDD[i][k][l] += sp.Rational(1, 2) * gammabarUU[i][m] * \
                                            (gammabarDD_dD[m][k][l] + gammabarDD_dD[m][l][k] - gammabarDD_dD[k][l][m])
# Step 5: det(gammabarDD) and its derivatives
def detgammabar_and_derivs():
    # Step 5.a: Declare as globals all expressions that may be used
    #           outside this function, declare BSSN gridfunctions
    #           if not defined already, and set DIM=3.
    global detgammabar,detgammabar_dD,detgammabar_dDD
    hDD, aDD, lambdaU, vetU, betU, trK, cf, alpha = declare_BSSN_gridfunctions_if_not_declared_already()
    DIM = 3
    detgbarOverdetghat = sp.sympify(1)
    detgbarOverdetghat_dD = ixp.zerorank1()
    detgbarOverdetghat_dDD = ixp.zerorank2()
    if par.parval_from_str(thismodule + "::detgbarOverdetghat_equals_one") == "False":
        print("Error: detgbarOverdetghat_equals_one=\"False\" is not fully implemented yet.")
        sys.exit(1)
    ## Approach for implementing detgbarOverdetghat_equals_one=False:
    #     detgbarOverdetghat = gri.register_gridfunctions("AUX", ["detgbarOverdetghat"])
    #     detgbarOverdetghatInitial = gri.register_gridfunctions("AUX", ["detgbarOverdetghatInitial"])
    #     detgbarOverdetghat_dD = ixp.declarerank1("detgbarOverdetghat_dD")
    #     detgbarOverdetghat_dDD = ixp.declarerank2("detgbarOverdetghat_dDD", "sym01")
    # Step 5.b: Define detgammabar, detgammabar_dD, and detgammabar_dDD (needed for \partial_t \bar{\Lambda}^i below)
    detgammabar = detgbarOverdetghat * rfm.detgammahat
    detgammabar_dD = ixp.zerorank1()
    for i in range(DIM):
        detgammabar_dD[i] = detgbarOverdetghat_dD[i] * rfm.detgammahat + detgbarOverdetghat * rfm.detgammahatdD[i]
    detgammabar_dDD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            detgammabar_dDD[i][j] = detgbarOverdetghat_dDD[i][j] * rfm.detgammahat + \
                                    detgbarOverdetghat_dD[i] * rfm.detgammahatdD[j] + \
                                    detgbarOverdetghat_dD[j] * rfm.detgammahatdD[i] + \
                                    detgbarOverdetghat * rfm.detgammahatdDD[i][j]
# Step 6: Quantities related to conformal traceless
#         extrinsic curvature AbarDD:
#         AbarUU, AbarUD, and trAbar
def AbarUU_AbarUD_trAbar_AbarDD_dD():
    # Step 6.a: Declare as globals all expressions that may be used
    #           outside this function, declare BSSN gridfunctions
    #           if not defined already, and set DIM=3.
    global AbarUU,AbarUD,trAbar,AbarDD_dD,AbarDD_dupD
    hDD, aDD, lambdaU, vetU, betU, trK, cf, alpha = declare_BSSN_gridfunctions_if_not_declared_already()
    DIM = 3
    # Define AbarDD and gammabarDD in terms of BSSN gridfunctions
    BSSN_basic_tensors()
    # Define gammabarUU in terms of BSSN gridfunctions
    gammabar__inverse_and_derivs()
    # Step 6.a.i: Compute Abar^{ij} in terms of Abar_{ij} and gammabar^{ij}
    AbarUU = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    # Abar^{ij} = gammabar^{ik} gammabar^{jl} Abar_{kl}
                    AbarUU[i][j] += gammabarUU[i][k] * gammabarUU[j][l] * AbarDD[k][l]
    # Step 6.a.ii: Compute Abar^i_j in terms of Abar_{ij} and gammabar^{ij}
    AbarUD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                # Abar^i_j = gammabar^{ik} Abar_{kj}
                AbarUD[i][j] += gammabarUU[i][k] * AbarDD[k][j]
    # Step 6.a.iii: Compute Abar^k_k = trace of Abar:
    trAbar = sp.sympify(0)
    for k in range(DIM):
        for j in range(DIM):
            # Abar^k_k = gammabar^{kj} Abar_{jk}
            trAbar += gammabarUU[k][j] * AbarDD[j][k]
            
    # Step 6.a.iv: Compute Abar_{ij,k}
    AbarDD_dD = ixp.zerorank3()
    AbarDD_dupD = ixp.zerorank3()
    aDD_dD   = ixp.declarerank3("aDD_dD"  ,"sym01")
    aDD_dupD = ixp.declarerank3("aDD_dupD","sym01")
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                AbarDD_dupD[i][j][k] = rfm.ReDDdD[i][j][k]*aDD[i][j] + rfm.ReDD[i][j]*aDD_dupD[i][j][k]
                AbarDD_dD[i][j][k]   = rfm.ReDDdD[i][j][k]*aDD[i][j] + rfm.ReDD[i][j]*aDD_dD[  i][j][k]
# Step 7: The conformal ("barred") Ricci tensor RbarDD
#         and associated quantities
def RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU():
    # Step 7.a: Declare as globals all expressions that may be used
    #           outside this function, declare BSSN gridfunctions
    #           if not defined already, and set DIM=3.
    global RbarDD,DGammaUDD,gammabarDD_dHatD,DGammaU
    hDD, aDD, lambdaU, vetU, betU, trK, cf, alpha = declare_BSSN_gridfunctions_if_not_declared_already()
    DIM = 3
    # GammabarUDD is used below, defined in
    #    gammabar__inverse_and_derivs()
    gammabar__inverse_and_derivs()
    # Step 7.a.i: Define \varepsilon_{ij} = epsDD[i][j]
    epsDD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            epsDD[i][j] = hDD[i][j] * rfm.ReDD[i][j]
    # Step 7.a.ii: Define epsDD_dD[i][j][k]
    hDD_dD = ixp.declarerank3("hDD_dD", "sym01")
    epsDD_dD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                epsDD_dD[i][j][k] = hDD_dD[i][j][k] * rfm.ReDD[i][j] + hDD[i][j] * rfm.ReDDdD[i][j][k]
    # Step 7.a.iii: Define epsDD_dDD[i][j][k][l]
    hDD_dDD = ixp.declarerank4("hDD_dDD", "sym01_sym23")
    epsDD_dDD = ixp.zerorank4()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    epsDD_dDD[i][j][k][l] = hDD_dDD[i][j][k][l] * rfm.ReDD[i][j] + \
                                            hDD_dD[i][j][k] * rfm.ReDDdD[i][j][l] + \
                                            hDD_dD[i][j][l] * rfm.ReDDdD[i][j][k] + \
                                            hDD[i][j] * rfm.ReDDdDD[i][j][k][l]
    # Step 7.a.iv: DhatgammabarDDdD[i][j][l] = \bar{\gamma}_{ij;\hat{l}}
    # \bar{\gamma}_{ij;\hat{l}} = \varepsilon_{i j,l}
    #                           - \hat{\Gamma}^m_{i l} \varepsilon_{m j}
    #                           - \hat{\Gamma}^m_{j l} \varepsilon_{i m}
    gammabarDD_dHatD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for l in range(DIM):
                gammabarDD_dHatD[i][j][l] = epsDD_dD[i][j][l]
                for m in range(DIM):
                    gammabarDD_dHatD[i][j][l] += - rfm.GammahatUDD[m][i][l] * epsDD[m][j] \
                                                 - rfm.GammahatUDD[m][j][l] * epsDD[i][m]
    # Step 7.a.v: \bar{\gamma}_{ij;\hat{l},k} = DhatgammabarDD_dHatD_dD[i][j][l][k]:
    #        \bar{\gamma}_{ij;\hat{l},k} = \varepsilon_{ij,lk}
    #                                      - \hat{\Gamma}^m_{i l,k} \varepsilon_{m j}
    #                                      - \hat{\Gamma}^m_{i l} \varepsilon_{m j,k}
    #                                      - \hat{\Gamma}^m_{j l,k} \varepsilon_{i m}
    #                                      - \hat{\Gamma}^m_{j l} \varepsilon_{i m,k}
    gammabarDD_dHatD_dD = ixp.zerorank4()
    for i in range(DIM):
        for j in range(DIM):
            for l in range(DIM):
                for k in range(DIM):
                    gammabarDD_dHatD_dD[i][j][l][k] = epsDD_dDD[i][j][l][k]
                    for m in range(DIM):
                        gammabarDD_dHatD_dD[i][j][l][k] += -rfm.GammahatUDDdD[m][i][l][k] * epsDD[m][j] \
                                                           - rfm.GammahatUDD[m][i][l] * epsDD_dD[m][j][k] \
                                                           - rfm.GammahatUDDdD[m][j][l][k] * epsDD[i][m] \
                                                           - rfm.GammahatUDD[m][j][l] * epsDD_dD[i][m][k]
    # Step 7.a.vi: \bar{\gamma}_{ij;\hat{l}\hat{k}} = DhatgammabarDD_dHatDD[i][j][l][k]
    #          \bar{\gamma}_{ij;\hat{l}\hat{k}} = \partial_k \hat{D}_{l} \varepsilon_{i j}
    #                                           - \hat{\Gamma}^m_{lk} \left(\hat{D}_{m} \varepsilon_{i j}\right)
    #                                           - \hat{\Gamma}^m_{ik} \left(\hat{D}_{l} \varepsilon_{m j}\right)
    #                                           - \hat{\Gamma}^m_{jk} \left(\hat{D}_{l} \varepsilon_{i m}\right)
    gammabarDD_dHatDD = ixp.zerorank4()
    for i in range(DIM):
        for j in range(DIM):
            for l in range(DIM):
                for k in range(DIM):
                    gammabarDD_dHatDD[i][j][l][k] = gammabarDD_dHatD_dD[i][j][l][k]
                    for m in range(DIM):
                        gammabarDD_dHatDD[i][j][l][k] += - rfm.GammahatUDD[m][l][k] * gammabarDD_dHatD[i][j][m] \
                                                         - rfm.GammahatUDD[m][i][k] * gammabarDD_dHatD[m][j][l] \
                                                         - rfm.GammahatUDD[m][j][k] * gammabarDD_dHatD[i][m][l]
    # Step 7.b: Second term of RhatDD: compute \hat{D}_{j} \bar{\Lambda}^{k} = LambarU_dHatD[k][j]
    lambdaU_dD = ixp.declarerank2("lambdaU_dD", "nosym")
    LambarU_dHatD = ixp.zerorank2()
    for j in range(DIM):
        for k in range(DIM):
            LambarU_dHatD[k][j] = lambdaU_dD[k][j] * rfm.ReU[k] + lambdaU[k] * rfm.ReUdD[k][j]
            for m in range(DIM):
                LambarU_dHatD[k][j] += rfm.GammahatUDD[k][m][j] * lambdaU[m] * rfm.ReU[m]
                
    # Step 7.c: Conformal Ricci tensor, part 3: The \Delta^{k} \Delta_{(i j) k}  
    #           + \bar{\gamma}^{k l}*(2 \Delta_{k(i}^{m} \Delta_{j) m l} 
    #           + \Delta_{i k}^{m} \Delta_{m j l}) terms
    # Step 7.c.i: Define \Delta^i_{jk} = \bar{\Gamma}^i_{jk} - \hat{\Gamma}^i_{jk} = DGammaUDD[i][j][k]
    DGammaUDD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                DGammaUDD[i][j][k] = GammabarUDD[i][j][k] - rfm.GammahatUDD[i][j][k]
    # Step 7.c.ii: Define \Delta^i = \bar{\gamma}^{jk} \Delta^i_{jk}
    DGammaU = ixp.zerorank1()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                DGammaU[i] += gammabarUU[j][k] * DGammaUDD[i][j][k]
    # Step 7.c.iii: Define \Delta_{ijk} = \bar{\gamma}_{im} \Delta^m_{jk}
    DGammaDDD = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for m in range(DIM):
                    DGammaDDD[i][j][k] += gammabarDD[i][m] * DGammaUDD[m][j][k]
    if par.parval_from_str(thismodule+"::LeaveRicciSymbolic") == "True":
        for i in range(len(gri.glb_gridfcs_list)):
            if "RbarDD00" in gri.glb_gridfcs_list[i].name:
                return
        RbarDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","RbarDD","sym01")
        return
    # Step 7.d: Summing the terms and defining \bar{R}_{ij}
    # Step 7.d.i: Add the first term to RbarDD:
    #         Rbar_{ij} += - \frac{1}{2} \bar{\gamma}^{k l} \hat{D}_{k} \hat{D}_{l} \bar{\gamma}_{i j}
    RbarDD = ixp.zerorank2()
    RbarDDpiece = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    RbarDD[i][j] += -sp.Rational(1, 2) * gammabarUU[k][l] * gammabarDD_dHatDD[i][j][l][k]
                    RbarDDpiece[i][j] += -sp.Rational(1, 2) * gammabarUU[k][l] * gammabarDD_dHatDD[i][j][l][k]
    # Step 7.d.ii: Add the second term to RbarDD:
    #         Rbar_{ij} += (1/2) * (gammabar_{ki} Lambar^k_{;\hat{j}} + gammabar_{kj} Lambar^k_{;\hat{i}})
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                RbarDD[i][j] += sp.Rational(1, 2) * ( gammabarDD[k][i] * LambarU_dHatD[k][j] +
                                                      gammabarDD[k][j] * LambarU_dHatD[k][i]  )
    # Step 7.d.iii: Add the remaining term to RbarDD:
    #      Rbar_{ij} += \Delta^{k} \Delta_{(i j) k} = 1/2 \Delta^{k} (\Delta_{i j k} + \Delta_{j i k})
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                RbarDD[i][j] += sp.Rational(1, 2) * DGammaU[k] * (DGammaDDD[i][j][k] + DGammaDDD[j][i][k])
    # Step 7.d.iv: Add the final term to RbarDD:
    #      Rbar_{ij} += \bar{\gamma}^{k l} (\Delta^{m}_{k i} \Delta_{j m l}
    #                   + \Delta^{m}_{k j} \Delta_{i m l}
    #                   + \Delta^{m}_{i k} \Delta_{m j l})
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    for m in range(DIM):
                        RbarDD[i][j] += gammabarUU[k][l] * (DGammaUDD[m][k][i] * DGammaDDD[j][m][l] +
                                                            DGammaUDD[m][k][j] * DGammaDDD[i][m][l] +
                                                            DGammaUDD[m][i][k] * DGammaDDD[m][j][l])
# Step 8: The unrescaled shift vector betaU spatial derivatives:
#         betaUdD & betaUdDD, written in terms of the
#         rescaled shift vector vetU
def betaU_derivs():
    # Step 8.i: Declare as globals all expressions that may be used
    #           outside this function, declare BSSN gridfunctions
    #           if not defined already, and set DIM=3.
    global betaU_dD,betaU_dupD,betaU_dDD
    hDD, aDD, lambdaU, vetU, betU, trK, cf, alpha = declare_BSSN_gridfunctions_if_not_declared_already()
    DIM = 3
    # Step 8.ii: Compute the unrescaled shift vector beta^i = ReU[i]*vet^i
    vetU_dD = ixp.declarerank2("vetU_dD", "nosym")
    vetU_dupD = ixp.declarerank2("vetU_dupD", "nosym")  # Needed for upwinded \beta^i_{,j}
    vetU_dDD = ixp.declarerank3("vetU_dDD", "sym12")  # Needed for \beta^i_{,j}
    betaU_dD = ixp.zerorank2()
    betaU_dupD = ixp.zerorank2()  # Needed for, e.g., \beta^i RHS
    betaU_dDD = ixp.zerorank3()  # Needed for, e.g., \bar{\Lambda}^i RHS
    for i in range(DIM):
        for j in range(DIM):
            betaU_dD[i][j] = vetU_dD[i][j] * rfm.ReU[i] + vetU[i] * rfm.ReUdD[i][j]
            betaU_dupD[i][j] = vetU_dupD[i][j] * rfm.ReU[i] + vetU[i] * rfm.ReUdD[i][j]  # Needed for \beta^i RHS
            for k in range(DIM):
                # Needed for, e.g., \bar{\Lambda}^i RHS:
                betaU_dDD[i][j][k] = vetU_dDD[i][j][k] * rfm.ReU[i] + vetU_dD[i][j] * rfm.ReUdD[i][k] + \
                                     vetU_dD[i][k] * rfm.ReUdD[i][j] + vetU[i] * rfm.ReUdDD[i][j][k]
# Step 9: Standard BSSN conformal factor phi,
#         and its partial and covariant derivatives,
#         all in terms of BSSN gridfunctions like cf
def phi_and_derivs():
    # Step 9.a: Declare as globals all expressions that may be used
    #           outside this function, declare BSSN gridfunctions
    #           if not defined already, and set DIM=3.
    global phi_dD,phi_dupD,phi_dDD,exp_m4phi,phi_dBarD,phi_dBarDD
    hDD, aDD, lambdaU, vetU, betU, trK, cf, alpha = declare_BSSN_gridfunctions_if_not_declared_already()
    DIM = 3
    # GammabarUDD is used below, defined in
    #    gammabar__inverse_and_derivs()
    gammabar__inverse_and_derivs()
    # Step 9.a.i: Define partial derivatives of \phi in terms of evolved quantity "cf":
    cf_dD = ixp.declarerank1("cf_dD")
    cf_dupD = ixp.declarerank1("cf_dupD")  # Needed for \partial_t \phi next.
    cf_dDD = ixp.declarerank2("cf_dDD", "sym01")
    phi_dD = ixp.zerorank1()
    phi_dupD = ixp.zerorank1()
    phi_dDD = ixp.zerorank2()
    exp_m4phi = sp.sympify(0)
    # Step 9.a.ii: Assuming cf=phi, define exp_m4phi, phi_dD,
    #              phi_dupD (upwind finite-difference version of phi_dD), and phi_DD
    if par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") == "phi":
        for i in range(DIM):
            phi_dD[i] = cf_dD[i]
            phi_dupD[i] = cf_dupD[i]
            for j in range(DIM):
                phi_dDD[i][j] = cf_dDD[i][j]
        exp_m4phi = sp.exp(-4 * cf)
    # Step 9.a.iii: Assuming cf=W=e^{-2 phi}, define exp_m4phi, phi_dD,
    #               phi_dupD (upwind finite-difference version of phi_dD), and phi_DD
    if par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") == "W":
        # \partial_i W = \partial_i (e^{-2 phi}) = -2 e^{-2 phi} \partial_i phi
        # -> \partial_i phi = -\partial_i cf / (2 cf)
        for i in range(DIM):
            phi_dD[i] = - cf_dD[i] / (2 * cf)
            phi_dupD[i] = - cf_dupD[i] / (2 * cf)
            for j in range(DIM):
                # \partial_j \partial_i phi = - \partial_j [\partial_i cf / (2 cf)]
                #                           = - cf_{,ij} / (2 cf) + \partial_i cf \partial_j cf / (2 cf^2)
                phi_dDD[i][j] = (- cf_dDD[i][j] + cf_dD[i] * cf_dD[j] / cf) / (2 * cf)
        exp_m4phi = cf * cf
    # Step 9.a.iv: Assuming cf=chi=e^{-4 phi}, define exp_m4phi, phi_dD,
    #              phi_dupD (upwind finite-difference version of phi_dD), and phi_DD
    if par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") == "chi":
        # \partial_i chi = \partial_i (e^{-4 phi}) = -4 e^{-4 phi} \partial_i phi
        # -> \partial_i phi = -\partial_i cf / (4 cf)
        for i in range(DIM):
            phi_dD[i] = - cf_dD[i] / (4 * cf)
            phi_dupD[i] = - cf_dupD[i] / (4 * cf)
            for j in range(DIM):
                # \partial_j \partial_i phi = - \partial_j [\partial_i cf / (4 cf)]
                #                           = - cf_{,ij} / (4 cf) + \partial_i cf \partial_j cf / (4 cf^2)
                phi_dDD[i][j] = (- cf_dDD[i][j] + cf_dD[i] * cf_dD[j] / cf) / (4 * cf)
        exp_m4phi = cf
    # Step 9.a.v: Error out if unsupported EvolvedConformalFactor_cf choice is made:
    cf_choice = par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf")
    if not (cf_choice == "phi" or cf_choice == "W" or cf_choice == "chi"):
        print("Error: EvolvedConformalFactor_cf == " + par.parval_from_str("BSSN.BSSN_quantities::EvolvedConformalFactor_cf") + " unsupported!")
        sys.exit(1)
    # Step 9.b: Define phi_dBarD = phi_dD (since phi is a scalar) and phi_dBarDD (covariant derivative)
    #          \bar{D}_i \bar{D}_j \phi = \phi_{;\bar{i}\bar{j}} = \bar{D}_i \phi_{,j}
    #                                   = \phi_{,ij} - \bar{\Gamma}^k_{ij} \phi_{,k}
    phi_dBarD = phi_dD
    phi_dBarDD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            phi_dBarDD[i][j] = phi_dDD[i][j]
            for k in range(DIM):
                phi_dBarDD[i][j] += - GammabarUDD[k][i][j] * phi_dD[k]

# As documented in the NRPy+ tutorial module
#   Tutorial-BSSN_stress_energy_source_terms.ipynb,
#   this module will construct expressions for
#   BSSN stress-energy source terms, in terms of
#   elements of T^{mu nu}.
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
import sympy as sp                         # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par             # NRPy+: Parameter interface
import indexedexp as ixp                   # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm             # NRPy+: Reference metric support
import BSSN.ADMBSSN_tofrom_4metric as AB4m # NRPy+: ADM/BSSN <-> 4-metric conversions
import sys                                 # Standard Python modules for multiplatform OS-level functions
thismodule = __name__
# Define BSSN source terms in terms of T^{mu nu} or T_{mu nu}
def stress_energy_source_terms_ito_T4UU_and_ADM_or_BSSN_metricvars(inputvars,custom_T4UU=None):
    # Step 1: Check if rfm.reference_metric() already called. If not, BSSN
    #         quantities are not yet defined, so cannot proceed!
    if rfm.have_already_called_reference_metric_function == False:
        print("BSSN_source_terms_ito_T4UU(): Must call reference_metric() first!")
        sys.exit(1)
    # Step 2.a: Define gamma4DD[mu][nu] = g_{mu nu} + n_{mu} n_{nu}
    alpha = sp.symbols("alpha", real=True)
    zero = sp.sympify(0)
    n4D = [-alpha, zero, zero, zero]
    AB4m.g4DD_ito_BSSN_or_ADM(inputvars)
    gamma4DD = ixp.zerorank2(DIM=4)
    for mu in range(4):
        for nu in range(4):
            gamma4DD[mu][nu] = AB4m.g4DD[mu][nu] + n4D[mu] * n4D[nu]
    # Step 2.b: If expression for components of T4UU not given, declare T4UU here
    if custom_T4UU == None:
        T4UU = ixp.declarerank2("T4UU","sym01",DIM=4)
    else:
        T4UU = custom_T4UU
    # Step 2.c: Define BSSN source terms
    global SDD,SD,S,rho
    # Step 2.c.i: S_{ij} = gamma_{i mu} gamma_{j nu} T^{mu nu}
    SDD = ixp.zerorank2()
    for i in range(3):
        for j in range(3):
            for mu in range(4):
                for nu in range(4):
                    SDD[i][j] += gamma4DD[i+1][mu] * gamma4DD[j+1][nu] * T4UU[mu][nu]
    # Step 2.c.ii: S_{i} = -gamma_{i mu} n_{nu} T^{mu nu}
    SD = ixp.zerorank1()
    for i in range(3):
        for mu in range(4):
            for nu in range(4):
                SD[i] += - gamma4DD[i+1][mu] * n4D[nu] * T4UU[mu][nu]
    # Step 2.c.iii: S = gamma^{ij} S_{ij}
    if inputvars == "ADM":
        gammaDD = ixp.declarerank2("gammaDD", "sym01")
        gammaUU, dummydet = ixp.symm_matrix_inverter3x3(gammaDD)  # Set gammaUU
    elif inputvars == "BSSN":
        import BSSN.ADM_in_terms_of_BSSN as AitoB  # NRPy+: ADM quantities in terms of BSSN quantities
        AitoB.ADM_in_terms_of_BSSN()
        gammaUU = AitoB.gammaUU
    S = zero
    for i in range(3):
        for j in range(3):
            S += gammaUU[i][j] * SDD[i][j]
    # Step 2.c.iv: rho = n_{mu} n_{nu} T^{mu nu}
    rho = zero
    for mu in range(4):
        for nu in range(4):
            rho += n4D[mu] * n4D[nu] * T4UU[mu][nu]
    return SDD,SD,S,rho
# Step 3: Add BSSN stress-energy source terms to BSSN RHSs
def BSSN_source_terms_for_BSSN_RHSs(custom_T4UU=None):
    global sourceterm_trK_rhs, sourceterm_a_rhsDD, sourceterm_lambda_rhsU, sourceterm_Lambdabar_rhsU
    # Step 3.a: Call BSSN_source_terms_ito_T4UU to get SDD, SD, S, & rho
    if custom_T4UU == "unrescaled BSSN source terms already given":
        SDD = ixp.declarerank2("SDD", "sym01")
        SD = ixp.declarerank1("SD")
        S = sp.symbols("S", real=True)
        rho = sp.symbols("rho", real=True)
    else:
        SDD,SD,S,rho = stress_energy_source_terms_ito_T4UU_and_ADM_or_BSSN_metricvars("BSSN", custom_T4UU)
    PI = par.Cparameters("REAL", thismodule, ["PI"], "3.14159265358979323846264338327950288")
    alpha = sp.symbols("alpha", real=True)
    # Step 3.b: trK_rhs
    sourceterm_trK_rhs = 4 * PI * alpha * (rho + S)
    # Step 3.c: Abar_rhsDD:
    # Step 3.c.i: Compute trace-free part of S_{ij}:
    import BSSN.BSSN_quantities as Bq
    Bq.BSSN_basic_tensors()  # Sets gammabarDD
    gammabarUU, dummydet = ixp.symm_matrix_inverter3x3(Bq.gammabarDD)  # Set gammabarUU
    tracefree_SDD = ixp.zerorank2()
    for i in range(3):
        for j in range(3):
            tracefree_SDD[i][j] = SDD[i][j]
    for i in range(3):
        for j in range(3):
            for k in range(3):
                for m in range(3):
                    tracefree_SDD[i][j] += -sp.Rational(1, 3) * Bq.gammabarDD[i][j] * gammabarUU[k][m] * SDD[k][m]
    # Step 3.c.ii: Define exp_m4phi = e^{-4 phi}
    Bq.phi_and_derivs()
    # Step 3.c.iii: Evaluate stress-energy part of AbarDD's RHS
    sourceterm_a_rhsDD = ixp.zerorank2()
    for i in range(3):
        for j in range(3):
            Abar_rhsDDij = -8 * PI * alpha * Bq.exp_m4phi * tracefree_SDD[i][j]
            sourceterm_a_rhsDD[i][j] = Abar_rhsDDij / rfm.ReDD[i][j]
    # Step 3.d: Stress-energy part of Lambdabar_rhsU = stressenergy_Lambdabar_rhsU
    sourceterm_Lambdabar_rhsU = ixp.zerorank1()
    for i in range(3):
        for j in range(3):
            sourceterm_Lambdabar_rhsU[i] += -16 * PI * alpha * gammabarUU[i][j] * SD[j]
    sourceterm_lambda_rhsU = ixp.zerorank1()
    for i in range(3):
        sourceterm_lambda_rhsU[i] = sourceterm_Lambdabar_rhsU[i] / rfm.ReU[i]
# Step 4: Add BSSN stress-energy source terms to BSSN constraints
def BSSN_source_terms_for_BSSN_constraints(custom_T4UU=None):
    global sourceterm_H, sourceterm_MU
    # Step 4.a: Call BSSN_source_terms_ito_T4UU to get SDD, SD, S, & rho
    if custom_T4UU == "unrescaled BSSN source terms already given":
        SDD = ixp.declarerank2("SDD", "sym01")
        SD = ixp.declarerank1("SD")
        S = sp.symbols("S", real=True)
        rho = sp.symbols("rho", real=True)
    else:
        SDD,SD,S,rho = stress_energy_source_terms_ito_T4UU_and_ADM_or_BSSN_metricvars("BSSN", custom_T4UU)
    PI = par.Cparameters("REAL", thismodule, ["PI"], "3.14159265358979323846264338327950288")
    # Step 4.b: Add source term to the Hamiltonian constraint H
    sourceterm_H = -16 * PI * rho
    # Step 4.c: Add source term to the momentum constraint M^i
    # Step 4.c.i: Compute gammaUU in terms of BSSN quantities
    import BSSN.ADM_in_terms_of_BSSN as AitoB
    AitoB.ADM_in_terms_of_BSSN()  # Provides gammaUU
    # Step 4.c.ii: Raise S_i
    SU = ixp.zerorank1()
    for i in range(3):
        for j in range(3):
            SU[i] += AitoB.gammaUU[i][j] * SD[j]
    # Step 4.c.iii: Add source term to momentum constraint & rescale:
    sourceterm_MU = ixp.zerorank1()
    for i in range(3):
        sourceterm_MU[i] = -8 * PI * SU[i] / rfm.ReU[i]

# This module sets up Brill-Lindquist initial data in terms of
# the variables used in BSSN_RHSs.py
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# # Setting up Two Black Hole Initial Data, in Curvilinear Coordinates
# 
# ## This module sets up initial data for two black holes at rest in spherical coordinates, converts the initial data/basis to the desired CoordSystem, and then rescales all quantities according to the BSSNCurvilinear prescription.
# 
# ### Brill-Lindquist initial data ([Brill & Lindquist, Phys. Rev. 131, 471, 1963](https://journals.aps.org/pr/abstract/10.1103/PhysRev.131.471); see also Eq. 1 of [Brandt & Br\"ugmann, arXiv:gr-qc/9711015v1](https://arxiv.org/pdf/gr-qc/9711015v1.pdf)) may be written in terms of the BSSN conformal factor and ADM extrinsic curvature as
# 
# $$\psi = e^{\phi} = 1 + \sum_{i=1}^N \frac{m_{(i)}}{2 \left|\vec{r}_{(i)} - \vec{r}\right|};\quad K_{ij}=0.$$
# 
# These data consist of $N$ nonspinning black holes initially at rest. This module restricts to the case of two such black holes, positioned along either the $x$ or $z$ axis. Here, we implement $N=2$.
# 
# **Inputs for $\psi$**:
# * The position and (bare) mass of black hole 1: $\left(x_{(1)},y_{(1)},z_{(1)}\right)$ and $m_{(1)}$, respectively
# * The position and (bare) mass of black hole 2: $\left(x_{(2)},y_{(2)},z_{(2)}\right)$ and $m_{(2)}$, respectively
# 
# **Additional variables needed for spacetime evolution**:
# * Desired coordinate system
# * Desired initial lapse $\alpha$ and shift $\beta^i$
# 
# **Transformation to curvilinear coordinates**:
# * Once the above variables have been set in Cartesian coordinates, we will apply the appropriate coordinate transformations and tensor rescalings ([described in the BSSN NRPy+ tutorial module](Tutorial-BSSNCurvilinear.ipynb))
# Step 1: Initialize core Python/NRPy+ modules
import sympy as sp             # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par # NRPy+: Parameter interface
import indexedexp as ixp       # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import BSSN.ADM_Exact_Spherical_or_Cartesian_to_BSSNCurvilinear as AtoB
thismodule = __name__
BH1_posn_x, BH1_posn_y, BH1_posn_z = par.Cparameters("REAL", thismodule,
                                                     ["BH1_posn_x", "BH1_posn_y", "BH1_posn_z"],
                                                     [0.0, 0.0, +0.5])
BH1_mass = par.Cparameters("REAL", thismodule, ["BH1_mass"], 1.0)
BH2_posn_x, BH2_posn_y, BH2_posn_z = par.Cparameters("REAL", thismodule,
                                                     ["BH2_posn_x", "BH2_posn_y", "BH2_posn_z"],
                                                     [0.0, 0.0, -0.5])
BH2_mass = par.Cparameters("REAL", thismodule, ["BH2_mass"], 1.0)
# ComputeADMGlobalsOnly == True will only set up the ADM global quantities.
#                       == False will perform the full ADM SphorCart->BSSN Curvi conversion
def BrillLindquist(ComputeADMGlobalsOnly = False):
    global Cartxyz,gammaCartDD, KCartDD, alphaCart, betaCartU, BCartU
    
    # Step 2: Setting up Brill-Lindquist initial data
    # Step 2.a: Set spatial dimension (must be 3 for BSSN)
    DIM = 3
    par.set_parval_from_str("grid::DIM",DIM)
    global Cartxyz, gammaCartDD, KCartDD, alphaCart, betaCartU, BCartU
    Cartxyz = ixp.declarerank1("Cartxyz")
    # Step 2.b: Set psi, the conformal factor:
    psi = sp.sympify(1)
    psi += BH1_mass / ( 2 * sp.sqrt((Cartxyz[0]-BH1_posn_x)**2 + (Cartxyz[1]-BH1_posn_y)**2 + (Cartxyz[2]-BH1_posn_z)**2) )
    psi += BH2_mass / ( 2 * sp.sqrt((Cartxyz[0]-BH2_posn_x)**2 + (Cartxyz[1]-BH2_posn_y)**2 + (Cartxyz[2]-BH2_posn_z)**2) )
    # Step 2.c: Set all needed ADM variables in Cartesian coordinates
    gammaCartDD = ixp.zerorank2()
    KCartDD     = ixp.zerorank2() # K_{ij} = 0 for these initial data
    for i in range(DIM):
        gammaCartDD[i][i] = psi**4
    alphaCart = 1/psi**2
    betaCartU = ixp.zerorank1() # We generally choose \beta^i = 0 for these initial data
    BCartU    = ixp.zerorank1() # We generally choose B^i = 0 for these initial data
    if ComputeADMGlobalsOnly == True:
        return
    
    cf,hDD,lambdaU,aDD,trK,alpha,vetU,betU = \
        AtoB.Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear("Cartesian",Cartxyz, 
                                                                    gammaCartDD,KCartDD,alphaCart,betaCartU,BCartU)
    import BSSN.BSSN_ID_function_string as bIDf
    global returnfunction
    returnfunction = bIDf.BSSN_ID_function_string(cf, hDD, lambdaU, aDD, trK, alpha, vetU, betU)

# This module implements the gammabar = gammahat constraint.
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step P1: import all needed modules from NRPy+:
from outputC import *             # NRPy+: Core C code output module
import finite_difference as fin   # NRPy+: Finite difference C code generation module
import grid as gri                # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import BSSN.BSSN_quantities as Bq # NRPy+: Computes useful BSSN quantities
import os                         # Standard Python modules for multiplatform OS-level functions
def Enforce_Detgammabar_Constraint_symb_expressions():
    # Set spatial dimension (must be 3 for BSSN)
    DIM = 3
    # Then we set the coordinate system for the numerical grid
    rfm.reference_metric()  # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.
    # We will need the h_{ij} quantities defined within BSSN_RHSs
    #    below when we enforce the gammahat=gammabar constraint
    # Step 1: All barred quantities are defined in terms of BSSN rescaled gridfunctions,
    #         which we declare here in case they haven't yet been declared elsewhere.
    Bq.declare_BSSN_gridfunctions_if_not_declared_already()
    hDD = Bq.hDD
    Bq.BSSN_basic_tensors()
    gammabarDD = Bq.gammabarDD
    # First define the Kronecker delta:
    KroneckerDeltaDD = ixp.zerorank2()
    for i in range(DIM):
        KroneckerDeltaDD[i][i] = sp.sympify(1)
    # The detgammabar in BSSN_RHSs is set to detgammahat when BSSN_RHSs::detgbarOverdetghat_equals_one=True (default),
    #    so we manually compute it here:
    dummygammabarUU, detgammabar = ixp.symm_matrix_inverter3x3(gammabarDD)
    # Next apply the constraint enforcement equation above.
    hprimeDD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            hprimeDD[i][j] = \
                (nrpyAbs(rfm.detgammahat) / detgammabar) ** (sp.Rational(1, 3)) * (KroneckerDeltaDD[i][j] + hDD[i][j]) \
                - KroneckerDeltaDD[i][j]
    enforce_detg_constraint_symb_expressions = [
        lhrh(lhs=gri.gfaccess("in_gfs", "hDD00"), rhs=hprimeDD[0][0]),
        lhrh(lhs=gri.gfaccess("in_gfs", "hDD01"), rhs=hprimeDD[0][1]),
        lhrh(lhs=gri.gfaccess("in_gfs", "hDD02"), rhs=hprimeDD[0][2]),
        lhrh(lhs=gri.gfaccess("in_gfs", "hDD11"), rhs=hprimeDD[1][1]),
        lhrh(lhs=gri.gfaccess("in_gfs", "hDD12"), rhs=hprimeDD[1][2]),
        lhrh(lhs=gri.gfaccess("in_gfs", "hDD22"), rhs=hprimeDD[2][2])]
    return enforce_detg_constraint_symb_expressions
def output_Enforce_Detgammabar_Constraint_Ccode(outdir="BSSN/", exprs="", Read_xxs=False):
    # Step 0: Check if outdir is string; error out if not.
    check_if_string__error_if_not(outdir,"outdir")
    desc = "Enforce det(gammabar) = det(gammahat) constraint."
    name = "enforce_detgammabar_constraint"
    params = "const rfm_struct *restrict rfmstruct,const paramstruct *restrict params, REAL *restrict in_gfs"
    loopopts = "AllPoints,Enable_rfm_precompute"
    if Read_xxs:
        params = "const paramstruct *restrict params, REAL *restrict xx[3], REAL *restrict in_gfs"
        loopopts = "AllPoints,Read_xxs"
    outCfunction(
        outfile=os.path.join(outdir, name + ".h"), desc=desc, name=name, params=params,
        body=fin.FD_outputC("returnstring", exprs,
                            params="outCverbose=False,preindent=1,includebraces=False").replace("IDX4", "IDX4S"),
        loopopts=loopopts)

// This function takes as input either (x,y,z) or (r,th,ph) and outputs
//   all ADM quantities in the Cartesian or Spherical basis, respectively.
void ID_ADM_SphorCart(const REAL xyz_or_rthph[3], 
                     REAL *gammaDD00,REAL *gammaDD01,REAL *gammaDD02,REAL *gammaDD11,REAL *gammaDD12,REAL *gammaDD22,
                     REAL *KDD00,REAL *KDD01,REAL *KDD02,REAL *KDD11,REAL *KDD12,REAL *KDD22,
                     REAL *alpha,
                     REAL *betaU0,REAL *betaU1,REAL *betaU2,
                     REAL *BU0,REAL *BU1,REAL *BU2) {
      const REAL r  = xyz_or_rthph[0];
      const REAL th = xyz_or_rthph[1];
      const REAL ph = xyz_or_rthph[2];
      const double tmp0 = (1.0/4.0)*M;
      const double tmp1 = pow(M, 2);
      const double tmp2 = pow(chi, 2)*tmp1;
      const double tmp3 = -tmp2;
      const double tmp4 = sqrt(tmp1 + tmp3);
      const double tmp5 = (1.0/4.0)*tmp4;
      const double tmp6 = cos(th);
      const double tmp7 = 1.0/r;
      const double tmp8 = (1.0/4.0)*tmp7*(M + tmp4) + 1;
      const double tmp9 = pow(r, 2)*pow(tmp8, 4);
      const double tmp10 = tmp2*pow(tmp6, 2) + tmp9;
      const double tmp11 = r*pow(tmp8, 2);
      const double tmp12 = -M + tmp11 + tmp4;
      const double tmp13 = tmp2 + tmp9;
      const double tmp14 = 2*tmp11;
      const double tmp15 = sin(th);
      const double tmp16 = pow(tmp15, 2);
      const double tmp17 = tmp16*tmp2;
      const double tmp18 = pow(tmp13, 2) - tmp17*(-M*tmp14 + tmp13);
      const double tmp19 = 1.0/tmp10;
      const double tmp20 = tmp16*tmp19;
      const double tmp21 = pow(tmp10*tmp18, -1.0/2.0);
      const double tmp22 = pow(M, 4);
      *gammaDD00 = tmp10*pow(r + tmp0 + tmp5, 2)/(pow(r, 3)*tmp12);
      *gammaDD01 = 0;
      *gammaDD02 = 0;
      *gammaDD11 = tmp10;
      *gammaDD12 = 0;
      *gammaDD22 = tmp18*tmp20;
      *KDD00 = 0;
      *KDD01 = 0;
      *KDD02 = chi*tmp1*tmp20*tmp21*tmp8*(-pow(chi, 4)*tmp22 + 3*pow(r, 4)*pow(tmp8, 8) - tmp17*(tmp3 + tmp9) + 2*tmp2*tmp9)/sqrt(r*tmp12);
      *KDD11 = 0;
      *KDD12 = -pow(chi, 3)*tmp14*pow(tmp15, 3)*tmp19*tmp21*tmp22*tmp6*sqrt(tmp12*tmp7)*(r - tmp0 - tmp5);
      *KDD22 = 0;
      *alpha = 1;
      *betaU0 = 0;
      *betaU1 = 0;
      *betaU2 = 0;
      *BU0 = 0;
      *BU1 = 0;
      *BU2 = 0;
}

# As documented in the NRPy+ tutorial module
#   Tutorial-Psi4.ipynb,
#   this module will construct a generic
#   expression for \psi_4
# Author: Zachariah B. Etienne
#         zachetie **at** gmail **dot* com
# Step 1.a: import all needed modules from NRPy+:
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
def Psi4(specify_tetrad=True):
    global psi4_im_pt, psi4_re_pt
    # Step 1.b: Given the chosen coordinate system, set up
    #           corresponding reference metric and needed
    #           reference metric quantities
    # The following function call sets up the reference metric
    #    and related quantities, including rescaling matrices ReDD,
    #    ReU, and hatted quantities.
    rfm.reference_metric()
    # Step 1.c: Set spatial dimension (must be 3 for BSSN, as BSSN is
    #           a 3+1-dimensional decomposition of the general
    #           relativistic field equations)
    DIM = 3
    # Step 1.d: Import all ADM quantities as written in terms of BSSN quantities
    import BSSN.ADM_in_terms_of_BSSN as AB
    AB.ADM_in_terms_of_BSSN()
    # Step 1.e: Set up tetrad vectors
    if specify_tetrad==True:
        import BSSN.Psi4_tetrads as BP4t
        BP4t.Psi4_tetrads()
        mre4U = BP4t.mre4U
        mim4U = BP4t.mim4U
        n4U   = BP4t.n4U
    else:
        # For code validation against NRPy+ psi_4 tutorial module (Tutorial-Psi4.ipynb);
        #   ensures a more complete code validation.
        mre4U = ixp.declarerank1("mre4U", DIM=4)
        mim4U = ixp.declarerank1("mim4U", DIM=4)
        n4U = ixp.declarerank1("n4U", DIM=4)
    # Step 2: Construct the (rank-4) Riemann curvature tensor associated with the ADM 3-metric:
    RDDDD = ixp.zerorank4()
    gammaDDdDD = AB.gammaDDdDD
    for i in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                for m in range(DIM):
                    RDDDD[i][k][l][m] = sp.Rational(1, 2) * \
                                        (gammaDDdDD[i][m][k][l] + gammaDDdDD[k][l][i][m] - gammaDDdDD[i][l][k][m] -
                                         gammaDDdDD[k][m][i][l])
    # ... then we add the term on the right:
    gammaDD = AB.gammaDD
    GammaUDD = AB.GammaUDD
    for i in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                for m in range(DIM):
                    for n in range(DIM):
                        for p in range(DIM):
                            RDDDD[i][k][l][m] += gammaDD[n][p] * \
                                                 (GammaUDD[n][k][l] * GammaUDD[p][i][m] - GammaUDD[n][k][m] * GammaUDD[p][i][l])
    # Step 3: Construct the (rank-4) tensor in term 1 of psi_4 (referring to Eq 5.1 in
    #   Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf
    rank4term1DDDD = ixp.zerorank4()
    KDD = AB.KDD
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    rank4term1DDDD[i][j][k][l] = RDDDD[i][j][k][l] + KDD[i][k] * KDD[l][j] - KDD[i][l] * KDD[k][j]
    # Step 4: Construct the (rank-3) tensor in term 2 of psi_4 (referring to Eq 5.1 in
    #   Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf
    rank3term2DDD = ixp.zerorank3()
    KDDdD = AB.KDDdD
    for j in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                rank3term2DDD[j][k][l] = sp.Rational(1, 2) * (KDDdD[j][k][l] - KDDdD[j][l][k])
    # ... then we construct the second term in this sum:
    #  \Gamma^{p}_{j[k} K_{l]p} = \frac{1}{2} (\Gamma^{p}_{jk} K_{lp}-\Gamma^{p}_{jl} K_{kp}):
    for j in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                for p in range(DIM):
                    rank3term2DDD[j][k][l] += sp.Rational(1, 2) * (
                                GammaUDD[p][j][k] * KDD[l][p] - GammaUDD[p][j][l] * KDD[k][p])
    # Finally, we multiply the term by $-8$:
    for j in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                rank3term2DDD[j][k][l] *= sp.sympify(-8)
    # Step 5: Construct the (rank-2) tensor in term 3 of psi_4 (referring to Eq 5.1 in 
    #   Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf
    # Step 5.1: Construct 3-Ricci tensor R_{ij} = gamma^{im} R_{ijml}
    RDD = ixp.zerorank2()
    gammaUU = AB.gammaUU
    for j in range(DIM):
        for l in range(DIM):
            for i in range(DIM):
                for m in range(DIM):
                    RDD[j][l] += gammaUU[i][m]*RDDDD[i][j][m][l]
    # Step 5.2: Construct K^p_l = gamma^{pi} K_{il}
    KUD = ixp.zerorank2()
    for p in range(DIM):
        for l in range(DIM):
            for i in range(DIM):
                KUD[p][l] += gammaUU[p][i]*KDD[i][l]
    # Step 5.3: Construct trK = gamma^{ij} K_{ij}
    trK = sp.sympify(0)
    for i in range(DIM):
        for j in range(DIM):
            trK += gammaUU[i][j]*KDD[i][j]
    # Next we put these terms together to construct the entire term in parentheses:
    # +4 \left(R_{jl} - K_{jp} K^p_l + K K_{jl} \right),
    rank2term3DD = ixp.zerorank2()
    for j in range(DIM):
        for l in range(DIM):
            rank2term3DD[j][l] = RDD[j][l] + trK*KDD[j][l]
            for p in range(DIM):
                rank2term3DD[j][l] += - KDD[j][p]*KUD[p][l]
    # Finally we multiply by +4:
    for j in range(DIM):
        for l in range(DIM):
            rank2term3DD[j][l] *= sp.sympify(4)
    # Step 6: Construct real & imaginary parts of psi_4
    #         by contracting constituent rank 2, 3, and 4
    #         tensors with input tetrads mre4U, mim4U, & n4U.
    def tetrad_product__Real_psi4(n, Mre, Mim, mu, nu, eta, delta):
        return +n[mu] * Mre[nu] * n[eta] * Mre[delta] - n[mu] * Mim[nu] * n[eta] * Mim[delta]
    def tetrad_product__Imag_psi4(n, Mre, Mim, mu, nu, eta, delta):
        return -n[mu] * Mre[nu] * n[eta] * Mim[delta] - n[mu] * Mim[nu] * n[eta] * Mre[delta]
    # We split psi_4 into three pieces, to expedite & possibly parallelize C code generation.
    psi4_re_pt = [sp.sympify(0),sp.sympify(0),sp.sympify(0)]
    psi4_im_pt = [sp.sympify(0),sp.sympify(0),sp.sympify(0)]
    # First term:
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    psi4_re_pt[0] += rank4term1DDDD[i][j][k][l] * tetrad_product__Real_psi4(n4U, mre4U, mim4U,
                                                                                            i + 1, j + 1,k + 1, l + 1)
                    psi4_im_pt[0] += rank4term1DDDD[i][j][k][l] * tetrad_product__Imag_psi4(n4U, mre4U, mim4U,
                                                                                            i + 1, j + 1, k + 1, l + 1)
    # Second term:
    for j in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                psi4_re_pt[1] += rank3term2DDD[j][k][l] * \
                           sp.Rational(1, 2) * (+tetrad_product__Real_psi4(n4U, mre4U, mim4U, 0, j + 1, k + 1, l + 1)
                                                - tetrad_product__Real_psi4(n4U, mre4U, mim4U, j + 1, 0, k + 1, l + 1))
                psi4_im_pt[1] += rank3term2DDD[j][k][l] * \
                           sp.Rational(1, 2) * (+tetrad_product__Imag_psi4(n4U, mre4U, mim4U, 0, j + 1, k + 1, l + 1)
                                                - tetrad_product__Imag_psi4(n4U, mre4U, mim4U, j + 1, 0, k + 1, l + 1))
    # Third term:
    for j in range(DIM):
        for l in range(DIM):
            psi4_re_pt[2] += rank2term3DD[j][l] * \
                       (sp.Rational(1, 4) * (+tetrad_product__Real_psi4(n4U, mre4U, mim4U, 0, j + 1, 0, l + 1)
                                             - tetrad_product__Real_psi4(n4U, mre4U, mim4U, j + 1, 0, 0, l + 1)
                                             - tetrad_product__Real_psi4(n4U, mre4U, mim4U, 0, j + 1, l + 1, 0)
                                             + tetrad_product__Real_psi4(n4U, mre4U, mim4U, j + 1, 0, l + 1, 0)))
            psi4_im_pt[2] += rank2term3DD[j][l] * \
                       (sp.Rational(1, 4) * (+tetrad_product__Imag_psi4(n4U, mre4U, mim4U, 0, j + 1, 0, l + 1)
                                             - tetrad_product__Imag_psi4(n4U, mre4U, mim4U, j + 1, 0, 0, l + 1)
                                             - tetrad_product__Imag_psi4(n4U, mre4U, mim4U, 0, j + 1, l + 1, 0)
                                             + tetrad_product__Imag_psi4(n4U, mre4U, mim4U, j + 1, 0, l + 1, 0)))

# As documented in the NRPy+ tutorial module
#   Tutorial-Psi4_tetrads.ipynb,
#   this module will construct tetrads
#   needed to compute \psi_4 (as well as other
#   Weyl scalars and invariants in principle)
# Authors: Zachariah B. Etienne
#          (zachetie **at** gmail **dot* com),
#          and Patrick Nelson
# Step 1.a: import all needed modules from NRPy+:
import sympy as sp
import NRPy_param_funcs as par
import indexedexp as ixp
import grid as gri
import finite_difference as fin
import reference_metric as rfm
# Step 1.b: Initialize TetradChoice parameter
thismodule = __name__
# Current option: QuasiKinnersley = choice made in Baker, Campanelli, and Lousto. PRD 65, 044001 (2002)
par.initialize_param(par.glb_param("char", thismodule, "TetradChoice", "QuasiKinnersley"))
par.initialize_param(par.glb_param("char", thismodule, "UseCorrectUnitNormal", "False"))
def Psi4_tetrads():
    global l4U, n4U, mre4U, mim4U
    # Step 1.c: Check if tetrad choice is implemented:
    if par.parval_from_str(thismodule+"::TetradChoice") != "QuasiKinnersley":
        print("ERROR: "+thismodule+"::TetradChoice = "+par.parval_from_str("TetradChoice")+" currently unsupported!")
        exit(1)
    # Step 1.d: Given the chosen coordinate system, set up
    #           corresponding reference metric and needed
    #           reference metric quantities
    # The following function call sets up the reference metric
    #    and related quantities, including rescaling matrices ReDD,
    #    ReU, and hatted quantities.
    rfm.reference_metric()
    # Step 1.e: Set spatial dimension (must be 3 for BSSN, as BSSN is
    #           a 3+1-dimensional decomposition of the general
    #           relativistic field equations)
    DIM = 3
    # Step 1.f: Import all ADM quantities as written in terms of BSSN quantities
    # import BSSN.ADM_in_terms_of_BSSN as AB
    # AB.ADM_in_terms_of_BSSN()
    # Step 2.a: Declare the Cartesian x,y,z as input parameters
    #           and v_1^a, v_2^a, and v_3^a tetrads,
    #           as well as detgamma and gammaUU from
    #           BSSN.ADM_in_terms_of_BSSN
    x,y,z = gri.register_gridfunctions("AUX",["x","y","z"])
# x, y, z = par.Cparameters("REAL", thismodule, ["x", "y", "z"])
    v1UCart = ixp.zerorank1()
    v2UCart = ixp.zerorank1()
    gammaDD = ixp.declarerank2("gammaDD","sym01")
    gammaUU,detgamma = ixp.symm_matrix_inverter3x3(gammaDD)
    # detgamma = AB.detgamma
    # gammaUU = AB.gammaUU
    # Step 2.b: Define v1U and v2U
    v1UCart = [-y, x, sp.sympify(0)]
    v2UCart = [x, y, z]
    v1U = ixp.zerorank1()
    v2U = ixp.zerorank1()
    for i in range(DIM):
        v1U[i] = v1UCart[i]
        v2U[i] = v2UCart[i]
    # # Step 2.c: Construct the Jacobian d x_Cart^i / d xx^j
    # Jac_dUCart_dDrfmUD = ixp.zerorank2()
    # for i in range(DIM):
    #     for j in range(DIM):
    #         Jac_dUCart_dDrfmUD[i][j] = sp.diff(rfm.xxCart[i], rfm.xx[j])
    #
    # # Step 2.d: Invert above Jacobian to get needed d xx^j / d x_Cart^i
    # Jac_dUrfm_dDCartUD, dummyDET = ixp.generic_matrix_inverter3x3(Jac_dUCart_dDrfmUD)
    #
    # # Step 2.e: Transform v1U and v2U from the Cartesian to the xx^i basis
    # v1U = ixp.zerorank1()
    # v2U = ixp.zerorank1()
    # for i in range(DIM):
    #     for j in range(DIM):
    #         v1U[i] += Jac_dUrfm_dDCartUD[i][j] * v1UCart[j]
    #         v2U[i] += Jac_dUrfm_dDCartUD[i][j] * v2UCart[j]
    # Step 2.f: Define the rank-3 version of the Levi-Civita symbol. Amongst
    #         other uses, this is needed for the construction of the approximate
    #         quasi-Kinnersley tetrad.
    def define_LeviCivitaSymbol_rank3(DIM=-1):
        if DIM == -1:
            DIM = par.parval_from_str("DIM")
        LeviCivitaSymbol = ixp.zerorank3()
        for i in range(DIM):
            for j in range(DIM):
                for k in range(DIM):
                    # From https://codegolf.stackexchange.com/questions/160359/levi-civita-symbol :
                    LeviCivitaSymbol[i][j][k] = (i - j) * (j - k) * (k - i) / 2
        return LeviCivitaSymbol
    # Step 2.g: Define v3U
    v3U = ixp.zerorank1()
    LeviCivitaSymbolDDD = define_LeviCivitaSymbol_rank3(DIM=3)
    for a in range(DIM):
        for b in range(DIM):
            for c in range(DIM):
                for d in range(DIM):
                    v3U[a] += sp.sqrt(detgamma) * gammaUU[a][d] * LeviCivitaSymbolDDD[d][b][c] * v1U[b] * v2U[c]
    # Step 2.h: Define omega_{ij}
    omegaDD = ixp.zerorank2()
    #gammaDD = AB.gammaDD
    def v_vectorDU(v1U,v2U,v3U,  i,a):
        if i==0:
            return v1U[a]
        elif i==1:
            return v2U[a]
        elif i==2:
            return v3U[a]
        else:
            print("ERROR: unknown vector!")
            exit(1)
    def update_omega(omegaDD, i,j, v1U,v2U,v3U,gammaDD):
        omegaDD[i][j] = sp.sympify(0)
        for a in range(DIM):
            for b in range(DIM):
                omegaDD[i][j] += v_vectorDU(v1U,v2U,v3U, i,a)*v_vectorDU(v1U,v2U,v3U, j,b)*gammaDD[a][b]
    # Step 2.i: Define e^a_i. Note that:
    #           omegaDD[0][0] = \omega_{11} above; 
    #           omegaDD[1][1] = \omega_{22} above, etc.
    # First e_1^a: Orthogonalize & normalize:
    e1U = ixp.zerorank1()
    update_omega(omegaDD, 0,0, v1U,v2U,v3U,gammaDD)
    for a in range(DIM):
        e1U[a] = v1U[a]/sp.sqrt(omegaDD[0][0])
    # Next e_2^a: First orthogonalize:
    e2U = ixp.zerorank1()
    update_omega(omegaDD, 0,1, e1U,v2U,v3U,gammaDD)
    for a in range(DIM):
        e2U[a] = (v2U[a] - omegaDD[0][1]*e1U[a])
    # Then normalize:
    update_omega(omegaDD, 1,1, e1U,e2U,v3U,gammaDD)
    for a in range(DIM):
        e2U[a] /= sp.sqrt(omegaDD[1][1])
    # Next e_3^a: First orthogonalize:
    e3U = ixp.zerorank1()
    update_omega(omegaDD, 0,2, e1U,e2U,v3U,gammaDD)
    update_omega(omegaDD, 1,2, e1U,e2U,v3U,gammaDD)
    for a in range(DIM):
        e3U[a] = (v3U[a] - omegaDD[0][2]*e1U[a] - omegaDD[1][2]*e2U[a])
    # Then normalize:
    update_omega(omegaDD, 2,2, e1U,e2U,e3U,gammaDD)
    for a in range(DIM):
        e3U[a] /= sp.sqrt(omegaDD[2][2])
    # Step 2.j: Construct l^mu, n^mu, and m^mu, based on r^mu, theta^mu, phi^mu, and u^mu:
    r4U     = ixp.zerorank1(DIM=4)
    u4U     = ixp.zerorank1(DIM=4)
    theta4U = ixp.zerorank1(DIM=4)
    phi4U   = ixp.zerorank1(DIM=4)
    for a in range(DIM):
        r4U[    a+1] = e2U[a]
        theta4U[a+1] = e3U[a]
        phi4U[  a+1] = e1U[a]
    # FIXME? assumes alpha=1, beta^i = 0
    if par.parval_from_str(thismodule+"::UseCorrectUnitNormal") == "False":
        u4U[0] = 1
    else:
        # Eq. 2.116 in Baumgarte & Shapiro:
        #  n^mu = {1/alpha, -beta^i/alpha}. Note that n_mu = {alpha,0}, so n^mu n_mu = -1.
        import BSSN.BSSN_quantities as Bq
        Bq.declare_BSSN_gridfunctions_if_not_declared_already()
        Bq.BSSN_basic_tensors()
        u4U[0] = 1/Bq.alpha
        for i in range(DIM):
            u4U[i+1] = -Bq.betaU[i]/Bq.alpha
    l4U = ixp.zerorank1(DIM=4)
    n4U = ixp.zerorank1(DIM=4)
    mre4U  = ixp.zerorank1(DIM=4)
    mim4U  = ixp.zerorank1(DIM=4)
    isqrt2 = 1/sp.sqrt(2)
    for mu in range(4):
        l4U[mu]   = isqrt2*(u4U[mu] + r4U[mu])
        n4U[mu]   = isqrt2*(u4U[mu] - r4U[mu])
        mre4U[mu] = isqrt2*theta4U[mu]
        mim4U[mu] = isqrt2*  phi4U[mu]

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Start-to-Finish validation of $\\psi_4$ in curvilinear coordinates against Cartesian formulation provided by [Patrick Nelson's Weyl scalars & invariants in Cartesian coordinates module](../../Tutorial-WeylScalarsInvariants-Cartesian.ipynb)\n",
    "\n",
    "### Author: Zach Etienne\n",
    "\n",
    "<font color='blue'>**This module exists as a modification of [the NRPy+ $\\psi_4$ in curvilinear coordinates module](../../Tutorial-Psi4.ipynb), writing all spacetime quantities in terms of ADM variables and their derivatives directly.**</font>\n",
    "\n",
    "## A Note on Notation\n",
    "\n",
    "As is standard in NRPy+, \n",
    "\n",
    "* Greek indices range from 0 to 3, inclusive, with the zeroth component denoting the temporal (time) component.\n",
    "* Latin indices range from 0 to 2, inclusive, with the zeroth component denoting the first spatial component.\n",
    "\n",
    "As a corollary, any expressions involving mixed Greek and Latin indices will need to offset one set of indices by one: A Latin index in a four-vector will be incremented and a Greek index in a three-vector will be decremented (however, the latter case does not occur in this tutorial module).\n",
    "\n",
    "<a id='toc'></a>\n",
    "\n",
    "# Introduction, Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This module constructs $\\psi_4$, a quantity that is immensely useful when extracting gravitational wave content from a numerical relativity simulation. $\\psi_4$ is related to the gravitational wave strain via\n",
    "\n",
    "$$\n",
    "\\psi_4 = \\ddot{h}_+ - i \\ddot{h}_\\times.\n",
    "$$\n",
    "\n",
    "We construct $\\psi_4$ from the standard ADM spatial metric $\\gamma_{ij}$ and extrinsic curvature $K_{ij}$, and their derivatives. The full expression is given by Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf):\n",
    "\n",
    "\\begin{align}\n",
    "\\psi_4 &= \\left[ {R}_{ijkl}+2K_{i[k}K_{l]j}\\right]\n",
    "{n}^i\\bar{m}^j{n}^k\\bar{m}^l  \\\\\n",
    "& -8\\left[ K_{j[k,l]}+{\\Gamma }_{j[k}^pK_{l]p}\\right]\n",
    "{n}^{[0}\\bar{m}^{j]}{n}^k\\bar{m}^l \\\\\n",
    "& +4\\left[ {R}_{jl}-K_{jp}K_l^p+KK_{jl}\\right]\n",
    "{n}^{[0}\\bar{m}^{j]}{n}^{[0}\\bar{m}^{l]},\n",
    "\\end{align}\n",
    "\n",
    "Note that $\\psi_4$ is imaginary, with the imaginary components originating from the tetrad vector $m^\\mu$. This module does not specify a tetrad; instead it only constructs the above expression leaving $m^\\mu$ and $n^\\mu$ unspecified. The [next module on tetrads defines these tetrad quantities](Tutorial-Psi4_tetrads.ipynb) (currently only a quasi-Kinnersley tetrad is supported).\n",
    "\n",
    "**This tutorial module is organized as follows:**\n",
    "\n",
    "1. [Step 1](#initializenrpy): Initialize needed NRPy+ modules\n",
    "1. [Step 2](#riemann): Constructing the 3-Riemann tensor $R_{ik\\ell m}$\n",
    "1. [Step 3](#termone): Constructing the rank-4 tensor in Term 1 of $\\psi_4$: $R_{ijkl} + 2 K_{i[k} K_{l]j}$\n",
    "1. [Step 4](#termtwo): Constructing the rank-3 tensor in Term 2 of $\\psi_4$: $-8 \\left(K_{j[k,l]} + \\Gamma^{p}_{j[k} K_{l]p} \\right)$\n",
    "1. [Step 5](#termthree): Constructing the rank-2 tensor in term 3 of $\\psi_4$: $+4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right)$\n",
    "1. [Step 6](#psifour): Constructing $\\psi_4$ through contractions of the above terms with arbitrary tetrad vectors $n^\\mu$ and $m^\\mu$\n",
    "1. [Step 7](#code_validation): Code Validation against BSSN.Psi4 NRPy+ module\n",
    "1. [Step 8](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize core NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "Let's start by importing all the needed modules from NRPy+:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1.a: import all needed modules from NRPy+:\n",
    "import sympy as sp\n",
    "from outputC import *\n",
    "import NRPy_param_funcs as par\n",
    "import indexedexp as ixp\n",
    "import grid as gri\n",
    "import finite_difference as fin\n",
    "import reference_metric as rfm\n",
    "\n",
    "# Step 1.b: Set the coordinate system for the numerical grid\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n",
    "\n",
    "# Step 1.c: Given the chosen coordinate system, set up \n",
    "#           corresponding reference metric and needed\n",
    "#           reference metric quantities\n",
    "# The following function call sets up the reference metric\n",
    "#    and related quantities, including rescaling matrices ReDD,\n",
    "#    ReU, and hatted quantities.\n",
    "rfm.reference_metric()\n",
    "\n",
    "# Step 1.d: Set spatial dimension (must be 3 for BSSN, as BSSN is \n",
    "#           a 3+1-dimensional decomposition of the general \n",
    "#           relativistic field equations)\n",
    "DIM = 3\n",
    "\n",
    "# Step 1.e: Import all ADM quantities as written in terms of BSSN quantities\n",
    "# import BSSN.ADM_in_terms_of_BSSN as AB\n",
    "# AB.ADM_in_terms_of_BSSN()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='riemann'></a>\n",
    "\n",
    "# Step 2: Constructing the 3-Riemann tensor $R_{ik\\ell m}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{riemann}$$\n",
    "\n",
    "Analogously to Christoffel symbols, the Riemann tensor is a measure of the curvature of an $N$-dimensional manifold. Thus the 3-Riemann tensor is not simply a projection of the 4-Riemann tensor (see e.g., Eq. 2.7 of [Campanelli *et al* (1998)](https://arxiv.org/pdf/gr-qc/9803058.pdf) for the relation between 4-Riemann and 3-Riemann), as $N$-dimensional Riemann tensors are meant to define a notion of curvature given only the associated $N$-dimensional metric. \n",
    "\n",
    "So, given the ADM 3-metric, the Riemann tensor in arbitrary dimension is given by the 3-dimensional version of Eq. 1.19 in Baumgarte & Shapiro's *Numerical Relativity*. I.e.,\n",
    "\n",
    "$$\n",
    "R^i_{jkl} = \\partial_k \\Gamma^{i}_{jl} - \\partial_l \\Gamma^{i}_{jk} + \\Gamma^i_{mk} \\Gamma^m_{jl} - \\Gamma^{i}_{ml} \\Gamma^{m}_{jk},\n",
    "$$\n",
    "where $\\Gamma^i_{jk}$ is the Christoffel symbol associated with the 3-metric $\\gamma_{ij}$:\n",
    "\n",
    "$$\n",
    "\\Gamma^l_{ij} = \\frac{1}{2} \\gamma^{lk} \\left(\\gamma_{ki,j} + \\gamma_{kj,i} - \\gamma_{ij,k} \\right) \n",
    "$$\n",
    "\n",
    "Notice that this equation for the Riemann tensor is equivalent to the equation given in the Wikipedia article on [Formulas in Riemannian geometry](https://en.wikipedia.org/w/index.php?title=List_of_formulas_in_Riemannian_geometry&oldid=882667524):\n",
    "\n",
    "$$\n",
    "R^\\ell{}_{ijk}=\n",
    "\\partial_j \\Gamma^\\ell{}_{ik}-\\partial_k\\Gamma^\\ell{}_{ij}\n",
    "+\\Gamma^\\ell{}_{js}\\Gamma_{ik}^s-\\Gamma^\\ell{}_{ks}\\Gamma^s{}_{ij},\n",
    "$$\n",
    "with the replacements $i\\to \\ell$, $j\\to i$, $k\\to j$, $l\\to k$, and $s\\to m$. Wikipedia also provides a simpler form in terms of second-derivatives of three-metric itself (using the definition of Christoffel symbol), so that we need not define derivatives of the Christoffel symbol:\n",
    "\n",
    "$$\n",
    "R_{ik\\ell m}=\\frac{1}{2}\\left(\n",
    "\\gamma_{im,k\\ell} \n",
    "+ \\gamma_{k\\ell,im}\n",
    "- \\gamma_{i\\ell,km}\n",
    "- \\gamma_{km,i\\ell} \\right)\n",
    "+\\gamma_{np} \\left(\n",
    "\\Gamma^n{}_{k\\ell} \\Gamma^p{}_{im} - \n",
    "\\Gamma^n{}_{km} \\Gamma^p{}_{i\\ell} \\right).\n",
    "$$\n",
    "\n",
    "First we construct the term on the left:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2: Construct the (rank-4) Riemann curvature tensor associated with the ADM 3-metric:\n",
    "RDDDD = ixp.zerorank4()\n",
    "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUX\",\"gammaDD\", \"sym01\") # The AUX or EVOL designation is *not*\n",
    "                                                                                # used in diagnostic modules.\n",
    "kDD = ixp.register_gridfunctions_for_single_rank2(\"AUX\",\"kDD\", \"sym01\")\n",
    "gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\")\n",
    "gammaDD_dDD = ixp.declarerank4(\"gammaDD_dDD\",\"sym01_sym23\")\n",
    "\n",
    "# gammaDD_dDD = AB.gammaDD_dDD\n",
    "\n",
    "for i in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                RDDDD[i][k][l][m] = sp.Rational(1,2) * \\\n",
    "                    (gammaDD_dDD[i][m][k][l] + gammaDD_dDD[k][l][i][m] - gammaDD_dDD[i][l][k][m] - gammaDD_dDD[k][m][i][l])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... then we add the term on the right:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ... then we add the term on the right:\n",
    "# Define the Christoffel symbols\n",
    "GammaUDD = ixp.zerorank3(DIM)\n",
    "gammaUU,gammadetdummy = ixp.symm_matrix_inverter3x3(gammaDD)\n",
    "for i in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                GammaUDD[i][k][l] += (sp.Rational(1,2))*gammaUU[i][m]*\\\n",
    "                                     (gammaDD_dD[m][k][l] + gammaDD_dD[m][l][k] - gammaDD_dD[k][l][m])\n",
    "\n",
    "for i in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                for n in range(DIM):\n",
    "                    for p in range(DIM):\n",
    "                        RDDDD[i][k][l][m] += gammaDD[n][p] * \\\n",
    "                            (GammaUDD[n][k][l]*GammaUDD[p][i][m] - GammaUDD[n][k][m]*GammaUDD[p][i][l])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='termone'></a>\n",
    "\n",
    "# Step 3: Constructing the rank-4 tensor in Term 1 of $\\psi_4$: $R_{ijkl} + 2 K_{i[k} K_{l]j}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{termone}$$\n",
    "\n",
    "Following Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf), the rank-4 tensor in the first term of $\\psi_4$ is given by\n",
    "\n",
    "$$\n",
    "R_{ijkl} + 2 K_{i[k} K_{l]j} = R_{ijkl} + K_{ik} K_{lj} - K_{il} K_{kj}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3: Construct the (rank-4) tensor in term 1 of psi_4 (referring to Eq 5.1 in \n",
    "#   Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf\n",
    "rank4term1 = ixp.zerorank4()\n",
    "# kDD = AB.kDD\n",
    "\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            for l in range(DIM):\n",
    "                rank4term1[i][j][k][l] = RDDDD[i][j][k][l] + kDD[i][k]*kDD[l][j] - kDD[i][l]*kDD[k][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='termtwo'></a>\n",
    "\n",
    "# Step 4: Constructing the rank-3 tensor in Term 2 of $\\psi_4$: $-8 \\left(K_{j[k,l]} + \\Gamma^{p}_{j[k} K_{l]p} \\right)$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{termtwo}$$\n",
    "\n",
    "Following Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf), the rank-3 tensor in the second term of $\\psi_4$ is given by\n",
    "\n",
    "$$\n",
    "-8 \\left(K_{j[k,l]} + \\Gamma^{p}_{j[k} K_{l]p} \\right)\n",
    "$$\n",
    "First let's construct the first term in this sum: $K_{j[k,l]} = \\frac{1}{2} (K_{jk,l} - K_{jl,k})$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4: Construct the (rank-3) tensor in term 2 of psi_4 (referring to Eq 5.1 in \n",
    "#   Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf\n",
    "rank3term2 = ixp.zerorank3()\n",
    "# kDD_dD = AB.kDD_dD\n",
    "kDD_dD = ixp.declarerank3(\"kDD_dD\",\"sym01\")\n",
    "\n",
    "for j in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            rank3term2[j][k][l] = sp.Rational(1,2)*(kDD_dD[j][k][l] - kDD_dD[j][l][k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... then we construct the second term in this sum: $\\Gamma^{p}_{j[k} K_{l]p} = \\frac{1}{2} (\\Gamma^{p}_{jk} K_{lp}-\\Gamma^{p}_{jl} K_{kp})$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ... then we construct the second term in this sum: \n",
    "#  \\Gamma^{p}_{j[k} K_{l]p} = \\frac{1}{2} (\\Gamma^{p}_{jk} K_{lp}-\\Gamma^{p}_{jl} K_{kp}):\n",
    "for j in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            for p in range(DIM):\n",
    "                rank3term2[j][k][l] += sp.Rational(1,2)*(GammaUDD[p][j][k]*kDD[l][p] - GammaUDD[p][j][l]*kDD[k][p])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we multiply the term by $-8$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finally, we multiply the term by $-8$:\n",
    "for j in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            rank3term2[j][k][l] *= sp.sympify(-8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='termthree'></a>\n",
    "\n",
    "# Step 5: Constructing the rank-2 tensor in term 3 of $\\psi_4$: $+4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right)$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{termthree}$$\n",
    "\n",
    "Following Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf), the rank-2 tensor in the third term of $\\psi_4$ is given by\n",
    "\n",
    "$$\n",
    "+4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right),\n",
    "$$\n",
    "where\n",
    "\\begin{align}\n",
    "R_{jl} &= R^i_{jil} \\\\\n",
    "&= \\gamma^{im} R_{ijml} \\\\\n",
    "K &= K^i_i \\\\\n",
    "&= \\gamma^{im} K_{im}\n",
    "\\end{align}\n",
    "\n",
    "Let's build the components of this term: $R_{jl}$, $K^p_l$, and $K$, as defined above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 5: Construct the (rank-2) tensor in term 3 of psi_4 (referring to Eq 5.1 in \n",
    "#   Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf\n",
    "\n",
    "# Step 5.1: Construct 3-Ricci tensor R_{ij} = gamma^{im} R_{ijml}\n",
    "RDD = ixp.zerorank2()\n",
    "for j in range(DIM):\n",
    "    for l in range(DIM):\n",
    "        for i in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                RDD[j][l] += gammaUU[i][m]*RDDDD[i][j][m][l]\n",
    "\n",
    "# Step 5.2: Construct K^p_l = gamma^{pi} K_{il}\n",
    "KUD = ixp.zerorank2()\n",
    "for p in range(DIM):\n",
    "    for l in range(DIM):\n",
    "        for i in range(DIM):\n",
    "            KUD[p][l] += gammaUU[p][i]*kDD[i][l]\n",
    "\n",
    "# Step 5.3: Construct trK = gamma^{ij} K_{ij}\n",
    "trK = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        trK += gammaUU[i][j]*kDD[i][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we put these terms together to construct the entire term:\n",
    "$$\n",
    "+4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right),\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Next we put these terms together to construct the entire term in parentheses:\n",
    "# +4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right),\n",
    "rank2term3 = ixp.zerorank2()\n",
    "for j in range(DIM):\n",
    "    for l in range(DIM):\n",
    "        rank2term3[j][l] = RDD[j][l] + trK*kDD[j][l]\n",
    "        for p in range(DIM):\n",
    "            rank2term3[j][l] += - kDD[j][p]*KUD[p][l]\n",
    "# Finally we multiply by +4:\n",
    "for j in range(DIM):\n",
    "    for l in range(DIM):\n",
    "        rank2term3[j][l] *= sp.sympify(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='psifour'></a>\n",
    "\n",
    "# Step 6: Constructing $\\psi_4$ through contractions of the above terms with an arbitrary tetrad vectors $m^\\mu$ and $n^\\mu$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{psifour}$$\n",
    "\n",
    "Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf) writes $\\psi_4$ (which is complex) as the contraction of each of the above terms with products of tetrad vectors:\n",
    "\n",
    "\\begin{align}\n",
    "\\psi_4 &= \\left[ {R}_{ijkl}+2K_{i[k}K_{l]j}\\right]\n",
    "{n}^i\\bar{m}^j{n}^k\\bar{m}^l  \\\\\n",
    "& -8\\left[ K_{j[k,l]}+{\\Gamma }_{j[k}^pK_{l]p}\\right]\n",
    "{n}^{[0}\\bar{m}^{j]}{n}^k\\bar{m}^l \\\\\n",
    "& +4\\left[ {R}_{jl}-K_{jp}K_l^p+KK_{jl}\\right]\n",
    "{n}^{[0}\\bar{m}^{j]}{n}^{[0}\\bar{m}^{l]},\n",
    "\\end{align}\n",
    "where $\\bar{m}^\\mu$ is the complex conjugate of $m^\\mu$, and $n^\\mu$ is real. The third term is given by\n",
    "\\begin{align}\n",
    "{n}^{[0}\\bar{m}^{j]}{n}^{[0}\\bar{m}^{l]}\n",
    "&= \\frac{1}{2}({n}^{0}\\bar{m}^{j} - {n}^{j}\\bar{m}^{0} )\\frac{1}{2}({n}^{0}\\bar{m}^{l} - {n}^{l}\\bar{m}^{0} )\\\\\n",
    "&= \\frac{1}{4}({n}^{0}\\bar{m}^{j} - {n}^{j}\\bar{m}^{0} )({n}^{0}\\bar{m}^{l} - {n}^{l}\\bar{m}^{0} )\\\\\n",
    "&= \\frac{1}{4}({n}^{0}\\bar{m}^{j}{n}^{0}\\bar{m}^{l} - {n}^{j}\\bar{m}^{0}{n}^{0}\\bar{m}^{l} - {n}^{0}\\bar{m}^{j}{n}^{l}\\bar{m}^{0} +  {n}^{j}\\bar{m}^{0}{n}^{l}\\bar{m}^{0})\n",
    "\\end{align}\n",
    "\n",
    "Only $m^\\mu$ is complex, so we can separate the real and imaginary parts of $\\psi_4$ by hand, defining $M^\\mu$ to now be the real part of $\\bar{m}^\\mu$ and $\\mathcal{M}^\\mu$ to be the imaginary part. All of the above products are of the form ${n}^\\mu\\bar{m}^\\nu{n}^\\eta\\bar{m}^\\delta$, so let's evalute the real and imaginary parts of this product once, for all such terms:\n",
    "\n",
    "\\begin{align}\n",
    "{n}^\\mu\\bar{m}^\\nu{n}^\\eta\\bar{m}^\\delta\n",
    "&= {n}^\\mu(M^\\nu - i \\mathcal{M}^\\nu){n}^\\eta(M^\\delta - i \\mathcal{M}^\\delta) \\\\\n",
    "&= \\left({n}^\\mu M^\\nu {n}^\\eta M^\\delta -\n",
    "{n}^\\mu \\mathcal{M}^\\nu {n}^\\eta \\mathcal{M}^\\delta \\right)+\n",
    "i \\left(\n",
    "-{n}^\\mu M^\\nu {n}^\\eta \\mathcal{M}^\\delta\n",
    "-{n}^\\mu \\mathcal{M}^\\nu {n}^\\eta M^\\delta\n",
    "\\right)\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initialize_param() minor warning: Did nothing; already initialized parameter reference_metric::M_PI\n",
      "initialize_param() minor warning: Did nothing; already initialized parameter reference_metric::xmin\n",
      "initialize_param() minor warning: Did nothing; already initialized parameter reference_metric::xmax\n",
      "initialize_param() minor warning: Did nothing; already initialized parameter reference_metric::ymin\n",
      "initialize_param() minor warning: Did nothing; already initialized parameter reference_metric::ymax\n",
      "initialize_param() minor warning: Did nothing; already initialized parameter reference_metric::zmin\n",
      "initialize_param() minor warning: Did nothing; already initialized parameter reference_metric::zmax\n"
     ]
    }
   ],
   "source": [
    "# mre4U = ixp.declarerank1(\"mre4U\",DIM=4)\n",
    "# mim4U = ixp.declarerank1(\"mim4U\",DIM=4)\n",
    "# n4U   = ixp.declarerank1(\"n4U\"  ,DIM=4)\n",
    "\n",
    "import BSSN.Psi4_tetrads as P4t\n",
    "P4t.Psi4_tetrads()\n",
    "mre4U = P4t.mre4U\n",
    "mim4U = P4t.mim4U\n",
    "n4U = P4t.n4U\n",
    "\n",
    "def tetrad_product__Real_psi4(n,Mre,Mim,  mu,nu,eta,delta):\n",
    "    return +n[mu]*Mre[nu]*n[eta]*Mre[delta] - n[mu]*Mim[nu]*n[eta]*Mim[delta]\n",
    "\n",
    "def tetrad_product__Imag_psi4(n,Mre,Mim,  mu,nu,eta,delta):\n",
    "    return -n[mu]*Mre[nu]*n[eta]*Mim[delta] - n[mu]*Mim[nu]*n[eta]*Mre[delta]\n",
    "\n",
    "\n",
    "psi4_re = sp.sympify(0)\n",
    "psi4_im = sp.sympify(0)\n",
    "# First term:\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            for l in range(DIM):\n",
    "                psi4_re += rank4term1[i][j][k][l]*tetrad_product__Real_psi4(n4U,mre4U,mim4U, i+1,j+1,k+1,l+1)\n",
    "                psi4_im += rank4term1[i][j][k][l]*tetrad_product__Imag_psi4(n4U,mre4U,mim4U, i+1,j+1,k+1,l+1)\n",
    "\n",
    "# Second term:\n",
    "for j in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            psi4_re += rank3term2[j][k][l] * \\\n",
    "                       sp.Rational(1,2)*(+tetrad_product__Real_psi4(n4U,mre4U,mim4U, 0,j+1,k+1,l+1)\n",
    "                                         -tetrad_product__Real_psi4(n4U,mre4U,mim4U, j+1,0,k+1,l+1) )\n",
    "            psi4_im += rank3term2[j][k][l] * \\\n",
    "                       sp.Rational(1,2)*(+tetrad_product__Imag_psi4(n4U,mre4U,mim4U, 0,j+1,k+1,l+1)\n",
    "                                         -tetrad_product__Imag_psi4(n4U,mre4U,mim4U, j+1,0,k+1,l+1) )\n",
    "# Third term:\n",
    "for j in range(DIM):\n",
    "    for l in range(DIM):\n",
    "        psi4_re += rank2term3[j][l] * \\\n",
    "                       (sp.Rational(1,4)*(+tetrad_product__Real_psi4(n4U,mre4U,mim4U, 0,j+1,0,l+1)\n",
    "                                          -tetrad_product__Real_psi4(n4U,mre4U,mim4U, j+1,0,0,l+1)\n",
    "                                          -tetrad_product__Real_psi4(n4U,mre4U,mim4U, 0,j+1,l+1,0)\n",
    "                                          +tetrad_product__Real_psi4(n4U,mre4U,mim4U, j+1,0,l+1,0)))\n",
    "        psi4_im += rank2term3[j][l] * \\\n",
    "                       (sp.Rational(1,4)*(+tetrad_product__Imag_psi4(n4U,mre4U,mim4U, 0,j+1,0,l+1)\n",
    "                                          -tetrad_product__Imag_psi4(n4U,mre4U,mim4U, j+1,0,0,l+1)\n",
    "                                          -tetrad_product__Imag_psi4(n4U,mre4U,mim4U, 0,j+1,l+1,0)\n",
    "                                          +tetrad_product__Imag_psi4(n4U,mre4U,mim4U, j+1,0,l+1,0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 6: Code validation against BSSN.Psi4 NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "As a code validation check, we verify agreement in the SymPy expressions for the RHSs of the BSSN equations between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ BSSN.Psi4 module.\n",
    "\n",
    "By default, we compare all quantities in Spherical coordinates, though other coordinate systems may be chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "STARTING NEW\n",
      "Wrote to file \"Psi4_new.h\"\n",
      "FINISHED NEW\n",
      "STARTING OLD\n",
      "Wrote to file \"Psi4_old.h\"\n",
      "FINISHED OLD\n"
     ]
    }
   ],
   "source": [
    "outCparams = \"preindent=1,outCfileaccess=w,outCverbose=False,includebraces=False\"\n",
    "print(\"STARTING NEW\")\n",
    "fin.FD_outputC(\"Psi4_new.h\", lhrh(lhs=\"psi4_real\", rhs=psi4_im), outCparams)\n",
    "print(\"FINISHED NEW\")\n",
    "\n",
    "gri.glb_gridfcs_list = []\n",
    "\n",
    "import WeylScal4NRPy.WeylScalars_Cartesian as W4\n",
    "W4.WeylScalars_Cartesian()\n",
    "print(\"STARTING OLD\")\n",
    "fin.FD_outputC(\"Psi4_old.h\", lhrh(lhs=\"psi4_real\", rhs=W4.psi4i), outCparams)\n",
    "print(\"FINISHED OLD\")\n",
    "\n",
    "# print(\"FullSimplify[\"+str(sp.mathematica_code(psi4_re-W4.psi4r))+\"]\")\n",
    "# with open(\"math.txt\",\"w\") as file:\n",
    "#     file.write(\"FullSimplify[\"+str(sp.mathematica_code(psi4_re-W4.psi4r))+\"]\")\n",
    "    \n",
    "# # Call the BSSN_RHSs() function from within the\n",
    "# #          BSSN/BSSN_RHSs.py module,\n",
    "# #          which should do exactly the same as in Steps 1-16 above.\n",
    "# print(\"vvv Ignore the minor warnings below. vvv\")\n",
    "# import BSSN.Psi4 as BP4\n",
    "# BP4.Psi4()\n",
    "# print(\"^^^ Ignore the minor warnings above. ^^^\\n\")\n",
    "\n",
    "# print(\"Consistency check between this tutorial and BSSN.Psi4 NRPy+ module: ALL SHOULD BE ZERO.\")\n",
    "\n",
    "# print(\"psi4_im - BP4.psi4_im = \" + str(psi4_im - BP4.psi4_im))\n",
    "# print(\"psi4_re - BP4.psi4_re = \" + str(psi4_re - BP4.psi4_re))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GOOD: 9.858968807130199e-16 : 1.678026634904849e-08 1.678026634904847e-08\r\n"
     ]
    }
   ],
   "source": [
    "!gcc -O2 psi4_tester.c -o psi4_tester -lm\n",
    "!./psi4_tester 4 4 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 7: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[NbConvertApp] Converting notebook Tutorial-Psi4.ipynb to latex\n",
      "[NbConvertApp] Writing 61056 bytes to Tutorial-Psi4.tex\n",
      "This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)\n",
      " restricted \\write18 enabled.\n",
      "entering extended mode\n",
      "This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)\n",
      " restricted \\write18 enabled.\n",
      "entering extended mode\n",
      "This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)\n",
      " restricted \\write18 enabled.\n",
      "entering extended mode\n"
     ]
    }
   ],
   "source": [
    "!jupyter nbconvert --to latex --template latex_nrpy_style.tplx --log-level='WARN' Tutorial-Psi4.ipynb\n",
    "!pdflatex -interaction=batchmode Tutorial-Psi4.tex\n",
    "!pdflatex -interaction=batchmode Tutorial-Psi4.tex\n",
    "!pdflatex -interaction=batchmode Tutorial-Psi4.tex\n",
    "!rm -f Tut*.out Tut*.aux Tut*.log"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

# This code calculates the Weyl scalars psi0, psi1, psi2, psi3, and psi4. It does so by following the paper
# Baker, Campanelli, and Lousto. PRD 65, 044001 (2002), gr-qc/0104063 and the example set by the Kranc-
# generated ETK thorn which can be found at https://bitbucket.org/einsteintoolkit/einsteinanalysis/src
# Step 1: import all needed modules from NRPy+:
import NRPy_param_funcs as par
import indexedexp as ixp
import grid as gri
import finite_difference as fin
from outputC import *
import sympy as sp
# Step 2: Initialize WeylScalars parameters
thismodule = __name__
# Current option: Approx_QuasiKinnersley = choice made in Baker, Campanelli, and Lousto. PRD 65, 044001 (2002)
par.initialize_param(par.glb_param("char *", thismodule, "TetradChoice", "Approx_QuasiKinnersley"))
# This controls what gets output. Acceptable values are "psi4_only", "all_psis", and "all_psis_and_invariants"
par.initialize_param(par.glb_param("char *", thismodule, "output_scalars", "all_psis_and_invariants"))
# Step 3: Define the rank-3 version of the Levi-Civita symbol. Amongst
#         other possible uses, this is needed for the construction of the approximate 
#         quasi-Kinnersley tetrad.
def define_LeviCivitaSymbol_rank3(DIM=-1):
    if DIM == -1:
        DIM = par.parval_from_str("DIM")
    LeviCivitaSymbol = ixp.zerorank3()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                # From https://codegolf.stackexchange.com/questions/160359/levi-civita-symbol :
                LeviCivitaSymbol[i][j][k] = (i - j) * (j - k) * (k - i) / 2
    return LeviCivitaSymbol
# Step 1: Call BSSNs. This module computes many different quantities related to the metric,
#         many of which will be essential here. We must first change to our desired coordinate
#         system, however.
def WeylScalars_Cartesian():
    # We do not need the barred or hatted quantities calculated when using Cartesian coordinates.
    # Instead, we declare the PHYSICAL metric and extrinsic curvature as grid functions.
    gammaDD = ixp.register_gridfunctions_for_single_rank2("AUX","gammaDD", "sym01")
    kDD = ixp.register_gridfunctions_for_single_rank2("AUX","kDD", "sym01")
    gammaUU, detgamma = ixp.symm_matrix_inverter3x3(gammaDD)
    
    output_scalars = par.parval_from_str("output_scalars")
    global psi4r,psi4i,psi3r,psi3i,psi2r,psi2i,psi1r,psi1i,psi0r,psi0i
#    if output_scalars is "all_psis_and_invariants":
#        psi4r,psi4i,psi3r,psi3i,psi2r,psi2i,psi1r,psi1i,psi0r,psi0i = sp.symbols("psi4r psi4i\
#                                                                                  psi3r psi3i\
#                                                                                  psi2r psi2i\
#                                                                                  psi1r psi1i\
#                                                                                  psi0r psi0i")
 #   elif output_scalars is "all_psis":
    psi4r,psi4i,psi3r,psi3i,psi2r,psi2i,psi1r,psi1i,psi0r,psi0i = gri.register_gridfunctions("AUX",["psi4r","psi4i",\
                                                                                                    "psi3r","psi3i",\
                                                                                                    "psi2r","psi2i",\
                                                                                                    "psi1r","psi1i",\
                                                                                                    "psi0r","psi0i"])
    # Step 2a: Set spatial dimension (must be 3 for BSSN)
    DIM = 3
    par.set_parval_from_str("grid::DIM",DIM)
    # Step 2b: Set the coordinate system to Cartesian
    x,y,z = gri.register_gridfunctions("AUX",["x","y","z"])
    # Step 2c: Set which tetrad is used; at the moment, only one supported option
    if par.parval_from_str("WeylScal4NRPy.WeylScalars_Cartesian::TetradChoice") == "Approx_QuasiKinnersley":
        # Step 3a: Choose 3 orthogonal vectors. Here, we choose one in the azimuthal 
        #          direction, one in the radial direction, and the cross product of the two. 
        # Eqs 5.6, 5.7 in https://arxiv.org/pdf/gr-qc/0104063.pdf:
        # v_1^a &= [-y,x,0] \\
        # v_2^a &= [x,y,z] \\
        # v_3^a &= {\rm det}(g)^{1/2} g^{ad} \epsilon_{dbc} v_1^b v_2^c,
        v1U = ixp.zerorank1()
        v2U = ixp.zerorank1()
        v3U = ixp.zerorank1()
        v1U[0] = -y
        v1U[1] = x
        v1U[2] = sp.sympify(0)
        v2U[0] = x
        v2U[1] = y
        v2U[2] = z
        LeviCivitaSymbol_rank3 = define_LeviCivitaSymbol_rank3()
        for a in range(DIM):
            for b in range(DIM):
                for c in range(DIM):
                    for d in range(DIM):
                        v3U[a] += sp.sqrt(detgamma) * gammaUU[a][d] * LeviCivitaSymbol_rank3[d][b][c] * v1U[b] *v2U[c]
        # Step 3b: Gram-Schmidt orthonormalization of the vectors.
        # The w_i^a vectors here are used to temporarily hold values on the way to the final vectors e_i^a
        # e_1^a &= \frac{v_1^a}{\omega_{11}} \\
        # e_2^a &= \frac{v_2^a - \omega_{12} e_1^a}{\omega_{22}} \\
        # e_3^a &= \frac{v_3^a - \omega_{13} e_1^a - \omega_{23} e_2^a}{\omega_{33}}, \\
        
        # Normalize the first vector
        w1U = ixp.zerorank1()
        for a in range(DIM):
            w1U[a] = v1U[a]
        omega11 = sp.sympify(0)
        for a in range(DIM):
            for b in range(DIM):
                omega11 += w1U[a] * w1U[b] * gammaDD[a][b]
        e1U = ixp.zerorank1()
        for a in range(DIM):
            e1U[a] = w1U[a] / sp.sqrt(omega11)
        # Subtract off the portion of the first vector along the second, then normalize
        omega12 = sp.sympify(0)
        for a in range(DIM):
            for b in range(DIM):
                omega12 += e1U[a] * v2U[b] * gammaDD[a][b]
        w2U = ixp.zerorank1()
        for a in range(DIM):
            w2U[a] = v2U[a] - omega12*e1U[a]
        omega22 = sp.sympify(0)
        for a in range(DIM):
            for b in range(DIM):
                omega22 += w2U[a] * w2U[b] *gammaDD[a][b]
        e2U = ixp.zerorank1()
        for a in range(DIM):
            e2U[a] = w2U[a] / sp.sqrt(omega22)
        # Subtract off the portion of the first and second vectors along the third, then normalize
        omega13 = sp.sympify(0)
        for a in range(DIM):
            for b in range(DIM):
                omega13 += e1U[a] * v3U[b] * gammaDD[a][b]
        omega23 = sp.sympify(0)
        for a in range(DIM):
            for b in range(DIM):
                omega23 += e2U[a] * v3U[b] * gammaDD[a][b]
        w3U = ixp.zerorank1()
        for a in range(DIM):
            w3U[a] = v3U[a] - omega13*e1U[a] - omega23*e2U[a]
        omega33 = sp.sympify(0)
        for a in range(DIM):
            for b in range(DIM):
                omega33 += w3U[a] * w3U[b] * gammaDD[a][b]
        e3U = ixp.zerorank1()
        for a in range(DIM):
            e3U[a] = w3U[a] / sp.sqrt(omega33)
        # Step 3c: Construct the tetrad itself.
        # Eqs. 5.6:
        # l^a &= \frac{1}{\sqrt{2}} e_2^a \\
        # n^a &= -\frac{1}{\sqrt{2}} e_2^a \\
        # m^a &= \frac{1}{\sqrt{2}} (e_3^a + i e_1^a) \\
        # \overset{*}{m}{}^a &= \frac{1}{\sqrt{2}} (e_3^a - i e_1^a)
        isqrt2 = 1/sp.sqrt(2)
        ltetU = ixp.zerorank1()
        ntetU = ixp.zerorank1()
        #mtetU = ixp.zerorank1()
        #mtetccU = ixp.zerorank1()
        remtetU = ixp.zerorank1() # SymPy does not like trying to take the real/imaginary parts of such a
        immtetU = ixp.zerorank1() # complicated expression as the Weyl scalars, so we will do it ourselves.
        for i in range(DIM):
            ltetU[i] = isqrt2 * e2U[i]
            ntetU[i] = -isqrt2 * e2U[i]
            remtetU[i] = isqrt2 * e3U[i]
            immtetU[i] = isqrt2 * e1U[i]
        nn = isqrt2
    else:
        print("Error: TetradChoice == "+par.parval_from_str("TetradChoice")+" unsupported!")
        exit(1)
    gammaDD_dD = ixp.declarerank3("gammaDD_dD","sym01")
    # Define the Christoffel symbols
    GammaUDD = ixp.zerorank3(DIM)
    for i in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                for m in range(DIM):
                    GammaUDD[i][k][l] += (sp.Rational(1,2))*gammaUU[i][m]*\
                                         (gammaDD_dD[m][k][l] + gammaDD_dD[m][l][k] - gammaDD_dD[k][l][m])
    # Step 4b: Declare and construct the Riemann curvature tensor:
    # R_{abcd} = \frac{1}{2} (\gamma_{ad,cb}+\gamma_{bc,da}-\gamma_{ac,bd}-\gamma_{bd,ac}) 
    #            + \gamma_{je} \Gamma^{j}_{bc}\Gamma^{e}_{ad} - \gamma_{je} \Gamma^{j}_{bd} \Gamma^{e}_{ac}
    gammaDD_dDD = ixp.declarerank4("gammaDD_dDD","sym01_sym23")
    RiemannDDDD = ixp.zerorank4()
    for a in range(DIM):
        for b in range(DIM):
            for c in range(DIM):
                for d in range(DIM):
                    RiemannDDDD[a][b][c][d] = (gammaDD_dDD[a][d][c][b] + \
                                               gammaDD_dDD[b][c][d][a] - \
                                               gammaDD_dDD[a][c][b][d] - \
                                               gammaDD_dDD[b][d][a][c]) / 2
                    for e in range(DIM):
                        for j in range(DIM):
                            RiemannDDDD[a][b][c][d] +=  gammaDD[j][e] * GammaUDD[j][b][c] * GammaUDD[e][a][d] - \
                                                        gammaDD[j][e] * GammaUDD[j][b][d] * GammaUDD[e][a][c]
    # Step 4c: We also need the extrinsic curvature tensor $K_{ij}$. 
    # In Cartesian coordinates, we already made the components gridfunctions.
    # We will, however, need to calculate the trace of K seperately:
    trK = sp.sympify(0)
    for i in range(DIM):
        for j in range(DIM):
            trK += gammaUU[i][j] * kDD[i][j]
    # Step 5: Build the formula for \psi_4.
    # Gauss equation: involving the Riemann tensor and extrinsic curvature.
    # GaussDDDD[i][j][k][l] =& R_{ijkl} + 2K_{i[k}K_{l]j}
    GaussDDDD = ixp.zerorank4()
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    GaussDDDD[i][j][k][l] = RiemannDDDD[i][j][k][l] + kDD[i][k]*kDD[l][j] - kDD[i][l]*kDD[k][j]
    # Codazzi equation: involving partial derivatives of the extrinsic curvature. 
    # We will first need to declare derivatives of kDD
    # CodazziDDD[j][k][l] =& -2 (K_{j[k,l]} + \Gamma^p_{j[k} K_{l]p})
    kDD_dD = ixp.declarerank3("kDD_dD","sym01")
    CodazziDDD = ixp.zerorank3()
    for j in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                CodazziDDD[j][k][l] = kDD_dD[j][l][k] - kDD_dD[j][k][l]
                for p in range(DIM):
                    CodazziDDD[j][k][l] += GammaUDD[p][j][l]*kDD[k][p] - GammaUDD[p][j][k]*kDD[l][p]
    # Another piece. While not associated with any particular equation,
    # this is still useful for organizational purposes.
    # RojoDD[j][l]} = & R_{jl} - K_{jp} K^p_l + KK_{jl} \\
    #               = & \gamma^{pd} R_{jpld} - K_{jp} K^p_l + KK_{jl}
    RojoDD = ixp.zerorank2()
    for j in range(DIM):
        for l in range(DIM):
            RojoDD[j][l] = trK*kDD[j][l]
            for p in range(DIM):
                for d in range(DIM):
                    RojoDD[j][l] += gammaUU[p][d]*RiemannDDDD[j][p][l][d] - kDD[j][p]*gammaUU[p][d]*kDD[d][l]
    # Now we can calculate $\psi_4$ itself! We assume l^0 = n^0 = \frac{1}{\sqrt{2}} 
    # and m^0 = \overset{*}{m}{}^0 = 0 to simplify these equations.
    # We calculate the Weyl scalars as defined in https://arxiv.org/abs/gr-qc/0104063
    # In terms of the above-defined quantites, the psis are defined as:
    # \psi_4 =&\ (\text{GaussDDDD[i][j][k][l]}) n^i \overset{*}{m}{}^j n^k \overset{*}{m}{}^l \\
    #         &+2 (\text{CodazziDDD[j][k][l]}) n^{0} \overset{*}{m}{}^{j} n^k \overset{*}{m}{}^l \\
    #         &+ (\text{RojoDD[j][l]}) n^{0} \overset{*}{m}{}^{j} n^{0} \overset{*}{m}{}^{l}.
    # \psi_3 =&\ (\text{GaussDDDD[i][j][k][l]}) l^i n^j \overset{*}{m}{}^k n^l \\
    #         &+ (\text{CodazziDDD[j][k][l]}) (l^{0} n^{j} \overset{*}{m}{}^k n^l - l^{j} n^{0} \overset{*}{m}{}^k n^l - l^k n^j\overset{*}{m}{}^l n^0) \\
    #         &- (\text{RojoDD[j][l]}) l^{0} n^{j} \overset{*}{m}{}^l n^0 - l^{j} n^{0} \overset{*}{m}{}^l n^0 \\
    # \psi_2 =&\ (\text{GaussDDDD[i][j][k][l]}) l^i m^j \overset{*}{m}{}^k n^l \\
    #         &+ (\text{CodazziDDD[j][k][l]}) (l^{0} m^{j} \overset{*}{m}{}^k n^l - l^{j} m^{0} \overset{*}{m}{}^k n^l - l^k m^l \overset{*}{m}{}^l n^0) \\
    #         &- (\text{RojoDD[j][l]}) l^0 m^j \overset{*}{m}{}^l n^0 \\
    # \psi_1 =&\ (\text{GaussDDDD[i][j][k][l]}) n^i l^j m^k l^l \\
    #         &+ (\text{CodazziDDD[j][k][l]}) (n^{0} l^{j} m^k l^l - n^{j} l^{0} m^k l^l - n^k l^l m^j l^0) \\
    #         &- (\text{RojoDD[j][l]}) (n^{0} l^{j} m^l l^0 - n^{j} l^{0} m^l l^0) \\
    # \psi_0 =&\ (\text{GaussDDDD[i][j][k][l]}) l^i m^j l^k m^l \\
    #         &+2 (\text{CodazziDDD[j][k][l]}) (l^0 m^j l^k m^l + l^k m^l l^0 m^j) \\
    #         &+ (\text{RojoDD[j][l]}) l^0 m^j l^0 m^j. \\
    psi4r = sp.sympify(0)
    psi4i = sp.sympify(0)
    psi3r = sp.sympify(0)
    psi3i = sp.sympify(0)
    psi2r = sp.sympify(0)
    psi2i = sp.sympify(0)
    psi1r = sp.sympify(0)
    psi1i = sp.sympify(0)
    psi0r = sp.sympify(0)
    psi0i = sp.sympify(0)
    for l in range(DIM):
        for j in range(DIM):
            psi4r += RojoDD[j][l] * nn * nn * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])
            psi4i += RojoDD[j][l] * nn * nn * (-remtetU[j]*immtetU[l]-immtetU[j]*remtetU[l])
            psi3r +=-RojoDD[j][l] * nn * nn * (ntetU[j]-ltetU[j]) * remtetU[l]
            psi3i += RojoDD[j][l] * nn * nn * (ntetU[j]-ltetU[j]) * immtetU[l]
            psi2r +=-RojoDD[j][l] * nn * nn * (remtetU[l]*remtetU[j]+immtetU[j]*immtetU[l])
            psi2i +=-RojoDD[j][l] * nn * nn * (immtetU[l]*remtetU[j]-remtetU[j]*immtetU[l])
            psi1r += RojoDD[j][l] * nn * nn * (ntetU[j]*remtetU[l]-ltetU[j]*remtetU[l])
            psi1i += RojoDD[j][l] * nn * nn * (ntetU[j]*immtetU[l]-ltetU[j]*immtetU[l])
            psi0r += RojoDD[j][l] * nn * nn * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])
            psi0i += RojoDD[j][l] * nn * nn * (remtetU[j]*immtetU[l]+immtetU[j]*remtetU[l])
    for l in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                psi4r += 2 * CodazziDDD[j][k][l] * ntetU[k] * nn * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])
                psi4i += 2 * CodazziDDD[j][k][l] * ntetU[k] * nn * (-remtetU[j]*immtetU[l]-immtetU[j]*remtetU[l])
                psi3r += 1 * CodazziDDD[j][k][l] * nn * ((ntetU[j]-ltetU[j])*remtetU[k]*ntetU[l]-remtetU[j]*ltetU[k]*ntetU[l])
                psi3i +=-1 * CodazziDDD[j][k][l] * nn * ((ntetU[j]-ltetU[j])*immtetU[k]*ntetU[l]-immtetU[j]*ltetU[k]*ntetU[l])
                psi2r += 1 * CodazziDDD[j][k][l] * nn * (ntetU[l]*(remtetU[j]*remtetU[k]+immtetU[j]*immtetU[k])-ltetU[k]*(remtetU[j]*remtetU[l]+immtetU[j]*immtetU[l]))
                psi2i += 1 * CodazziDDD[j][k][l] * nn * (ntetU[l]*(immtetU[j]*remtetU[k]-remtetU[j]*immtetU[k])-ltetU[k]*(remtetU[j]*immtetU[l]-immtetU[j]*remtetU[l]))
                psi1r += 1 * CodazziDDD[j][k][l] * nn * (ltetU[j]*remtetU[k]*ltetU[l]-remtetU[j]*ntetU[k]*ltetU[l]-ntetU[j]*remtetU[k]*ltetU[l])
                psi1i += 1 * CodazziDDD[j][k][l] * nn * (ltetU[j]*immtetU[k]*ltetU[l]-immtetU[j]*ntetU[k]*ltetU[l]-ntetU[j]*immtetU[k]*ltetU[l])
                psi0r += 2 * CodazziDDD[j][k][l] * nn * ltetU[k]*(remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])
                psi0i += 2 * CodazziDDD[j][k][l] * nn * ltetU[k]*(remtetU[j]*immtetU[l]+immtetU[j]*remtetU[l])
    for l in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for i in range(DIM):
                    psi4r += GaussDDDD[i][j][k][l] * ntetU[i] * ntetU[k] * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])
                    psi4i += GaussDDDD[i][j][k][l] * ntetU[i] * ntetU[k] * (-remtetU[j]*immtetU[l]-immtetU[j]*remtetU[l])
                    psi3r += GaussDDDD[i][j][k][l] * ltetU[i] * ntetU[j] * remtetU[k] * ntetU[l]
                    psi3i +=-GaussDDDD[i][j][k][l] * ltetU[i] * ntetU[j] * immtetU[k] * ntetU[l]
                    psi2r += GaussDDDD[i][j][k][l] * ltetU[i] * ntetU[l] * (remtetU[j]*remtetU[k]+immtetU[j]*immtetU[k])
                    psi2i += GaussDDDD[i][j][k][l] * ltetU[i] * ntetU[l] * (immtetU[j]*remtetU[k]-remtetU[j]*immtetU[k])
                    psi1r += GaussDDDD[i][j][k][l] * ntetU[i] * ltetU[j] * remtetU[k] * ltetU[l]
                    psi1i += GaussDDDD[i][j][k][l] * ntetU[i] * ltetU[j] * immtetU[k] * ltetU[l]
                    psi0r += GaussDDDD[i][j][k][l] * ltetU[i] * ltetU[k] * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])
                    psi0i += GaussDDDD[i][j][k][l] * ltetU[i] * ltetU[k] * (remtetU[j]*immtetU[l]+immtetU[j]*remtetU[l])