{
 "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",
    "## This notebook presents the construction of $\\psi_4$, a complex scalar for gravitational wave analysis. Using the ADM spatial metric, extrinsic curvature, and arbitrary tetrad vectors, a detailed process is outlined to form $\\psi_4$ following [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf).\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",
    "## Introduction\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)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\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 shutil, os\n",
    "import sys#TylerK: Add sys to get cmdline_helper from NRPy top directory; remove this line and next when debugged\n",
    "sys.path.append('../../')\n",
    "\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 evaluate 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": [],
   "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_re), 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.psi4r), 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": 14,
   "metadata": {},
   "outputs": [],
   "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": [
      "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": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Psi4-Cartesian_validation\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}