{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Calculating Spin-Weighted Spherical Harmonics\n", "## Authors: Zach Etienne & Brandon Clark\n", "\n", "## This notebook presents a Python code, designed to calculate spin-weighted spherical harmonics utilizing Sympy and the Goldberg function [Goldberg et al. (1967)](https://aip.scitation.org/doi/10.1063/1.1705135). The implementation is verified against a recognized Mathematica notebook and related NRPy+ module, and is then output as a C code.\n", "\n", "**Notebook Status:** Validated \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": [ "\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": [ "\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": { "execution": { "iopub.execute_input": "2021-03-07T17:25:18.808178Z", "iopub.status.busy": "2021-03-07T17:25:18.806632Z", "iopub.status.idle": "2021-03-07T17:25:19.145007Z", "shell.execute_reply": "2021-03-07T17:25:19.144211Z" } }, "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": [ "\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": { "execution": { "iopub.execute_input": "2021-03-07T17:25:19.154983Z", "iopub.status.busy": "2021-03-07T17:25:19.153138Z", "iopub.status.idle": "2021-03-07T17:25:19.157790Z", "shell.execute_reply": "2021-03-07T17:25:19.157186Z" } }, "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": [ "\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": { "execution": { "iopub.execute_input": "2021-03-07T17:25:19.162011Z", "iopub.status.busy": "2021-03-07T17:25:19.161346Z", "iopub.status.idle": "2021-03-07T17:25:19.163515Z", "shell.execute_reply": "2021-03-07T17:25:19.164031Z" }, "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": [ "\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": { "execution": { "iopub.execute_input": "2021-03-07T17:25:19.176134Z", "iopub.status.busy": "2021-03-07T17:25:19.173573Z", "iopub.status.idle": "2021-03-07T17:25:23.194116Z", "shell.execute_reply": "2021-03-07T17:25:23.193385Z" } }, "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": [ "\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": { "execution": { "iopub.execute_input": "2021-03-07T17:25:23.200726Z", "iopub.status.busy": "2021-03-07T17:25:23.199972Z", "iopub.status.idle": "2021-03-07T17:25:26.026593Z", "shell.execute_reply": "2021-03-07T17:25:26.027096Z" } }, "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": [ "\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": { "execution": { "iopub.execute_input": "2021-03-07T17:25:26.031970Z", "iopub.status.busy": "2021-03-07T17:25:26.031201Z", "iopub.status.idle": "2021-03-07T17:25:29.330415Z", "shell.execute_reply": "2021-03-07T17:25:29.329397Z" } }, "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 (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 }