Tutorial-WeylScalarsInvariants-Cartesian.ipynb

{
 "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",
    "# Using NRPy+ to Construct SymPy expressions for Weyl scalars and invariants in Cartesian coordinates\n",
    "\n",
    "## Author: Patrick Nelson & Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module adopts the prescription of [Baker, Campanelli, and Lousto. PRD 65, 044001 (2002)](https://arxiv.org/abs/gr-qc/0104063) (henceforth, the \"BCL paper\") to construct the Weyl scalars $\\psi_0$, $\\psi_1$, $\\psi_2$, $\\psi_3$, and $\\psi_4$ from an approximate Kinnersley tetrad. \n",
    "\n",
    "### It also constructs the corresponding Weyl invariants adopting the same approach as Einstein Toolkit's [Kranc-generated](http://kranccode.org/) [WeylScal4 diagnostic module](https://bitbucket.org/einsteintoolkit/einsteinanalysis/src). We will also follow that thorn's approach for other parts of this code.\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** The numerical implementation of expressions constructed in this module have been validated against a trusted code (the WeylScal4 Einstein Toolkit thorn).\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [WeylScal4NRPy/WeylScalarInvariants_Cartesian.py](../edit/WeylScal4NRPy/WeylScalarInvariants_Cartesian.py)\n",
    "* [WeylScal4NRPy/WeylScalars_Cartesian.py](../edit/WeylScal4NRPy/WeylScalars_Cartesian.py)\n",
    "\n",
    "## Introduction: \n",
    "As this module is meant for Cartesian coordinates, all quantities are already rescaled. Further, we assume physical (as opposed to conformal) quantities including the 3-metric $\\gamma_{ij}$ and 3-extrinsic curvature $K_{ij}$ are provided as input gridfunctions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#initializenrpy): Set core NRPy+ parameters for numerical grids and reference metric\n",
    "1. [Step 2](#register_gridfunctions): Defining the Levi-Civita Symbol\n",
    "1. [Step 3](#qktetrad): Defining the Approximate Quasi-Kinnersley Tetrad\n",
    "1. [Step 4](#tensor): Building the Riemann and extrinsic curvature tensors\n",
    "1. [Step 5](#psi4): Putting it all together and calculating $\\psi_4$\n",
    "    1. [Step 5.a](#code_validation1): Code Validation against `WeylScal4NRPy.WeylScalars_Cartesian` NRPy+ Module\n",
    "1. [Step 6](#invariant_scalars): The Invariant Scalars\n",
    "    1. [Step 6.a](#code_validation2): Code Validation against `WeylScal4NRPy.WeylScalarInvariants_Cartesian` NRPy+ Module\n",
    "1. [Step 7](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Set core NRPy+ parameters for numerical grids and reference metric \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:22.099249Z",
     "iopub.status.busy": "2021-03-07T17:20:22.097843Z",
     "iopub.status.idle": "2021-03-07T17:20:22.427220Z",
     "shell.execute_reply": "2021-03-07T17:20:22.426373Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: import all needed modules from NRPy+:\n",
    "import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import grid as gri                # NRPy+: Functions having to do with numerical grids\n",
    "import NRPy_param_funcs as par    # NRPy+: Parameter interface\n",
    "import sys                        # Standard Python modules for multiplatform OS-level functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:22.434793Z",
     "iopub.status.busy": "2021-03-07T17:20:22.433355Z",
     "iopub.status.idle": "2021-03-07T17:20:22.435722Z",
     "shell.execute_reply": "2021-03-07T17:20:22.436212Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2: Initialize WeylScalar parameters\n",
    "thismodule = __name__\n",
    "# Currently only one option: Approx_QuasiKinnersley = choice made in Baker, Campanelli, and Lousto. PRD 65, 044001 (2002)\n",
    "par.initialize_param(par.glb_param(\"char\", thismodule, \"TetradChoice\", \"Approx_QuasiKinnersley\"))\n",
    "# The default value will output all psis\n",
    "par.initialize_param(par.glb_param(\"bool\", thismodule, \"output_all_psis\", False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='register_gridfunctions'></a>\n",
    "\n",
    "# Step 2: Declare input and output variable gridfunctions \\[Back to [top](#toc)\\]\n",
    "$$\\label{register_gridfunctions}$$\n",
    "\n",
    "We declare gridfunctions quantities defined here depend on, including:\n",
    "\n",
    "* the physical metric $\\gamma_{ij}$, \n",
    "* the extrinsic curvature $K_{ij}$, \n",
    "* the Cartesian coordinates $(x,y,z)$,\n",
    "\n",
    "Also, the output gridfunctions are as follows:\n",
    "\n",
    "* the real and imaginary components of $\\psi_4$, and\n",
    "* the Weyl curvature invariants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:22.448012Z",
     "iopub.status.busy": "2021-03-07T17:20:22.447099Z",
     "iopub.status.idle": "2021-03-07T17:20:22.449953Z",
     "shell.execute_reply": "2021-03-07T17:20:22.449412Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.a: Set spatial dimension (must be 3 for BSSN)\n",
    "par.set_parval_from_str(\"grid::DIM\",3)\n",
    "\n",
    "# Step 3.b: declare the additional gridfunctions (i.e., functions whose values are declared\n",
    "#          at every grid point, either inside or outside of our SymPy expressions) needed\n",
    "#         for this thorn:\n",
    "#           * the physical metric $\\gamma_{ij}$,\n",
    "#           * the extrinsic curvature $K_{ij}$,\n",
    "#           * the real and imaginary components of $\\psi_4$, and\n",
    "#           * the Weyl curvature invariants:\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",
    "x,y,z = gri.register_gridfunctions(\"AUX\",[\"x\",\"y\",\"z\"])\n",
    "psi4r,psi4i,psi3r,psi3i,psi2r,psi2i,psi1r,psi1i,psi0r,psi0i = gri.register_gridfunctions(\"AUX\",[\"psi4r\",\"psi4i\",\n",
    "                                                                                                \"psi3r\",\"psi3i\",\n",
    "                                                                                                \"psi2r\",\"psi2i\",\n",
    "                                                                                                \"psi1r\",\"psi1i\",\n",
    "                                                                                                \"psi0r\",\"psi0i\"])\n",
    "curvIr,curvIi,curvJr,curvJi,J1curv,J2curv,J3curv,J4curv = gri.register_gridfunctions(\"AUX\",[\"curvIr\",\"curvIi\",\n",
    "                                                                                            \"curvJr\",\"curvJi\",\n",
    "                                                                                            \"J1curv\",\"J2curv\",\n",
    "                                                                                            \"J3curv\",\"J4curv\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='qktetrad'></a>\n",
    "\n",
    "# Step 3: Defining the Approximate Quasi-Kinnersley Tetrad \\[Back to [top](#toc)\\]\n",
    "$$\\label{qktetrad}$$\n",
    "\n",
    "To define the Weyl scalars, first, a tetrad must be chosen. Below, for compatibility with the [WeylScal4 diagnostic module](https://bitbucket.org/einsteintoolkit/einsteinanalysis/src), we implement the approximate  quasi-Kinnersley tetrad of the [BCL paper](https://arxiv.org/abs/gr-qc/0104063).\n",
    "\n",
    "We begin with the vectors given in eqs. 5.6 and 5.7 of the [BCL paper](https://arxiv.org/abs/gr-qc/0104063),\n",
    "\\begin{align}\n",
    "    v_1^a &= [-y,x,0] \\\\\n",
    "    v_2^a &= [x,y,z] \\\\\n",
    "    v_3^a &= {\\rm det}(g)^{1/2} g^{ad} \\epsilon_{dbc} v_1^b v_2^c,\n",
    "\\end{align}\n",
    "and carry out the Gram-Schmidt orthonormalization process. Note that these vectors are initially orthogonal to each other; one is in the $\\phi$ direction, one is in $r$, and the third is the cross product of the first two. The vectors $w_i^a$ are placeholders in the code; the final product of the orthonormalization is the vectors $e_i^a$. So,\n",
    "\\begin{align}\n",
    "e_1^a &= \\frac{v_1^a}{\\sqrt{\\omega_{11}}} \\\\\n",
    "e_2^a &= \\frac{v_2^a - \\omega_{12} e_1^a}{\\sqrt{\\omega_{22}}} \\\\\n",
    "e_3^a &= \\frac{v_3^a - \\omega_{13} e_1^a - \\omega_{23} e_2^a}{\\sqrt{\\omega_{33}}}, \\\\\n",
    "\\end{align}\n",
    "where $\\omega_{ij} = v_i^a v_j^b \\gamma_{ab}$ needs to be updated between steps (to save resources, we can get away with only calculating components as needed) and uses $e_i^a$ instead of $v_i^a$ if it has been calculated. Recall that $\\gamma_{ab}$ was declared as a gridfunction above.\n",
    "\n",
    "Once we have orthogonal, normalized vectors, we can construct the tetrad itself, again drawing on eqs. 5.6. We could draw on SymPy's built-in tools for complex numbers to build the complex vectors $m^a$ and $(m^*)^a$; however, the final expressions for the Weyl scalars are complex enough that `sp.re()` and `sp.im()` are prohibitively time-consuming. To get around this, we will define the real and imaginary components of $\\overset{*}{m}{}^a$, and do the complex algebra by hand. Thus,\n",
    "\\begin{align}\n",
    "    l^a &= \\frac{1}{\\sqrt{2}} e_2^a \\\\\n",
    "    n^a &= -\\frac{1}{\\sqrt{2}} e_2^a \\\\\n",
    "    m^a &= \\frac{1}{\\sqrt{2}} (e_3^a + i e_1^a) \\\\\n",
    "    \\overset{*}{m}{}^a &= \\frac{1}{\\sqrt{2}} (e_3^a - i e_1^a).\n",
    "\\end{align}\n",
    "\n",
    "In coding this procedure, we will follow the code from [WeylScal4](https://bitbucket.org/einsteintoolkit/einsteinanalysis/src) very closely. We will also assume that $l^0 = n^0 = \\frac{1}{\\sqrt{2}}$ and that $m^0 = \\overset{*}{m}{}^0 = 0$ (again, following the example of the Kranc-generated WeylScal4). This last assumption in particular will significantly reduce the terms needed to find $\\psi_4$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:22.476068Z",
     "iopub.status.busy": "2021-03-07T17:20:22.475269Z",
     "iopub.status.idle": "2021-03-07T17:20:23.805413Z",
     "shell.execute_reply": "2021-03-07T17:20:23.804804Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 4: Set which tetrad is used; at the moment, only one supported option\n",
    "\n",
    "# The tetrad depends in general on the inverse 3-metric gammaUU[i][j]=\\gamma^{ij}\n",
    "#          and the determinant of the 3-metric (detgamma), which are defined in\n",
    "#          the following line of code from gammaDD[i][j]=\\gamma_{ij}.\n",
    "tmpgammaUU, detgamma = ixp.symm_matrix_inverter3x3(gammaDD)\n",
    "detgamma = sp.simplify(detgamma)\n",
    "gammaUU = ixp.zerorank2()\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        gammaUU[i][j] = sp.simplify(tmpgammaUU[i][j])\n",
    "\n",
    "if par.parval_from_str(\"TetradChoice\") == \"Approx_QuasiKinnersley\":\n",
    "    # Eqs 5.6 in https://arxiv.org/pdf/gr-qc/0104063.pdf\n",
    "    xmoved = x# - xorig\n",
    "    ymoved = y# - yorig\n",
    "    zmoved = z# - zorig\n",
    "\n",
    "    # Step 5.a: Choose 3 orthogonal vectors. Here, we choose one in the azimuthal\n",
    "    #          direction, one in the radial direction, and the cross product of the two.\n",
    "    # Eqs 5.7\n",
    "    v1U = ixp.zerorank1()\n",
    "    v2U = ixp.zerorank1()\n",
    "    v3U = ixp.zerorank1()\n",
    "    v1U[0] = -ymoved\n",
    "    v1U[1] = xmoved# + offset\n",
    "    v1U[2] = sp.sympify(0)\n",
    "    v2U[0] = xmoved# + offset\n",
    "    v2U[1] = ymoved\n",
    "    v2U[2] = zmoved\n",
    "    LeviCivitaSymbol_rank3 = ixp.LeviCivitaSymbol_dim3_rank3()\n",
    "    for a in range(3):\n",
    "        for b in range(3):\n",
    "            for c in range(3):\n",
    "                for d in range(3):\n",
    "                    v3U[a] += sp.sqrt(detgamma) * gammaUU[a][d] * LeviCivitaSymbol_rank3[d][b][c] * v1U[b] *v2U[c]\n",
    "\n",
    "    for a in range(3):\n",
    "        v3U[a] = sp.simplify(v3U[a])\n",
    "\n",
    "    # Step 5.b: Gram-Schmidt orthonormalization of the vectors.\n",
    "    # The w_i^a vectors here are used to temporarily hold values on the way to the final vectors e_i^a\n",
    "\n",
    "    # e_1^a &= \\frac{v_1^a}{\\omega_{11}}\n",
    "    # e_2^a &= \\frac{v_2^a - \\omega_{12} e_1^a}{\\omega_{22}}\n",
    "    # e_3^a &= \\frac{v_3^a - \\omega_{13} e_1^a - \\omega_{23} e_2^a}{\\omega_{33}},\n",
    "\n",
    "    # Normalize the first vector\n",
    "    w1U = ixp.zerorank1()\n",
    "    for a in range(3):\n",
    "        w1U[a] = v1U[a]\n",
    "    omega11 = sp.sympify(0)\n",
    "    for a in range(3):\n",
    "        for b in range(3):\n",
    "            omega11 += w1U[a] * w1U[b] * gammaDD[a][b]\n",
    "    e1U = ixp.zerorank1()\n",
    "    for a in range(3):\n",
    "        e1U[a] = w1U[a] / sp.sqrt(omega11)\n",
    "\n",
    "    # Subtract off the portion of the first vector along the second, then normalize\n",
    "    omega12 = sp.sympify(0)\n",
    "    for a in range(3):\n",
    "        for b in range(3):\n",
    "            omega12 += e1U[a] * v2U[b] * gammaDD[a][b]\n",
    "    w2U = ixp.zerorank1()\n",
    "    for a in range(3):\n",
    "        w2U[a] = v2U[a] - omega12*e1U[a]\n",
    "    omega22 = sp.sympify(0)\n",
    "    for a in range(3):\n",
    "        for b in range(3):\n",
    "            omega22 += w2U[a] * w2U[b] *gammaDD[a][b]\n",
    "    e2U = ixp.zerorank1()\n",
    "    for a in range(3):\n",
    "        e2U[a] = w2U[a] / sp.sqrt(omega22)\n",
    "\n",
    "    # Subtract off the portion of the first and second vectors along the third, then normalize\n",
    "    omega13 = sp.sympify(0)\n",
    "    for a in range(3):\n",
    "        for b in range(3):\n",
    "            omega13 += e1U[a] * v3U[b] * gammaDD[a][b]\n",
    "    omega23 = sp.sympify(0)\n",
    "    for a in range(3):\n",
    "        for b in range(3):\n",
    "            omega23 += e2U[a] * v3U[b] * gammaDD[a][b]\n",
    "    w3U = ixp.zerorank1()\n",
    "    for a in range(3):\n",
    "        w3U[a] = v3U[a] - omega13*e1U[a] - omega23*e2U[a]\n",
    "    omega33 = sp.sympify(0)\n",
    "    for a in range(3):\n",
    "        for b in range(3):\n",
    "            omega33 += w3U[a] * w3U[b] * gammaDD[a][b]\n",
    "    e3U = ixp.zerorank1()\n",
    "    for a in range(3):\n",
    "        e3U[a] = w3U[a] / sp.sqrt(omega33)\n",
    "\n",
    "    # Step 5.c: Construct the tetrad itself.\n",
    "    # Eqs. 5.6:\n",
    "    # l^a &= \\frac{1}{\\sqrt{2}} e_2^a\n",
    "    # n^a &= -\\frac{1}{\\sqrt{2}} e_2^a\n",
    "    # m^a &= \\frac{1}{\\sqrt{2}} (e_3^a + i e_1^a)\n",
    "    # \\overset{*}{m}{}^a &= \\frac{1}{\\sqrt{2}} (e_3^a - i e_1^a)\n",
    "    isqrt2 = 1/sp.sqrt(2)\n",
    "    ltetU = ixp.zerorank1()\n",
    "    ntetU = ixp.zerorank1()\n",
    "    #mtetU = ixp.zerorank1()\n",
    "    #mtetccU = ixp.zerorank1()\n",
    "    remtetU = ixp.zerorank1() # SymPy did not like trying to take the real/imaginary parts of such a\n",
    "    immtetU = ixp.zerorank1() # complicated expression, so we do it ourselves.\n",
    "    for i in range(3):\n",
    "        ltetU[i] = isqrt2 * e2U[i]\n",
    "        ntetU[i] = -isqrt2 * e2U[i]\n",
    "        remtetU[i] = isqrt2 * e3U[i]\n",
    "        immtetU[i] = isqrt2 * e1U[i]\n",
    "    nn = isqrt2\n",
    "\n",
    "else:\n",
    "    print(\"Error: TetradChoice == \"+par.parval_from_str(\"TetradChoice\")+\" unsupported!\")\n",
    "    sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='tensor'></a>\n",
    "\n",
    "# Step 4: Building the Riemann and extrinsic curvature tensors \\[Back to [top](#toc)\\]\n",
    "$$\\label{tensor}$$\n",
    "\n",
    "Now that we have the tetrad in place, we can contract it with the Weyl tensor to obtain the Weyl scalars. Naturally, we must first construct the Weyl tensor to do so; we will not do this directly, instead following the example of the [BCL paper](https://arxiv.org/abs/gr-qc/0104063) and the original WeylScal4. We will first build the Christoffel symbols,\n",
    "$$\\Gamma^i_{kl} = \\frac{1}{2} \\gamma^{im} (\\gamma_{mk,l} + \\gamma_{ml,k} - \\gamma_{kl,m}).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:23.830244Z",
     "iopub.status.busy": "2021-03-07T17:20:23.823927Z",
     "iopub.status.idle": "2021-03-07T17:20:23.832494Z",
     "shell.execute_reply": "2021-03-07T17:20:23.833060Z"
    }
   },
   "outputs": [],
   "source": [
    "#Step 5: Declare and construct the second derivative of the metric.\n",
    "gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\")\n",
    "\n",
    "# Define the Christoffel symbols\n",
    "GammaUDD = ixp.zerorank3(3)\n",
    "for i in range(3):\n",
    "    for k in range(3):\n",
    "        for l in range(3):\n",
    "            for m in range(3):\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])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also need the Riemann curvature tensor,\n",
    "$$\n",
    "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},\n",
    "$$\n",
    "since several terms in our expression for $\\psi_4$ are contractions of this tensor.\n",
    "To do this, we need second derivatives of the metric tensor, $\\gamma_{ab,cd}$, using the finite differencing functionality in NRPy+. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:23.887458Z",
     "iopub.status.busy": "2021-03-07T17:20:23.872039Z",
     "iopub.status.idle": "2021-03-07T17:20:24.073173Z",
     "shell.execute_reply": "2021-03-07T17:20:24.072581Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 6.a: Declare and construct the Riemann curvature tensor:\n",
    "gammaDD_dDD = ixp.declarerank4(\"gammaDD_dDD\",\"sym01_sym23\")\n",
    "RiemannDDDD = ixp.zerorank4()\n",
    "for a in range(3):\n",
    "    for b in range(3):\n",
    "        for c in range(3):\n",
    "            for d in range(3):\n",
    "                RiemannDDDD[a][b][c][d] += (gammaDD_dDD[a][d][c][b] +\n",
    "                                            gammaDD_dDD[b][c][d][a] -\n",
    "                                            gammaDD_dDD[a][c][b][d] -\n",
    "                                            gammaDD_dDD[b][d][a][c]) * sp.Rational(1,2)\n",
    "for a in range(3):\n",
    "    for b in range(3):\n",
    "        for c in range(3):\n",
    "            for d in range(3):\n",
    "                for e in range(3):\n",
    "                    for j in range(3):\n",
    "                        RiemannDDDD[a][b][c][d] +=  gammaDD[j][e] * GammaUDD[j][b][c] * GammaUDD[e][a][d] - \\\n",
    "                                                    gammaDD[j][e] * GammaUDD[j][b][d] * GammaUDD[e][a][c]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will need the trace of the extrinsic curvature tensor $K_{ij}$, which can be computed as usual: $K = \\gamma^{ij} K_{ij}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:24.082113Z",
     "iopub.status.busy": "2021-03-07T17:20:24.081453Z",
     "iopub.status.idle": "2021-03-07T17:20:24.083680Z",
     "shell.execute_reply": "2021-03-07T17:20:24.084232Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 6.b: We also need the extrinsic curvature tensor $K_{ij}$. This can be built from quantities from BSSN_RHSs\n",
    "#           For now, we assume this is a gridfunction (We assume the ADM formalism for now).\n",
    "#extrinsicKDD = ixp.zerorank2()\n",
    "#for i in range(3):\n",
    "#    for j in range(3):\n",
    "#        extrinsicKDD[i][j] = (bssn.AbarDD[i][j] + sp.Rational(1,3)*gammaDD[i][j]*bssn.trK)/bssn.exp_m4phi\n",
    "# We will, however, need to calculate the trace of K seperately:\n",
    "trK = sp.sympify(0)\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        trK += gammaUU[i][j] * kDD[i][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='psi4'></a>\n",
    "\n",
    "# Step 5: Putting it all together and calculating $\\psi_4$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{psi4}$$\n",
    "\n",
    "We do not need to explicitly build the Weyl tensor itself, because the [BCL paper](https://arxiv.org/abs/gr-qc/0104063) shows that, for the Weyl tensor $C_{ijkl}$,\n",
    "\\begin{align}\n",
    "\\psi_4 =&\\  C_{ijkl} n^i \\overset{*}{m}{}^j n^k \\overset{*}{m}{}^l\\\\\n",
    "=&\\ (R_{ijkl} + 2K_{i[k}K_{l]j}) n^i \\overset{*}{m}{}^j n^k \\overset{*}{m}{}^l \\\\\n",
    "&- 8 (K_{j[k,l]} + \\Gamma^p_{j[k} K_{l]p}) n^{[0} \\overset{*}{m}{}^{j]} n^k \\overset{*}{m}{}^l \\\\\n",
    "&+ 4 (R_{jl} - K_{jp} K^p_l + KK_{jl}) n^{[0} \\overset{*}{m}{}^{j]} n^{[0} \\overset{*}{m}{}^{l]}.\n",
    "\\end{align}\n",
    "\n",
    "Note here the brackets around pairs of indices. This indicates the antisymmetric part of a tensor; that is, for some arbitrary tensor $A_{ij}$, $A_{[ij]} = \\frac{1}{2}(A_{ij}-A_{ji})$. This applies identically for indices belonging to separate tensors as well as superscripts in place of subscripts.\n",
    "\n",
    "The other Weyl scalars from the [BCL paper](https://arxiv.org/abs/gr-qc/0104063), appendix A, that we may want to consider are \n",
    "\\begin{align}\n",
    "\\psi_3 =&\\ (R_{ijkl} + 2K_{i[k}K_{l]j}) l^i n^j \\overset{*}{m}{}^k n^l \\\\\n",
    "&- 4 (K_{j[k,l]} + \\Gamma^p_{j[k} K_{l]p}) (l^{[0} n^{j]} \\overset{*}{m}{}^k n^l + l^k n^j\\overset{*}{m}{}^{[0} n^{l]}) \\\\\n",
    "&+ 4 (R_{jl} - K_{jp} K^p_l + KK_{jl}) l^{[0} n^{j]} \\overset{*}{m}{}^{[0} n^{l]} \\\\\n",
    "\\psi_2 =&\\ (R_{ijkl} + 2K_{i[k}K_{l]j}) l^i m^j \\overset{*}{m}{}^k n^l \\\\\n",
    "&- 4 (K_{j[k,l]} + \\Gamma^p_{j[k} K_{l]p}) (l^{[0} m^{j]} \\overset{*}{m}{}^k n^l + l^k m^l \\overset{*}{m}{}^{[0} n^{j]}) \\\\\n",
    "&+ 4 (R_{jl} - K_{jp} K^p_l + KK_{jl}) l^{[0} m^{j]} \\overset{*}{m}{}^{[0} n^{l]} \\\\\n",
    "\\psi_1 =&\\ (R_{ijkl} + 2K_{i[k}K_{l]j}) n^i l^j m^k l^l \\\\\n",
    "&- 4 (K_{j[k,l]} + \\Gamma^p_{j[k} K_{l]p}) (n^{[0} l^{j]} m^k l^l + n^k l^l m^{[0} l^{j]}) \\\\\n",
    "&+ 4 (R_{jl} - K_{jp} K^p_l + KK_{jl}) n^{[0} l^{j]} m^{[0} l^{l]} \\\\\n",
    "\\psi_0 =&\\ (R_{ijkl} + 2K_{i[k}K_{l]j}) l^i m^j l^k m^l \\\\\n",
    "&- 8 (K_{j[k,l]} + \\Gamma^p_{j[k} K_{l]p}) l^{[0} m^{j]} l^k m^l \\\\\n",
    "&+ 4 (R_{jl} - K_{jp} K^p_l + KK_{jl}) l^{[0} m^{j]} l^{[0} m^{l]}. \\\\\n",
    "\\end{align}\n",
    "\n",
    "To make it easier to track the construction of this expression, we will break it down into three parts, by first defining each of the parenthetical terms above separately. This is effectively identical to the procedure used in the Mathematica notebook that generates the original [WeylScal4](https://bitbucket.org/einsteintoolkit/einsteinanalysis/src). That is, let \n",
    "\\begin{align}\n",
    "\\text{GaussDDDD[i][j][k][l]} =& R_{ijkl} + 2K_{i[k}K_{l]j},\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:24.111350Z",
     "iopub.status.busy": "2021-03-07T17:20:24.110693Z",
     "iopub.status.idle": "2021-03-07T17:20:24.113187Z",
     "shell.execute_reply": "2021-03-07T17:20:24.113753Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 7: Build the formula for \\psi_4.\n",
    "# Gauss equation: involving the Riemann tensor and extrinsic curvature.\n",
    "GaussDDDD = ixp.zerorank4()\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        for k in range(3):\n",
    "            for l in range(3):\n",
    "                GaussDDDD[i][j][k][l] += RiemannDDDD[i][j][k][l] + kDD[i][k]*kDD[l][j] - kDD[i][l]*kDD[k][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\\text{CodazziDDD[j][k][l]} =& -2 (K_{j[k,l]} + \\Gamma^p_{j[k} K_{l]p}),\n",
    "\\end{align}                                                       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:24.140129Z",
     "iopub.status.busy": "2021-03-07T17:20:24.139446Z",
     "iopub.status.idle": "2021-03-07T17:20:24.141671Z",
     "shell.execute_reply": "2021-03-07T17:20:24.142166Z"
    }
   },
   "outputs": [],
   "source": [
    "# Codazzi equation: involving partial derivatives of the extrinsic curvature.\n",
    "# We will first need to declare derivatives of kDD\n",
    "kDD_dD = ixp.declarerank3(\"kDD_dD\",\"sym01\")\n",
    "CodazziDDD = ixp.zerorank3()\n",
    "for j in range(3):\n",
    "    for k in range(3):\n",
    "        for l in range(3):\n",
    "            CodazziDDD[j][k][l] += kDD_dD[j][l][k] - kDD_dD[j][k][l]\n",
    "\n",
    "for j in range(3):\n",
    "    for k in range(3):\n",
    "        for l in range(3):\n",
    "            for p in range(3):\n",
    "                CodazziDDD[j][k][l] += GammaUDD[p][j][l]*kDD[k][p] - GammaUDD[p][j][k]*kDD[l][p]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and\n",
    "\\begin{align}\n",
    "\\text{RojoDD[j][l]} = & R_{jl} - K_{jp} K^p_l + KK_{jl} \\\\\n",
    "= & \\gamma^{pd} R_{jpld} - K_{jp} K^p_l + KK_{jl},\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:24.189256Z",
     "iopub.status.busy": "2021-03-07T17:20:24.181118Z",
     "iopub.status.idle": "2021-03-07T17:20:24.191334Z",
     "shell.execute_reply": "2021-03-07T17:20:24.191836Z"
    }
   },
   "outputs": [],
   "source": [
    "# Another piece. While not associated with any particular equation,\n",
    "# this is still useful for organizational purposes.\n",
    "RojoDD = ixp.zerorank2()\n",
    "for j in range(3):\n",
    "    for l in range(3):\n",
    "        RojoDD[j][l] += trK*kDD[j][l]\n",
    "\n",
    "for j in range(3):\n",
    "    for l in range(3):\n",
    "        for p in range(3):\n",
    "            for d in range(3):\n",
    "                RojoDD[j][l] += gammaUU[p][d]*RiemannDDDD[j][p][l][d] - kDD[j][p]*gammaUU[p][d]*kDD[d][l]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where these quantities are so named because of their relation to the Gauss-Codazzi equations. Then, we simply contract these with the tetrad we chose earlier to arrive at an expression for $\\psi_4$. The barred Christoffel symbols and barred Ricci tensor have already been calculated by [BSSN/BSSN_RHSs.py](../edit/BSSN/BSSN_RHSs.py), so we use those values. So, our expression for $\\psi_4$ has become \n",
    "\\begin{align}\n",
    "\\psi_4 =&\\ (\\text{GaussDDDD[i][j][k][l]}) n^i \\overset{*}{m}{}^j n^k \\overset{*}{m}{}^l \\\\\n",
    "&+2 (\\text{CodazziDDD[j][k][l]}) n^{0} \\overset{*}{m}{}^{j} n^k \\overset{*}{m}{}^l \\\\\n",
    "&+ (\\text{RojoDD[j][l]}) n^{0} \\overset{*}{m}{}^{j} n^{0} \\overset{*}{m}{}^{l}.\n",
    "\\end{align}\n",
    "\n",
    "Likewise, we can rewrite the other Weyl scalars:\n",
    "\\begin{align}\n",
    "\\psi_3 =&\\ (\\text{GaussDDDD[i][j][k][l]}) l^i n^j \\overset{*}{m}{}^k n^l \\\\\n",
    "&+ (\\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) \\\\\n",
    "&- (\\text{RojoDD[j][l]}) l^{0} n^{j} \\overset{*}{m}{}^l n^0 - l^{j} n^{0} \\overset{*}{m}{}^l n^0 \\\\\n",
    "\\psi_2 =&\\ (\\text{GaussDDDD[i][j][k][l]}) l^i m^j \\overset{*}{m}{}^k n^l \\\\\n",
    "&+ (\\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) \\\\\n",
    "&- (\\text{RojoDD[j][l]}) l^0 m^j \\overset{*}{m}{}^l n^0 \\\\\n",
    "\\psi_1 =&\\ (\\text{GaussDDDD[i][j][k][l]}) n^i l^j m^k l^l \\\\\n",
    "&+ (\\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) \\\\\n",
    "&- (\\text{RojoDD[j][l]}) (n^{0} l^{j} m^l l^0 - n^{j} l^{0} m^l l^0) \\\\\n",
    "\\psi_0 =&\\ (\\text{GaussDDDD[i][j][k][l]}) l^i m^j l^k m^l \\\\\n",
    "&+2 (\\text{CodazziDDD[j][k][l]}) (l^0 m^j l^k m^l + l^k m^l l^0 m^j) \\\\\n",
    "&+ (\\text{RojoDD[j][l]}) l^0 m^j l^0 m^j. \\\\\n",
    "\\end{align}\n",
    "\n",
    "We will start by setting the scalars to SymPy's $0$ (this is done so that Python knows that the scalars are symbolic, not numeric, avoiding some potential bugs later on) and then performing the needed contractions of `RojoDD[j][l]`. Recall that the tetrad vectors were defined above and that we just built `RojoDD[j][l]` from the Ricci tensor and extrinsic curvature.\n",
    "\n",
    "The relevant terms here are:\n",
    "\\begin{align}\n",
    "\\psi_4:&\\ \\ \\ (\\text{RojoDD[j][l]}) n^{0} \\overset{*}{m}{}^{j} n^{0} \\overset{*}{m}{}^{l} \\\\\n",
    "\\psi_3:&\\ -(\\text{RojoDD[j][l]}) (l^{0} n^{j} \\overset{*}{m}{}^l n^0 - l^{j} n^{0} \\overset{*}{m}{}^l n^0) \\\\\n",
    "\\psi_2:&\\ - (\\text{RojoDD[j][l]}) l^0 m^j \\overset{*}{m}{}^l n^0 \\\\\n",
    "\\psi_1:&\\ - (\\text{RojoDD[j][l]}) (n^{0} l^{j} m^l l^0 - n^{j} l^{0} m^l l^0) \\\\\n",
    "\\psi_0:&\\ \\ \\ (\\text{RojoDD[j][l]}) l^0 m^j l^0 m^j. \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:24.266572Z",
     "iopub.status.busy": "2021-03-07T17:20:24.230544Z",
     "iopub.status.idle": "2021-03-07T17:20:24.501201Z",
     "shell.execute_reply": "2021-03-07T17:20:24.501716Z"
    }
   },
   "outputs": [],
   "source": [
    "# Now we can calculate $\\psi_4$ itself!\n",
    "psi4r = sp.sympify(0)\n",
    "psi4i = sp.sympify(0)\n",
    "psi3r = sp.sympify(0)\n",
    "psi3i = sp.sympify(0)\n",
    "psi2r = sp.sympify(0)\n",
    "psi2i = sp.sympify(0)\n",
    "psi1r = sp.sympify(0)\n",
    "psi1i = sp.sympify(0)\n",
    "psi0r = sp.sympify(0)\n",
    "psi0i = sp.sympify(0)\n",
    "for l in range(3):\n",
    "    for j in range(3):\n",
    "        psi4r += RojoDD[j][l] * nn * nn * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])\n",
    "        psi4i += RojoDD[j][l] * nn * nn * (-remtetU[j]*immtetU[l]-immtetU[j]*remtetU[l])\n",
    "        psi3r +=-RojoDD[j][l] * nn * nn * (ntetU[j]-ltetU[j]) * remtetU[l]\n",
    "        psi3i += RojoDD[j][l] * nn * nn * (ntetU[j]-ltetU[j]) * immtetU[l]\n",
    "        psi2r +=-RojoDD[j][l] * nn * nn * (remtetU[l]*remtetU[j]+immtetU[j]*immtetU[l])\n",
    "        psi2i +=-RojoDD[j][l] * nn * nn * (immtetU[l]*remtetU[j]-remtetU[j]*immtetU[l])\n",
    "        psi1r += RojoDD[j][l] * nn * nn * (ntetU[j]*remtetU[l]-ltetU[j]*remtetU[l])\n",
    "        psi1i += RojoDD[j][l] * nn * nn * (ntetU[j]*immtetU[l]-ltetU[j]*immtetU[l])\n",
    "        psi0r += RojoDD[j][l] * nn * nn * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])\n",
    "        psi0i += RojoDD[j][l] * nn * nn * (remtetU[j]*immtetU[l]+immtetU[j]*remtetU[l])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will add the contractions of `CodazziDDD[j][k][l]` to the Weyl Scalars. Again, we use the null tetrad we constructed and the tensor `CodazziDDD[j][k][l]` we constructed from the extrinsic curvature and Christoffel symbols.\n",
    "\n",
    "The relevant terms here are:\n",
    "\\begin{align}\n",
    "\\psi_4:&\\ 2 (\\text{CodazziDDD[j][k][l]}) n^{0} \\overset{*}{m}{}^{j} n^k \\overset{*}{m}{}^l \\\\\n",
    "\\psi_3:&\\ \\ \\ (\\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) \\\\\n",
    "\\psi_2:&\\ \\ \\ (\\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) \\\\\n",
    "\\psi_1:&\\ \\ \\ (\\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) \\\\\n",
    "\\psi_0:&\\ 2 (\\text{CodazziDDD[j][k][l]}) (l^0 m^j l^k m^l + l^k m^l l^0 m^j). \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:24.520355Z",
     "iopub.status.busy": "2021-03-07T17:20:24.515420Z",
     "iopub.status.idle": "2021-03-07T17:20:24.908912Z",
     "shell.execute_reply": "2021-03-07T17:20:24.909416Z"
    }
   },
   "outputs": [],
   "source": [
    "for l in range(3):\n",
    "    for j in range(3):\n",
    "        for k in range(3):\n",
    "            psi4r += 2 * CodazziDDD[j][k][l] * ntetU[k] * nn * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])\n",
    "            psi4i += 2 * CodazziDDD[j][k][l] * ntetU[k] * nn * (-remtetU[j]*immtetU[l]-immtetU[j]*remtetU[l])\n",
    "            psi3r += 1 * CodazziDDD[j][k][l] * nn * ((ntetU[j]-ltetU[j])*remtetU[k]*ntetU[l]-remtetU[j]*ltetU[k]*ntetU[l])\n",
    "            psi3i +=-1 * CodazziDDD[j][k][l] * nn * ((ntetU[j]-ltetU[j])*immtetU[k]*ntetU[l]-immtetU[j]*ltetU[k]*ntetU[l])\n",
    "            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]))\n",
    "            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]))\n",
    "            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])\n",
    "            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])\n",
    "            psi0r += 2 * CodazziDDD[j][k][l] * nn * ltetU[k]*(remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])\n",
    "            psi0i += 2 * CodazziDDD[j][k][l] * nn * ltetU[k]*(remtetU[j]*immtetU[l]+immtetU[j]*remtetU[l])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we will add the contractions of `GaussDDDD[i][j][k][l]` (from the Riemann tensor and extrinsic curvature, above) with the null tetrad.\n",
    "\n",
    "The relevant terms here are:\n",
    "\\begin{align}\n",
    "\\psi_4:&\\ (\\text{GaussDDDD[i][j][k][l]}) n^i \\overset{*}{m}{}^j n^k \\overset{*}{m}{}^l \\\\\n",
    "\\psi_3:&\\ (\\text{GaussDDDD[i][j][k][l]}) l^i n^j \\overset{*}{m}{}^k n^l \\\\\n",
    "\\psi_2:&\\ (\\text{GaussDDDD[i][j][k][l]}) l^i m^j \\overset{*}{m}{}^k n^l \\\\\n",
    "\\psi_1:&\\ (\\text{GaussDDDD[i][j][k][l]}) n^i l^j m^k l^l \\\\\n",
    "\\psi_0:&\\ (\\text{GaussDDDD[i][j][k][l]}) l^i m^j l^k m^l. \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:24.969298Z",
     "iopub.status.busy": "2021-03-07T17:20:24.933358Z",
     "iopub.status.idle": "2021-03-07T17:20:25.649156Z",
     "shell.execute_reply": "2021-03-07T17:20:25.648575Z"
    }
   },
   "outputs": [],
   "source": [
    "for l in range(3):\n",
    "    for j in range(3):\n",
    "        for k in range(3):\n",
    "            for i in range(3):\n",
    "                psi4r += GaussDDDD[i][j][k][l] * ntetU[i] * ntetU[k] * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])\n",
    "                psi4i += GaussDDDD[i][j][k][l] * ntetU[i] * ntetU[k] * (-remtetU[j]*immtetU[l]-immtetU[j]*remtetU[l])\n",
    "                psi3r += GaussDDDD[i][j][k][l] * ltetU[i] * ntetU[j] * remtetU[k] * ntetU[l]\n",
    "                psi3i +=-GaussDDDD[i][j][k][l] * ltetU[i] * ntetU[j] * immtetU[k] * ntetU[l]\n",
    "                psi2r += GaussDDDD[i][j][k][l] * ltetU[i] * ntetU[l] * (remtetU[j]*remtetU[k]+immtetU[j]*immtetU[k])\n",
    "                psi2i += GaussDDDD[i][j][k][l] * ltetU[i] * ntetU[l] * (immtetU[j]*remtetU[k]-remtetU[j]*immtetU[k])\n",
    "                psi1r += GaussDDDD[i][j][k][l] * ntetU[i] * ltetU[j] * remtetU[k] * ltetU[l]\n",
    "                psi1i += GaussDDDD[i][j][k][l] * ntetU[i] * ltetU[j] * immtetU[k] * ltetU[l]\n",
    "                psi0r += GaussDDDD[i][j][k][l] * ltetU[i] * ltetU[k] * (remtetU[j]*remtetU[l]-immtetU[j]*immtetU[l])\n",
    "                psi0i += GaussDDDD[i][j][k][l] * ltetU[i] * ltetU[k] * (remtetU[j]*immtetU[l]+immtetU[j]*remtetU[l])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation1'></a>\n",
    "\n",
    "## Step 5.a: Code Validation against `WeylScal4NRPy.WeylScalars_Cartesian` NRPy+ Module  \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation1}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for Weyl invariants between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [WeylScal4NRPy.WeylScalars_Cartesian](../edit/WeylScal4NRPy/WeylScalars_Cartesian.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:25.656888Z",
     "iopub.status.busy": "2021-03-07T17:20:25.656157Z",
     "iopub.status.idle": "2021-03-07T17:20:30.434523Z",
     "shell.execute_reply": "2021-03-07T17:20:30.435063Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between WeylScalars_Cartesian tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n",
      "psi4r - weyl.psi4r = 0\n",
      "psi4i - weyl.psi4i = 0\n",
      "psi3r - weyl.psi3r = 0\n",
      "psi3i - weyl.psi3i = 0\n",
      "psi2r - weyl.psi2r = 0\n",
      "psi2i - weyl.psi2i = 0\n",
      "psi1r - weyl.psi1r = 0\n",
      "psi1i - weyl.psi1i = 0\n",
      "psi0r - weyl.psi0r = 0\n",
      "psi0i - weyl.psi0i = 0\n"
     ]
    }
   ],
   "source": [
    "#psi4rb,psi4ib,psi3rb,psi3ib,psi2rb,psi2ib,psi1rb,psi1ib,psi0rb,psi0ib = psi4r,psi4i,psi3r,psi3i,psi2r,psi2i,psi1r,psi1i,psi0r,psi0i\n",
    "gri.glb_gridfcs_list = []\n",
    "import WeylScal4NRPy.WeylScalars_Cartesian as weyl\n",
    "par.set_parval_from_str(\"WeylScal4NRPy.WeylScalars_Cartesian::output_scalars\",\"all_psis\")\n",
    "weyl.WeylScalars_Cartesian()\n",
    "\n",
    "print(\"Consistency check between WeylScalars_Cartesian tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n",
    "\n",
    "print(\"psi4r - weyl.psi4r = \" + str(psi4r - weyl.psi4r))\n",
    "print(\"psi4i - weyl.psi4i = \" + str(psi4i - weyl.psi4i))\n",
    "print(\"psi3r - weyl.psi3r = \" + str(psi3r - weyl.psi3r))\n",
    "print(\"psi3i - weyl.psi3i = \" + str(psi3i - weyl.psi3i))\n",
    "print(\"psi2r - weyl.psi2r = \" + str(psi2r - weyl.psi2r))\n",
    "print(\"psi2i - weyl.psi2i = \" + str(psi2i - weyl.psi2i))\n",
    "print(\"psi1r - weyl.psi1r = \" + str(psi1r - weyl.psi1r))\n",
    "print(\"psi1i - weyl.psi1i = \" + str(psi1i - weyl.psi1i))\n",
    "print(\"psi0r - weyl.psi0r = \" + str(psi0r - weyl.psi0r))\n",
    "print(\"psi0i - weyl.psi0i = \" + str(psi0i - weyl.psi0i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='invariant_scalars'></a>\n",
    "\n",
    "# Step 6: The Invariant Scalars \\[Back to [top](#toc)\\]\n",
    "$$\\label{invariant_scalars}$$\n",
    "\n",
    "We may also wish to compute the invariant scalars, whose value does not depend on the choice of the null tetrad. While they are defined using the Weyl tensor, they can also be expressed in terms of the Weyl scalars. We will use those expressions for simplicity.\n",
    "\n",
    "Following after the method used in the Kranc code, we will read in the already-computed values of the Weyl scalars to find the invariants instead of trying to make NRPy output a very large expression in terms of the metric and extrinsic curvature.\n",
    "\n",
    "We will start with the invariants $I$ and $J$, as defined in equations (2.3a) and (2.3b) of [arXiv:gr-qc/0407013](https://arxiv.org/abs/gr-qc/0407013). They are as follows.\n",
    "\\begin{align}\n",
    "I &= 3 \\psi_2^2 - 4 \\psi_1 \\psi_3 + \\psi_4 \\psi_0 \\\\\n",
    "J &= \n",
    "\\begin{vmatrix}\n",
    "\\psi_4 & \\psi_3 & \\psi_2 \\\\\n",
    "\\psi_3 & \\psi_2 & \\psi_1 \\\\\n",
    "\\psi_2 & \\psi_1 & \\psi_0 \\\\\n",
    "\\end{vmatrix}\n",
    "\\end{align}\n",
    "Here, since we can work in terms of the Weyl scalars themselves, we will use SymPy's built-in tools for handling complex numbers, which will not get overwhelmed as they did when computing the Weyl scalars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:30.509912Z",
     "iopub.status.busy": "2021-03-07T17:20:30.473793Z",
     "iopub.status.idle": "2021-03-07T17:20:30.576881Z",
     "shell.execute_reply": "2021-03-07T17:20:30.577489Z"
    }
   },
   "outputs": [],
   "source": [
    "gri.glb_gridfcs_list = []\n",
    "psi4r,psi4i,psi3r,psi3i,psi2r,psi2i,psi1r,psi1i,psi0r,psi0i = gri.register_gridfunctions(\"AUX\",[\"psi4r\",\"psi4i\",\n",
    "                                                                                                \"psi3r\",\"psi3i\",\n",
    "                                                                                                \"psi2r\",\"psi2i\",\n",
    "                                                                                                \"psi1r\",\"psi1i\",\n",
    "                                                                                                \"psi0r\",\"psi0i\"])\n",
    "\n",
    "psi4 = psi4r + sp.I * psi4i\n",
    "psi3 = psi3r + sp.I * psi3i\n",
    "psi2 = psi2r + sp.I * psi2i\n",
    "psi1 = psi1r + sp.I * psi1i\n",
    "psi0 = psi0r + sp.I * psi0i\n",
    "\n",
    "curvIr = sp.re(3*psi2*psi2 - 4*psi1*psi3 + psi4*psi0)\n",
    "curvIi = sp.im(3*psi2*psi2 - 4*psi1*psi3 + psi4*psi0)\n",
    "curvJr = sp.re(psi4 * (psi2*psi0 - psi1*psi1) - \\\n",
    "               psi3 * (psi3*psi0 - psi1*psi2) +\\\n",
    "               psi2 * (psi3*psi1 - psi2*psi2) )\n",
    "curvJi = sp.im(psi4 * (psi2*psi0 - psi1*psi1) - \\\n",
    "               psi3 * (psi3*psi0 - psi1*psi2) +\\\n",
    "               psi2 * (psi3*psi1 - psi2*psi2) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now code the invariants $J_1$, $J_2$, $J_3$, and $J_4$, as found in equations B5-B8 of [arXiv:0704.1756](https://arxiv.org/abs/0704.1756). As with the other invariants, we will simply read in the values of the gridfunctions that they already calculated (that is, the Weyl scalars). These equations are based directly on those used in the Mathematica notebook that generates `WeylScal4` (available at [this](https://bitbucket.org/einsteintoolkit/einsteinanalysis/src) repository), modified so that Python can interpret them. Those equations were generated in turn using `xTensor` from equations B5-B8."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:30.687426Z",
     "iopub.status.busy": "2021-03-07T17:20:30.653187Z",
     "iopub.status.idle": "2021-03-07T17:20:30.952947Z",
     "shell.execute_reply": "2021-03-07T17:20:30.952189Z"
    }
   },
   "outputs": [],
   "source": [
    "J1curv =-16*(3*psi2i**2-3*psi2r**2-4*psi1i*psi3i+4*psi1r*psi3r+psi0i*psi4i-psi0r*psi4r)\n",
    "\n",
    "J2curv = 96*(-3*psi2i**2*psi2r+psi2r**3+2*psi1r*psi2i*psi3i+2*psi1i*psi2r*psi3i-psi0r*psi3i**2+\n",
    "            2*psi1i*psi2i*psi3r-2*psi1r*psi2r*psi3r-2*psi0i*psi3i*psi3r+psi0r*psi3r**2-\n",
    "            2*psi1i*psi1r*psi4i+psi0r*psi2i*psi4i+psi0i*psi2r*psi4i-psi1i**2*psi4r+psi1r**2*psi4r+\n",
    "            psi0i*psi2i*psi4r-psi0r*psi2r*psi4r)\n",
    "\n",
    "J3curv = 64*(9*psi2i**4-54*psi2i**2*psi2r**2+9*psi2r**4-24*psi1i*psi2i**2*psi3i+48*psi1r*psi2i*psi2r*psi3i+\n",
    "            24*psi1i*psi2r**2*psi3i+16*psi1i**2*psi3i**2-16*psi1r**2*psi3i**2+\n",
    "            24*psi1r*psi2i**2*psi3r+48*psi1i*psi2i*psi2r*psi3r-24*psi1r*psi2r**2*psi3r-64*psi1i*psi1r*psi3i*psi3r-\n",
    "            16*psi1i**2*psi3r**2+16*psi1r**2*psi3r**2+6*psi0i*psi2i**2*psi4i-12*psi0r*psi2i*psi2r*psi4i-\n",
    "            6*psi0i*psi2r**2*psi4i-8*psi0i*psi1i*psi3i*psi4i+8*psi0r*psi1r*psi3i*psi4i+8*psi0r*psi1i*psi3r*psi4i+\n",
    "            8*psi0i*psi1r*psi3r*psi4i+psi0i**2*psi4i**2-psi0r**2*psi4i**2-6*psi0r*psi2i**2*psi4r-\n",
    "            12*psi0i*psi2i*psi2r*psi4r+6*psi0r*psi2r**2*psi4r+8*psi0r*psi1i*psi3i*psi4r+8*psi0i*psi1r*psi3i*psi4r+\n",
    "            8*psi0i*psi1i*psi3r*psi4r-8*psi0r*psi1r*psi3r*psi4r-4*psi0i*psi0r*psi4i*psi4r-psi0i**2*psi4r**2+\n",
    "            psi0r**2*psi4r**2)\n",
    "\n",
    "J4curv = -640*(-15*psi2i**4*psi2r+30*psi2i**2*psi2r**3-3*psi2r**5+10*psi1r*psi2i**3*psi3i+\n",
    "              30*psi1i*psi2i**2*psi2r*psi3i-30*psi1r*psi2i*psi2r**2*psi3i-10*psi1i*psi2r**3*psi3i-\n",
    "              16*psi1i*psi1r*psi2i*psi3i**2-3*psi0r*psi2i**2*psi3i**2-8*psi1i**2*psi2r*psi3i**2+\n",
    "              8*psi1r**2*psi2r*psi3i**2-6*psi0i*psi2i*psi2r*psi3i**2+3*psi0r*psi2r**2*psi3i**2+\n",
    "              4*psi0r*psi1i*psi3i**3+4*psi0i*psi1r*psi3i**3+10*psi1i*psi2i**3*psi3r-\n",
    "              30*psi1r*psi2i**2*psi2r*psi3r-30*psi1i*psi2i*psi2r**2*psi3r+10*psi1r*psi2r**3*psi3r-\n",
    "              16*psi1i**2*psi2i*psi3i*psi3r+16*psi1r**2*psi2i*psi3i*psi3r-6*psi0i*psi2i**2*psi3i*psi3r+\n",
    "              32*psi1i*psi1r*psi2r*psi3i*psi3r+12*psi0r*psi2i*psi2r*psi3i*psi3r+6*psi0i*psi2r**2*psi3i*psi3r+\n",
    "              12*psi0i*psi1i*psi3i**2*psi3r-12*psi0r*psi1r*psi3i**2*psi3r+16*psi1i*psi1r*psi2i*psi3r**2+\n",
    "              3*psi0r*psi2i**2*psi3r**2+8*psi1i**2*psi2r*psi3r**2-8*psi1r**2*psi2r*psi3r**2+\n",
    "              6*psi0i*psi2i*psi2r*psi3r**2-3*psi0r*psi2r**2*psi3r**2-12*psi0r*psi1i*psi3i*psi3r**2-\n",
    "              12*psi0i*psi1r*psi3i*psi3r**2-4*psi0i*psi1i*psi3r**3+4*psi0r*psi1r*psi3r**3-\n",
    "              6*psi1i*psi1r*psi2i**2*psi4i+2*psi0r*psi2i**3*psi4i-6*psi1i**2*psi2i*psi2r*psi4i+\n",
    "              6*psi1r**2*psi2i*psi2r*psi4i+6*psi0i*psi2i**2*psi2r*psi4i+6*psi1i*psi1r*psi2r**2*psi4i-\n",
    "              6*psi0r*psi2i*psi2r**2*psi4i-2*psi0i*psi2r**3*psi4i+12*psi1i**2*psi1r*psi3i*psi4i-\n",
    "              4*psi1r**3*psi3i*psi4i-2*psi0r*psi1i*psi2i*psi3i*psi4i-2*psi0i*psi1r*psi2i*psi3i*psi4i-\n",
    "              2*psi0i*psi1i*psi2r*psi3i*psi4i+2*psi0r*psi1r*psi2r*psi3i*psi4i-2*psi0i*psi0r*psi3i**2*psi4i+\n",
    "              4*psi1i**3*psi3r*psi4i-12*psi1i*psi1r**2*psi3r*psi4i-2*psi0i*psi1i*psi2i*psi3r*psi4i+\n",
    "              2*psi0r*psi1r*psi2i*psi3r*psi4i+2*psi0r*psi1i*psi2r*psi3r*psi4i+2*psi0i*psi1r*psi2r*psi3r*psi4i-\n",
    "              2*psi0i**2*psi3i*psi3r*psi4i+2*psi0r**2*psi3i*psi3r*psi4i+2*psi0i*psi0r*psi3r**2*psi4i-\n",
    "              psi0r*psi1i**2*psi4i**2-2*psi0i*psi1i*psi1r*psi4i**2+psi0r*psi1r**2*psi4i**2+\n",
    "              2*psi0i*psi0r*psi2i*psi4i**2+psi0i**2*psi2r*psi4i**2-psi0r**2*psi2r*psi4i**2-3*psi1i**2*psi2i**2*psi4r+\n",
    "              3*psi1r**2*psi2i**2*psi4r+2*psi0i*psi2i**3*psi4r+12*psi1i*psi1r*psi2i*psi2r*psi4r-\n",
    "              6*psi0r*psi2i**2*psi2r*psi4r+3*psi1i**2*psi2r**2*psi4r-3*psi1r**2*psi2r**2*psi4r-\n",
    "              6*psi0i*psi2i*psi2r**2*psi4r+2*psi0r*psi2r**3*psi4r+4*psi1i**3*psi3i*psi4r-12*psi1i*psi1r**2*psi3i*psi4r-\n",
    "              2*psi0i*psi1i*psi2i*psi3i*psi4r+2*psi0r*psi1r*psi2i*psi3i*psi4r+2*psi0r*psi1i*psi2r*psi3i*psi4r+\n",
    "              2*psi0i*psi1r*psi2r*psi3i*psi4r-psi0i**2*psi3i**2*psi4r+psi0r**2*psi3i**2*psi4r-\n",
    "              12*psi1i**2*psi1r*psi3r*psi4r+4*psi1r**3*psi3r*psi4r+2*psi0r*psi1i*psi2i*psi3r*psi4r+\n",
    "              2*psi0i*psi1r*psi2i*psi3r*psi4r+2*psi0i*psi1i*psi2r*psi3r*psi4r-2*psi0r*psi1r*psi2r*psi3r*psi4r+\n",
    "              4*psi0i*psi0r*psi3i*psi3r*psi4r+psi0i**2*psi3r**2*psi4r-psi0r**2*psi3r**2*psi4r-\n",
    "              2*psi0i*psi1i**2*psi4i*psi4r+4*psi0r*psi1i*psi1r*psi4i*psi4r+2*psi0i*psi1r**2*psi4i*psi4r+\n",
    "              2*psi0i**2*psi2i*psi4i*psi4r-2*psi0r**2*psi2i*psi4i*psi4r-4*psi0i*psi0r*psi2r*psi4i*psi4r+\n",
    "              psi0r*psi1i**2*psi4r**2+2*psi0i*psi1i*psi1r*psi4r**2-psi0r*psi1r**2*psi4r**2-\n",
    "              2*psi0i*psi0r*psi2i*psi4r**2-psi0i**2*psi2r*psi4r**2+psi0r**2*psi2r*psi4r**2)\n",
    "\n",
    "#cse_output = sp.cse(psi0i,sp.numbered_symbols(\"tmp\"))\n",
    "#for commonsubexpression in cse_output:\n",
    "#    print(\"hello?\",commonsubexpression)\n",
    "#for commonsubexpression in cse_output[0]:\n",
    "#    print((str(commonsubexpression[0])+\" = \"+str(commonsubexpression[1])+\";\").replace(\"**\",\"^\").replace(\"_d\",\"d\"))\n",
    "#for i,result in enumerate(cse_output[1]):\n",
    "#   print((\"psi0iPy = \"+str(result)+\";\").replace(\"**\",\"^\").replace(\"_d\",\"d\"))\n",
    "# These replace commands are used to allow us to validate against Einstein Toolkit's WeylScal4 thorn in Mathematica.\n",
    "# Specifically, the first changes exponentiation to Mathematica's format, and the second strips the underscores\n",
    "# that have a very specific meaning in Mathematica and thus cannot be used in variable names."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation2'></a>\n",
    "\n",
    "## Step 6.a: Code Validation against `WeylScal4NRPy.WeylScalarInvariants_Cartesian` NRPy+ Module  \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation2}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for Weyl invariants between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [WeylScal4NRPy.WeylScalarInvariants_Cartesian](../edit/WeylScal4NRPy/WeylScalarInvariants_Cartesian.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:30.985206Z",
     "iopub.status.busy": "2021-03-07T17:20:30.959156Z",
     "iopub.status.idle": "2021-03-07T17:20:30.988926Z",
     "shell.execute_reply": "2021-03-07T17:20:30.988380Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between ScalarInvariants_Cartesian tutorial and NRPy+ module for invariant scalars: PASSED.\n"
     ]
    }
   ],
   "source": [
    "# Reset the list of gridfunctions, as registering a gridfunction\n",
    "#   twice will spawn an error.\n",
    "gri.glb_gridfcs_list = []\n",
    "\n",
    "import WeylScal4NRPy.WeylScalarInvariants_Cartesian as invar\n",
    "invar.WeylScalarInvariants_Cartesian()\n",
    "\n",
    "num_failures = 0\n",
    "if curvIr - invar.curvIr != 0: num_failures += 1\n",
    "if curvIi - invar.curvIi != 0: num_failures += 1\n",
    "if curvJr - invar.curvJr != 0: num_failures += 1\n",
    "if curvJi - invar.curvJi != 0: num_failures += 1\n",
    "if J1curv - invar.J1curv != 0: num_failures += 1\n",
    "if J2curv - invar.J2curv != 0: num_failures += 1\n",
    "if J3curv - invar.J3curv != 0: num_failures += 1\n",
    "if J4curv - invar.J4curv != 0: num_failures += 1\n",
    "\n",
    "import sys\n",
    "if num_failures == 0:\n",
    "    print(\"ScalarInvariants_Cartesian tutorial vs NRPy+ module: TESTS PASSED.\")\n",
    "else:\n",
    "    print(\"ScalarInvariants_Cartesian tutorial vs NRPy+ module: TESTS FAILED, WITH \"+str(num_failures)+\" FAILURES\")\n",
    "    sys.exit(1)"
   ]
  },
  {
   "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}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-WeylScalarsInvariants-Cartesian.pdf](Tutorial-WeylScalarsInvariants-Cartesian.pdf). (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:20:30.993547Z",
     "iopub.status.busy": "2021-03-07T17:20:30.992805Z",
     "iopub.status.idle": "2021-03-07T17:20:35.491201Z",
     "shell.execute_reply": "2021-03-07T17:20:35.490407Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-WeylScalarsInvariants-Cartesian.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-WeylScalarsInvariants-Cartesian.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-WeylScalarsInvariants-Cartesian\")"
   ]
  }
 ],
 "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
}