{
"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
}