Tutorial-SpinWeighted_Spherical_Harmonics.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",
    "# Calculating Spin-Weighted Spherical Harmonics\n",
    "## Authors: Zach Etienne & Brandon Clark\n",
    "\n",
    "[comment]: <> (Abstract: TODO)\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). In addition, its results have been validated against a [trusted Mathematica notebook](https://demonstrations.wolfram.com/versions/source.jsp?id=SpinWeightedSphericalHarmonics&version=0012).\n",
    "\n",
    "### NRPy+ Source Code for this module: [SpinWeight_minus2_SphHarmonics/SpinWeight_minus2_SphHarmonics.py](../edit/SpinWeight_minus2_SphHarmonics/SpinWeight_minus2_SphHarmonics.py)\n",
    "\n",
    "## Introduction:\n",
    "This tutorial notebook defines a Python function for computing spin-weighted spherical harmonics using Sympy. Spin-weight $s=-2$ spherical harmonics are the natural basis for decomposing gravitational wave data.\n",
    "\n",
    "The tutorial contains code necessary to validate the resulting expressions assuming $s=-2$ against a trusted Mathematica notebook (validated for all $(\\ell,m)$ up to $\\ell=8$. Finally it outputs a C code capable of computing $_{-2}Y_{\\ell m} (\\theta, \\phi)$ for all $(\\ell,m)$ for $\\ell=0$ up to `maximum_l`."
   ]
  },
  {
   "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): Initialize needed Python/NRPy+ modules\n",
    "1. [Step 2](#gbf): Defining the Goldberg function\n",
    "1. [Step 3](#math_code_validation): Code Validation against Mathematica script\n",
    "1. [Step 4](#ccode): Generate C-code function for computing s=-2 spin-weighted spherical harmonics, using NRPy+\n",
    "1. [Step 5](#code_validation): Code Validation against SpinWeight_minus2_SphHarmonics/SpinWeight_minus2_SphHarmonics NRPy+ module\n",
    "1. [Step 6](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize needed Python/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: Initialize needed Python/NRPy+ modules\n",
    "from outputC import outputC   # NRPy+: Core C code output module\n",
    "import sympy as sp            # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import os, sys                # Standard Python modules for multiplatform OS-level functions\n",
    "\n",
    "# Step 1.a: Set maximum l to which we will validate the spin-weighted spherical harmonics with s=-2:\n",
    "maximum_l = 4 # Note that we have validated against Mathematica up to and including l=8 -- perfect agreement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='gbf'></a>\n",
    "\n",
    "# Step 2: Defining the Goldberg function [Back to [top](#toc)\\]\n",
    "$$\\label{gbf}$$\n",
    "\n",
    "One way to calculate the spin-weighted spherical harmonics is using the following formula\n",
    "from [Goldberg et al. (1967)](https://aip.scitation.org/doi/10.1063/1.1705135):\n",
    "\n",
    "$$ _sY_{\\ell m} (\\theta, \\phi) = \\left(-1\\right)^m \\sqrt{ \\frac{(\\ell+m)! (\\ell-m)! (2\\ell+1)} {4\\pi (\\ell+s)! (\\ell-s)!} } \\sin^{2\\ell} \\left( \\frac{\\theta}{2} \\right) \\times\\sum_{r=0}^{\\ell-s} {\\ell-s \\choose r} {\\ell+s \\choose r+s-m} \\left(-1\\right)^{\\ell-r-s} e^{i m \\phi} \\cot^{2r+s-m} \\left( \\frac{\\theta} {2} \\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2: Defining the Goldberg function\n",
    "\n",
    "# Step 2.a: Declare SymPy symbols:\n",
    "th, ph = sp.symbols('th ph',real=True)\n",
    "\n",
    "# Step 2.b: Define the Goldberg formula for spin-weighted spherical harmonics\n",
    "#           (https://aip.scitation.org/doi/10.1063/1.1705135);\n",
    "#           referenced & described in Wikipedia Spin-weighted spherical harmonics article:\n",
    "#           https://en.wikipedia.org/w/index.php?title=Spin-weighted_spherical_harmonics&oldid=853425244\n",
    "def Y(s, l, m, th, ph, GenerateMathematicaCode=False):\n",
    "    Sum = 0\n",
    "    for r in range(l-s + 1):\n",
    "        if GenerateMathematicaCode == True:\n",
    "            # Mathematica needs expression to be in terms of cotangent, so that code validation below\n",
    "            #    yields identity with existing Mathematica notebook on spin-weighted spherical harmonics.\n",
    "            Sum +=  sp.binomial(l-s, r)*sp.binomial(l+s, r+s-m)*(-1)**(l-r-s)*sp.exp(sp.I*m*ph)*sp.cot(th/2)**(2*r+s-m)\n",
    "        else:\n",
    "            # SymPy C code generation cannot handle the cotangent function, so define cot(th/2) as 1/tan(th/2):\n",
    "            Sum +=  sp.binomial(l-s, r)*sp.binomial(l+s, r+s-m)*(-1)**(l-r-s)*sp.exp(sp.I*m*ph)/sp.tan(th/2)**(2*r+s-m)\n",
    "\n",
    "    return (-1)**m*sp.simplify(sp.sqrt(sp.factorial(l+m)*sp.factorial(l-m)*(2*l+1)/(4*sp.pi*sp.factorial(l+s)*sp.factorial(l-s)))*sp.sin(th/2)**(2*l)*Sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='math_code_validation'></a>\n",
    "\n",
    "# Step 3: Code Validation against Mathematica script [Back to [top](#toc)\\]\n",
    "$$\\label{math_code_validation}$$\n",
    "\n",
    "To validate the code we wish to compare it with an existent [Mathematica notebook](https://demonstrations.wolfram.com/versions/source.jsp?id=SpinWeightedSphericalHarmonics&version=0012). We will validate the code using a spin-value of $s=-2$ and $\\ell = 8,7,6,5,4,3,2,1,0$ while leaving $m$, $\\theta$, and $\\phi$ unknown. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Step 3: Code Validation against Mathematica notebook:\n",
    "#         https://demonstrations.wolfram.com/versions/source.jsp?id=SpinWeightedSphericalHarmonics&version=0012\n",
    "\n",
    "# # For the l=0 case m=0, otherwise there is a divide-by-zero in the Y() function above.\n",
    "# print(\"FullSimplify[Y[-2, 0, 0, th, ph]-\"+str(sp.mathematica_code(sp.simplify(Y(-2, 0, 0, th, ph,GenerateMathematicaCode=True))))+\"] \\n\") # Agrees with Mathematica notebook for l = 0\n",
    "\n",
    "# # Check the other cases\n",
    "# for l in range(1,maximum_l+1): # Agrees with Mathematica notebook for  l = 1, 2, 4, 5, 6, 7, 8;\n",
    "#     print(\"FullSimplify[Y[-2, \"+str(l)+\", m, th, ph]-(\"+\n",
    "#           str(sp.mathematica_code(sp.simplify(Y(-2, l, m, th, ph, GenerateMathematicaCode=True)))).replace(\"binomial\",\"Binomial\").replace(\"factorial\",\"Factorial\")+\")] \\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ccode'></a>\n",
    "\n",
    "# Step 4: Generate C-code function for computing s=-2 spin-weighted spherical harmonics, using NRPy+ [Back to [top](#toc)\\]\n",
    "$$\\label{ccode}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4: Generating C Code function for computing\n",
    "#         s=-2 spin-weighted spherical harmonics,\n",
    "#         using NRPy+'s outputC() function.\n",
    "\n",
    "outCparams = \"preindent=3,outCfileaccess=a,outCverbose=False,includebraces=True\"\n",
    "\n",
    "with open(os.path.join(\"SpinWeight_minus2_SphHarmonics\",\"SpinWeight_minus2_SphHarmonics.h\"), \"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "void SpinWeight_minus2_SphHarmonics(const int l, const int m, const REAL th, const REAL ph,\n",
    "                                   REAL *reYlmswm2_l_m, REAL *imYlmswm2_l_m) {\n",
    "if(l<0 || l>\"\"\"+str(maximum_l)+\"\"\" || m<-l || m>+l) {\n",
    "    printf(\"ERROR: SpinWeight_minus2_SphHarmonics handles only l=[0,\"\"\"+str(maximum_l)+\"\"\"] and only m=[-l,+l] is defined.\\\\n\");\n",
    "    printf(\"       You chose l=%d and m=%d, which is out of these bounds.\\\\n\",l,m);\n",
    "    exit(1);\n",
    "}\\n\"\"\")\n",
    "\n",
    "    file.write(\"switch(l) {\\n\")\n",
    "    for l in range(maximum_l+1): # Output values up to and including l=8.\n",
    "        file.write(\"    case \"+str(l)+\":\\n\")\n",
    "        file.write(\"        switch(m) {\\n\")\n",
    "        for m in range(-l,l+1):\n",
    "            file.write(\"            case \"+str(m)+\":\\n\")\n",
    "            Y_m2_lm = Y(-2, l, m, th, ph)\n",
    "            Cstring = outputC([sp.re(Y_m2_lm),sp.im(Y_m2_lm)],[\"*reYlmswm2_l_m\",\"*imYlmswm2_l_m\"],\n",
    "                              \"returnstring\",outCparams)\n",
    "            file.write(Cstring)\n",
    "            file.write(\"                  return;\\n\")\n",
    "        file.write(\"        }  /* End switch(m) */\\n\")\n",
    "    file.write(\"    } /* End switch(l) */\\n\")\n",
    "    file.write(\"} /* End function SpinWeight_minus2_SphHarmonics() */\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# [Step 5](#code_validation): Code Validation against `SpinWeight_minus2_SphHarmonics.SpinWeight_minus2_SphHarmonics` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "As additional validation, we verify agreement in the SymPy expressions for the spin-weight -2 spherical harmonics expressions between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [`SpinWeight_minus2_SphHarmonics.SpinWeight_minus2_SphHarmonics`](../edit/SpinWeight_minus2_SphHarmonics/SpinWeight_minus2_SphHarmonics.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "### BEGIN VALIDATION TESTS ###\n",
      "VALIDATION TEST PASSED on file: SpinWeight_minus2_SphHarmonics/SpinWeight_minus2_SphHarmonics.h\n",
      "### END VALIDATION TESTS ###\n"
     ]
    }
   ],
   "source": [
    "import SpinWeight_minus2_SphHarmonics.SpinWeight_minus2_SphHarmonics as swm2\n",
    "swm2.SpinWeight_minus2_SphHarmonics(maximum_l=4,filename=os.path.join(\"SpinWeight_minus2_SphHarmonics\",\"SpinWeight_minus2_SphHarmonics-NRPymodule.h\"))\n",
    "\n",
    "print(\"\\n\\n### BEGIN VALIDATION TESTS ###\")\n",
    "import filecmp\n",
    "fileprefix = os.path.join(\"SpinWeight_minus2_SphHarmonics\",\"SpinWeight_minus2_SphHarmonics\")\n",
    "if filecmp.cmp(fileprefix+\"-NRPymodule.h\",fileprefix+\".h\") == False:\n",
    "    print(\"VALIDATION TEST FAILED ON file: \"+fileprefix+\".h\"+\".\")\n",
    "    sys.exit(1)\n",
    "\n",
    "print(\"VALIDATION TEST PASSED on file: \"+fileprefix+\".h\")\n",
    "print(\"### END VALIDATION TESTS ###\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 6: Output this notebook to $\\LaTeX$-formatted PDF \\[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-SpinWeighted_Spherical_Harmonics.pdf](Tutorial-SpinWeighted_Spherical_Harmonics.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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-SpinWeighted_Spherical_Harmonics.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-SpinWeighted_Spherical_Harmonics.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-SpinWeighted_Spherical_Harmonics\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}