### Install all recommended software

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Installs the package with all recommended dependencies.

```bash
pip install '.[dev, testing, extended]'
```

--------------------------------

### Install pybhpt with pip (developer)

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Installs the package from the source directory. Optionally, install with development and testing dependencies.

```bash
pip install .
```

```bash
pip install '.[dev, testing]'
```

--------------------------------

### Install target

Source: https://github.com/znasipak/pybhpt/blob/main/CMakeLists.txt

Installs the built cybhpt_full module to the current directory.

```cmake
# Optional install
install(TARGETS cybhpt_full DESTINATION "./")
```

--------------------------------

### Install pybhpt with pip

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Quick installation command for tagged releases of pybhpt from PyPI.

```bash
python3 -m pip install pybhpt
```

--------------------------------

### Install pybhpt with pip

Source: https://github.com/znasipak/pybhpt/blob/main/README.md

Installs the pybhpt package from the current directory.

```bash
pip install .
```

--------------------------------

### Install pybhpt with dev and testing extras

Source: https://github.com/znasipak/pybhpt/blob/main/README.md

Installs pybhpt with basic development and testing dependencies.

```bash
pip install '.[dev, testing]'
```

--------------------------------

### Install compiler for Linux

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Installs the necessary compiler for Linux using conda.

```bash
conda install gcc_linux-64 gxx_linux-64
```

--------------------------------

### Install compiler for macOS arm/silicon

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Installs the necessary compiler for macOS arm/silicon using conda.

```bash
conda install clang_osx-arm64 clangxx_osx-arm64
```

--------------------------------

### CMake project setup and dependency detection

Source: https://github.com/znasipak/pybhpt/blob/main/CMakeLists.txt

Initializes the CMake project with CXX language, finds Python, NumPy, and GSL, and adds Conda environment paths if CONDA_PREFIX is set.

```cmake
cmake_minimum_required(VERSION 3.15)
project(pybhpt LANGUAGES CXX)

# Python & NumPy
find_package(Python COMPONENTS Interpreter Development.Module NumPy REQUIRED)

# Conda detection
if(DEFINED ENV{CONDA_PREFIX})
  message(STATUS "Using Conda environment at $ENV{CONDA_PREFIX}")
  list(APPEND CMAKE_PREFIX_PATH "$ENV{CONDA_PREFIX}")
  list(APPEND CMAKE_LIBRARY_PATH "$ENV{CONDA_PREFIX}/lib")
  list(APPEND CMAKE_INCLUDE_PATH "$ENV{CONDA_PREFIX}/include")
endif()

# ------------------------- GSL Detection -------------------------
find_package(GSL REQUIRED)

# ------------------------- Boost Headers Only -------------------------
set(BOOST_INCLUDE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/extern/boost")
```

--------------------------------

### Install compiler for macOS Intel

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Installs the necessary compiler for macOS Intel using conda.

```bash
conda install clang_osx-64 clangxx_osx-64
```

--------------------------------

### Install Jupyter dependencies

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Installs ipykernel and matplotlib for running code in Jupyter notebooks.

```bash
pip install ipykernel matplotlib
```

--------------------------------

### Import required libraries

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/swsh.ipynb

Imports the pybhpt.swsh module, NumPy, and Matplotlib for the tutorial examples in this notebook.

```python
from pybhpt.swsh import SpinWeightedSpheroidalHarmonic, swsh_eigenvalue
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
```

--------------------------------

### Install Jupyter dependencies and register pybhpt-env kernel

Source: https://github.com/znasipak/pybhpt/blob/main/README.md

Install the required Jupyter dependencies and register the conda environment as a Jupyter kernel. The dependencies are already included in environment-extended.yml or the optional pip install dependencies.

```bash
pip install ipykernel matplotlib
```

```bash
python -m ipykernel install --user --name=pybhpt-env
```

--------------------------------

### Run unit tests with pytest

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Runs the unit tests to ensure the package installed properly.

```bash
pytest .
```

--------------------------------

### Flux summation example with output

Source: https://github.com/znasipak/pybhpt/blob/main/notebooks/flux-example.ipynb

Computes gravitational wave fluxes for a given spin, eccentricity, semi-latus rectum, and inclination using the FluxSummation class. Prints timing and the resulting dE/dt, dL/dt, and dQ/dt values.

```python
from datetime import datetime
spin = 0.9
e0 = 0.1
p0 = 6
x0 =  np.cos(40/180*np.pi)
print(p0, e0)
print(datetime.now().strftime("%H:%M:%S"))
summer = FluxSummation(spin, p0, e0, x0, tol=1e-10)
summer.sum(print_output_n_sum = True, print_output_k_sum = True, print_output_m_sum = True)
print(datetime.now().strftime("%H:%M:%S"))
print('dE/dt = ', 2.*summer.totalfluxes[0])
print('dL/dt = ', 2.*summer.totalfluxes[1])
print('dQ/dt = ', 2.*summer.totalfluxes[2])
```

```text
6 0.1
16:28:27
Flux reference: [0.00035882 0.00549491 0.02729088]
m summation results
7 1 0 0 [-7.97310493e-09 -1.22098052e-07 -1.63347770e-07]
m summation results
8 2 0 0 [0.00014164 0.00216904 0.00290183]
m summation results
9 3 0 0 [1.69899187e-05 2.60179190e-04 3.48078367e-04]
m summation results
10 4 0 0 [2.29180078e-06 3.50960401e-05 4.69529186e-05]
m summation results
10 5 0 0 [2.97788882e-07 4.56026137e-06 6.10090426e-06]
m summation results
12 6 0 0 [3.43094501e-08 5.25405982e-07 7.02909622e-07]
m summation results
12 7 0 0 [3.13083297e-09 4.79447607e-08 6.41424627e-08]
m summation results
13 8 0 0 [1.67457318e-10 2.56439775e-09 3.43075625e-09]
m summation results
12 9 0 0 [2.85699826e-13 4.37513272e-12 5.85323158e-12]
m summation results
14 10 0 0 [3.79494604e-12 5.81148152e-11 7.77483778e-11]
m summation results
15 11 0 0 [2.99842773e-12 4.59171413e-11 6.14298305e-11]
```

--------------------------------

### k-mode summation loop with convergence testing

Source: https://github.com/znasipak/pybhpt/blob/main/notebooks/flux-example.ipynb

Performs the sum over k-modes, starting at mpeak and increasing, then decreasing from mpeak-1. It records the max m-modes used, checks for peak k-mode, adds the k-mode to the n-mode sum, and tests convergence using convergenceTestModes.

```python
# perform sum with increasing m, starting at mpeak
self.kfluxes = np.zeros(3)
self.m_mode_summation(mpeak, increment = 1, include_zero = include_zero, print_output = print_output_m_sum)

# set minimum m value based on whether we include m = 0 for this value of (k,n)
mmin = 1
if include_zero:
    mmin = 0
self.mmaxmodenums[self.k] = [mmin, self.m - 1] # record the max m-modes used for the sum

# perform sum with decreasing m, starting at mpeak - 1 and stopping either
# when the sum converges or m = mmin
if mpeak > 0:
    self.m_mode_summation(mpeak - 1, increment = -1, include_zero = include_zero, print_output = print_output_m_sum)

# check to see where the k-modes are peaking as we sum over k
self.check_for_peak_k_mode()
if self.m + 1 < mmin:
    self.mmaxmodenums[self.k][0] = self.m + 1

# add the latest k-mode to the n-mode sum
self.nfluxes += self.kfluxes

# check if the k-modes are passing the convergence tests
convergedK = convergenceTestModes(self.kfluxes, self.fluxref, previous_kfluxes, previous_kfluxes_2, flux_position=self.flux_positions, tol=self.tol)
if convergedK is True:
    convergedCount += 1

# store most recent k-mode flux values
previous_kfluxes_2 = previous_kfluxes
previous_kfluxes = self.kfluxes

if print_output_k_sum:
    print("k summation results (mmax, k, n)")
    print(self.mmaxmodenums[self.k][1], self.k, self.n, self.kfluxes)

mpeak = self.mpeakmodenums[self.k]
self.k += increment
```

--------------------------------

### Link OpenMP and GSL libraries

Source: https://github.com/znasipak/pybhpt/blob/main/CMakeLists.txt

Links OpenMP if found, otherwise defines PYBHPT_NO_OPENMP. Links GSL libraries and sets RPATH for installation.

```cmake
# OpenMP
find_package(OpenMP)
if(OpenMP_CXX_FOUND)
  target_link_libraries(cybhpt_full PRIVATE OpenMP::OpenMP_CXX)
else()
  message(STATUS "OpenMP not found, proceeding without it")
  target_compile_definitions(cybhpt_full PRIVATE PYBHPT_NO_OPENMP)
endif()

# GSL
target_link_libraries(cybhpt_full PRIVATE GSL::gsl GSL::gslcblas)
set_target_properties(cybhpt_full PROPERTIES
  INSTALL_RPATH "$ORIGIN/../../opt/gsl-install/lib"
  BUILD_WITH_INSTALL_RPATH TRUE
)
```

--------------------------------

### Install Jupyter kernel for pybhpt-env

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Makes the conda environment accessible within Jupyter by installing a kernel named pybhpt-env.

```bash
python -m ipykernel install --user --name=pybhpt-env
```

--------------------------------

### Output of total fluxes

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/fluxes.ipynb

Shows the output of printing the total fluxes for the circular orbit example.

```text
Total fluxes (E, Lz, Q):  [0.00122957 0.01496423 0.        ]
```

--------------------------------

### Get Teukolsky amplitudes for In and Up

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/teuk.ipynb

Retrieves the Teukolsky amplitudes for the 'In' and 'Up' boundary conditions. The output shows the complex amplitudes.

```python
Psis.amplitude('In'), Psis.amplitude('Up')
```

```text
Result:
((-7.722257033945276e-16+2.1438981566064308e-15j),
 (-1.0298101162765425e-13-8.926429485261111e-13j))
```

--------------------------------

### Install conda compilers for pybhpt source compilation

Source: https://github.com/znasipak/pybhpt/blob/main/README.md

Install the necessary compiler into the conda environment to fix cmake compiler identification issues when compiling pybhpt from source. The commands are for macOS Intel, macOS arm/silicon, and Linux respectively.

```bash
conda install clang_osx-64 clangxx_osx-64
```

```bash
conda install clang_osx-arm64 clangxx_osx-arm64
```

```bash
conda install gcc_linux-64 gxx_linux-64
```

--------------------------------

### Output of flux printing

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/fluxes.ipynb

Shows the output of printing the horizon, infinity, and total fluxes for the example mode. The values are very small, as expected for a high-order mode.

```text
[-3.606222842872166e-36, -3.5075219557527026e-35, -1.979059439027494e-34]
[6.753687090778148e-25, 6.56884135710292e-24, 3.7063566971910574e-23]
[6.753687090742086e-25, 6.568841357067845e-24, 3.706356697171267e-23]
```

--------------------------------

### m-mode summation loop with adaptive sampling

Source: https://github.com/znasipak/pybhpt/blob/main/notebooks/flux-example.ipynb

Sums over m-modes starting from minit with a given increment. It increases the geometric sampling rate for high m-modes, sets flux positions based on m, and tests convergence using convergenceTestModes.

```python
def m_mode_summation(self, minit, increment = 1, include_zero = False, print_output = False):
    # increment sets how much to change m. Usually set to +1 or -1
    mmin = 1
    if include_zero:
        mmin = 0
    convergedM = False
    self.m = minit
    previous_mfluxes = np.zeros(3)
    previous_mfluxes_2 = np.zeros(3)
    while (self.m < self.mlimit and self.m >= mmin and convergedM is False):
        # increase sampling rate for high m-modes
        if 8*self.m > self.nsamples:
            self.increase_geo_sampling(2*self.nsamples)
        # perform sum over l-modes
        if self.m == 0:
            self.flux_positions = 0 # we ignore angular momentum for m = 0 because it vanishes
        else:
            self.flux_positions = [0, 1, 2]
            
        self.mfluxes = np.zeros(3)
        self.l_mode_summation()
        self.check_for_peak_m_mode()

        # add this m-mode contribution to n-mode and test for convergence
        self.kfluxes += self.mfluxes
        convergedM = convergenceTestModes(self.mfluxes, self.fluxref, previous_mfluxes, previous_mfluxes_2, flux_position=self.flux_positions, tol=self.tol)
        previous_mfluxes_2 = previous_mfluxes
        previous_mfluxes = self.mfluxes

        if print_output:
            print("m summation results")
            print(self.l - 1, self.m, self.k, self.n, self.mfluxes)
        self.m += increment

    if self.m < -1:
        print(self.m, self.k, self.n, mmin, include_zero)
    self.flux_positions = [0, 1, 2]
```

--------------------------------

### Activate Boost submodule

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Initializes and updates the Boost submodule under extern/boost.

```bash
git submodule update --init --recursive extern/boost
```

--------------------------------

### Build Boost headers

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Builds the Boost headers by running bootstrap.sh and b2 headers in the extern/boost directory.

```bash
cd extern/boost
chmod +x bootstrap.sh
./bootstrap.sh
./b2 headers
```

--------------------------------

### Import pybhpt modules and dependencies

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/waveform.ipynb

Imports TeukolskyMode, KerrGeodesic, SpinWeightedSpheroidalHarmonic, numpy, and matplotlib for waveform construction and plotting.

```python
from pybhpt.teuk import TeukolskyMode
from pybhpt.geo import KerrGeodesic
from pybhpt.swsh import SpinWeightedSpheroidalHarmonic
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
```

--------------------------------

### Initialize KerrGeodesic with parameters

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/tutorial.ipynb

Creates a KerrGeodesic object with spin a=0.9, semi-latus rectum p=8, eccentricity e=0.2, inclination x=0.9, and 2^9 samples.

```python
a, p, e, x, nsamples = (0.9, 8., 0.2, 0.9, 2**9)
geo = KerrGeodesic(a, p, e, x, nsamples)
```

--------------------------------

### Run all benchmarks

Source: https://github.com/znasipak/pybhpt/blob/main/docs/performance/index.md

Runs the full benchmark suite covering geodesic, swsh, radial, teuk, and flux computations. Use this command to reproduce or regenerate all committed performance data.

```bash
python -m benchmarks.run all          # geodesic, swsh, radial, teuk (+flux)
```

--------------------------------

### Get Mino-time frequencies

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/geodesics.ipynb

Retrieve the Mino-time frequencies of the geodesic. The output is an array of four values.

```python
geo.minofrequencies
```

--------------------------------

### Set include paths and C++ standard

Source: https://github.com/znasipak/pybhpt/blob/main/CMakeLists.txt

Adds NumPy, Boost, GSL, and project include directories. Sets C++14 standard and adds -march compiler option if supported.

```cmake
target_include_directories(cybhpt_full PRIVATE
  ${Python_NumPy_INCLUDE_DIRS}
  ${BOOST_INCLUDE_DIR}
  ${GSL_INCLUDE_DIRS}
  ${CMAKE_CURRENT_SOURCE_DIR}/cpp/include
)

# C++ flags
set_property(TARGET cybhpt_full PROPERTY CXX_STANDARD 14)
set(PYBHPT_MARCH native CACHE STRING "Value of the -march compiler option")
include(CheckCXXCompilerFlag)
set(MARCH_OPT "-march=${PYBHPT_MARCH}")
check_cxx_compiler_flag("${MARCH_OPT}" COMPILER_SUPPORTS_MARCH)
if(COMPILER_SUPPORTS_MARCH)
  target_compile_options(cybhpt_full PRIVATE "${MARCH_OPT}")
endif()
```

--------------------------------

### Get radial roots of the geodesic

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/geodesics.ipynb

Access the roots of the radial equation for the geodesic. The output is an array of four values.

```python
geo.radialroots
```

--------------------------------

### Run quick SWSh benchmark subset

Source: https://github.com/znasipak/pybhpt/blob/main/docs/performance/index.md

Runs a fast subset of the SWSh benchmarks, useful for quick checks during development. The --quick flag limits the parameter space to reduce runtime.

```bash
python -m benchmarks.run swsh --quick # fast subset
```

--------------------------------

### Get polar roots of the geodesic

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/geodesics.ipynb

Access the roots of the polar equation for the geodesic. The output is an array of two values.

```python
geo.polarroots
```

--------------------------------

### Run regression benchmark script

Source: https://github.com/znasipak/pybhpt/blob/main/benchmarks/regression/REPORT.md

Run the regression benchmark script to build each branch (main in a worktree; symlink extern/boost) and generate per-side summary statistics in summary.csv. Raw per-mode dumps (~28 MB each) are regenerable with this command and are not committed.

```bash
bash benchmarks/regression/run_regression.sh
```

--------------------------------

### Construct a KerrGeodesic with orbital parameters

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/geodesics.ipynb

Instantiate a KerrGeodesic object with spin a=0.9, semilatus rectum p=8.0, eccentricity e=0.6, inclination cosine x=0.9, and nsamples=512. The optional nsamples parameter controls the number of pre-sampled phase-space points stored for later computations.

```python
a, p, e, x, nsamples = (0.9, 8., 0.6, 0.9, 2**9)
geo = KerrGeodesic(a, p, e, x, nsamples)
```

--------------------------------

### Get fundamental frequencies and time frequencies

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/geodesics.ipynb

Access the fundamental frequencies with respect to coordinate time t. Both geo.frequencies and geo.timefrequencies return the same array of three values.

```python
geo.frequencies, geo.timefrequencies
```

--------------------------------

### Instantiate SpinWeightedSpheroidalHarmonic

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/swsh.ipynb

Instantiates the SpinWeightedSpheroidalHarmonic class with spin s, j, m, and gamma.

```python
s = -1
j = 4
m = 1
gamma = 0.43
Ssjm = SpinWeightedSpheroidalHarmonic(s, j, m, gamma)
```

--------------------------------

### Import pybhpt modules and list available gauges

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/tutorial.ipynb

Imports the main pybhpt classes and prints the available Hertz gauges. The output shows the list ['IRG', 'ORG', 'SRG0', 'SRG4', 'ARG0', 'ARG4'].

```python
from pybhpt.geo import KerrGeodesic
from pybhpt.teuk import TeukolskyMode
from pybhpt.hertz import HertzMode
from pybhpt.flux import FluxMode
from pybhpt.hertz import available_gauges
import numpy as np
print(available_gauges)
```

```text
['IRG', 'ORG', 'SRG0', 'SRG4', 'ARG0', 'ARG4']
```

--------------------------------

### Create conda environment with environment.yml

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Creates and activates a conda environment named pybhpt-env with necessary dependencies.

```bash
conda env create -f environment.yml
conda activate pybhpt-env
```

--------------------------------

### Generate all relevant modes

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/waveform.ipynb

Produces all relevant modes and stores them in the list `modes`. This can take several minutes depending on the tolerance and source nature. Uses n_mode_sum and checks convergence for each m and j.

```python
total_power_sq = 0.
nlimit = 120
klimit = 40

modes = []
jerror = 1
j = 2
jmax = 15
while (j <= jmax) and (jerror > tol):
    m = j
    j_total_power_sq = 0.
    mconvergence = 0
    while m >= 0 and mconvergence < 2:     
        print(j, m)

        m_total_power_sq = 0.
        modes, m_total_power_sq, total_power_sq = n_mode_sum(j, m, geo, modes, m_total_power_sq, total_power_sq, tol=tol, nlimit=nlimit, klimit=klimit)
        j_total_power_sq += m_total_power_sq

        merror = np.sqrt(m_total_power_sq/total_power_sq)
        if merror < tol:
            mconvergence += 1
        else:
            mconvergence = 0
        m -= 1
        # print(np.sqrt(m_total_power_sq))
    jerror = np.sqrt(j_total_power_sq/total_power_sq)
    # print(j, jerror, tol, jerror > tol, np.sqrt(total_power_sq))
    j += 1
```

--------------------------------

### Import TeukolskyMode and KerrGeodesic

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/teuk.ipynb

Imports the TeukolskyMode and KerrGeodesic classes along with numpy and matplotlib for the tutorial.

```python
from pybhpt.teuk import TeukolskyMode
from pybhpt.geo import KerrGeodesic
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
```

--------------------------------

### Import pybhpt and plotting libraries

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/fluxes.ipynb

Imports the TeukolskyMode, FluxMode, and KerrGeodesic classes from pybhpt, along with numpy and matplotlib for numerical and plotting work.

```python
from pybhpt.teuk import TeukolskyMode
from pybhpt.flux import FluxMode
from pybhpt.geo import KerrGeodesic
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
```

--------------------------------

### Specify time samples and viewing angles

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/waveform.ipynb

Sets the polar viewing angles th and phi, and creates a time grid u_grid spanning 4 orbital periods with 800 samples.

```python
th = 0.4
phi = 0.2
u_grid = np.linspace(0, 4*2*np.pi/geo.frequencies[0], 800)
```

--------------------------------

### Print flux values at infinity, horizon, and total

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/tutorial.ipynb

Prints the fluxes at infinity, horizon, and total for energy, angular momentum, and Carter constant. The output shows the computed values.

```python
print(f"[Edot^I, Ldot^I, Qdot^I] = {fluxes.infinityfluxes}")
print(f"[Edot^H, Ldot^H, Qdot^H] = {fluxes.horizonfluxes}")
print(f"[Edot, Ldot, Qdot] = {fluxes.totalfluxes}")
```

```text
[Edot^I, Ldot^I, Qdot^I] = [4.3416196240555994e-10, 4.190550504091079e-09, 1.5989065786113275e-08]
[Edot^H, Ldot^H, Qdot^H] = [-1.5687105316225854e-10, -1.5141263579703742e-09, -5.777156467256893e-09]
[Edot, Ldot, Qdot] = [2.772909092433014e-10, 2.6764241461207047e-09, 1.0211909318856381e-08]
```

--------------------------------

### Estimate amplitude precision for In and Up

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/teuk.ipynb

Estimates the fractional numerical error in the amplitudes for the 'In' and 'Up' boundary conditions. The output shows the error estimates.

```python
Psis.precision('In'), Psis.precision('Up')
```

```text
Result:
(8.635931505264635e-13, 5.361783459428641e-11)
```

--------------------------------

### Define Python extension module cybhpt_full

Source: https://github.com/znasipak/pybhpt/blob/main/CMakeLists.txt

Creates a Python module named cybhpt_full from the generated Cython wrapper and multiple C++ source files. Sets WITH_SOABI for Python ABI.

```cmake
python_add_library(cybhpt_full MODULE
  "${CMAKE_CURRENT_BINARY_DIR}/cython/redshift_wrap.cpp"
  cpp/src/utils.cpp
  cpp/src/specialfunc.cpp
  cpp/src/kerr.cpp
  cpp/src/geo.cpp
  cpp/src/cf.cpp
  cpp/src/swsh.cpp
  cpp/src/monodromy.cpp
  cpp/src/hypergeo_f.cpp
  cpp/src/bessel.cpp
  cpp/src/hypergeo_u.cpp
  cpp/src/nusolver.cpp
  cpp/src/mst.cpp
  cpp/src/gsn_asymp.cpp
  cpp/src/radialsolver.cpp
  cpp/src/sourceintegration.cpp
  cpp/src/teukolsky.cpp
  cpp/src/hertz.cpp
  cpp/src/metriccoeffs.cpp
  cpp/src/fluxes.cpp
  cpp/src/metric.cpp
  cpp/src/redshift.cpp
  WITH_SOABI
)
```

--------------------------------

### Import pybhpt.geo.KerrGeodesic and plotting libraries

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/geodesics.ipynb

Import the KerrGeodesic class from pybhpt.geo along with numpy and matplotlib for numerical and plotting operations.

```python
from pybhpt.geo import KerrGeodesic
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
```

--------------------------------

### Create RadialTeukolsky object with parameters

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/radial.ipynb

Set Teukolsky parameters (s, j, m, a, omega) and construct a radial grid from the event horizon plus 2 to 100 with 3000 points. Instantiates a RadialTeukolsky object with these parameters and the grid.

```python
s, j, m, a, omega = -2, 12, 3, 0.99, 2.2
r_hor = 1 + np.sqrt(1 - a**2)
r_grid = np.linspace(r_hor + 2, 100, 3000)
Rt = RadialTeukolsky(s, j, m, a, omega, r_grid)

```

--------------------------------

### Construct Hertz potentials for all gauges

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/tutorial.ipynb

Builds HertzMode objects for each available gauge, solves them, and evaluates the resulting potentials on a radial grid spanning from inside the inner turning point to r=30.

```python
rmin, rmax = geo.radialpoints[[0, -1]]
rinner = np.linspace(3., rmin - 0.001, 200)
rupper = np.linspace(rmax + 0.001, 30, 200)
r = np.concatenate((rinner, rupper))

phi = {}; PhiOfR = {}
for gauge in available_gauges:
    phi[gauge] = HertzMode(teuk, gauge)
    phi[gauge].solve()
    PhiOfR[gauge] = phi[gauge](r)
```

--------------------------------

### Clone pybhpt repository

Source: https://github.com/znasipak/pybhpt/blob/main/docs/installation.md

Clones the pybhpt repository from GitHub and changes into the directory.

```bash
git clone https://github.com/znasipak/pybhpt.git
cd pybhpt
```

--------------------------------

### Build waveform from mode amplitudes

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/waveform.ipynb

Constructs the waveform h by summing over modes, using SpinWeightedSpheroidalHarmonic and exponential factors. Includes both positive and negative m contributions.

```python
h = 0.j
for mode in modes:
    j, m, k, n, omega, amp = mode
    Slm = SpinWeightedSpheroidalHarmonic(-2, j, m, a*omega)(th)
    h += amp * Slm * np.exp(1j * (m * phi - omega * u_grid))
    Slm = SpinWeightedSpheroidalHarmonic(-2, j, -m, -a*omega)(th)
    h += (-1)**(j + k) * np.conj(amp) * Slm * np.exp(-1j * (m * phi - omega * u_grid))
```

--------------------------------

### Compute total fluxes by summing over j and m modes

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/fluxes.ipynb

Loops over spheroidal mode numbers j from 2 to jmax and azimuthal mode numbers m from 1 to j, solving each TeukolskyMode with auto_solve=True, computing the FluxMode, and accumulating the total, horizon, and infinity fluxes. The final fluxes are doubled to account for both positive and negative m modes, and the per-j contributions are stored.

```python
jmax = 18
total_fluxes = np.zeros(3)
horizon_fluxes = np.zeros(3)
infinity_fluxes = np.zeros(3)
total_fluxes_jmodes = np.zeros((jmax - 1, 3))
for j in range(2, jmax+1):
    for m in range(1, j+1):
        Psis = TeukolskyMode(-2, j, m, 0, 0, geo, auto_solve=True)
        flux = FluxMode(geo, Psis)
        total_fluxes += flux.totalfluxes
        total_fluxes_jmodes[j-2] += flux.totalfluxes
        horizon_fluxes += flux.horizonfluxes
        infinity_fluxes += flux.infinityfluxes
total_fluxes *= 2
horizon_fluxes *= 2
infinity_fluxes *= 2
total_fluxes_jmodes *= 2
```

--------------------------------

### Check orbital parameters properties against geo.orbitalparameters

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/geodesics.ipynb

Verify that the individual orbital parameter properties match the tuple returned by geo.orbitalparameters. The output is an array of booleans, all True.

```python
(geo.blackholespin, geo.semilatusrectum, geo.eccentricity, geo.inclination) == geo.orbitalparameters
```

--------------------------------

### Regenerate performance figures

Source: https://github.com/znasipak/pybhpt/blob/main/docs/performance/index.md

Generates the performance figures from the committed CSV data and writes them to docs/_static/figures/. Run this after updating the CSVs to refresh the documentation figures.

```bash
python -m benchmarks.plots            # figures -> docs/_static/figures/
```

--------------------------------

### Instantiate TeukolskyMode with mode numbers

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/fluxes.ipynb

Creates a TeukolskyMode object with spin-weight s=-2, spheroidal polar mode number j=12, azimuthal mode number m=3, polar mode number k=1, and radial mode number n=2, using the previously defined geodesic.

```python
s, j, m, k, n = -2, 12, 3, 1, 2
Psis = TeukolskyMode(s, j, m, k, n, geo)
```

--------------------------------

### Compute geodesic evolution vs coordinate time

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/geodesics.ipynb

Compute the geodesic evolution as a function of coordinate time t. First solve for Mino time lambda(t) using geo.mino_of_t, then evaluate the geodesic at those lambda values.

```python
t_grid = np.linspace(0, 1500, 1000)
la_t_grid = geo.mino_of_t(t_grid)
geo_t_grid = geo(la_t_grid)
```

--------------------------------

### Cython code generation custom command

Source: https://github.com/znasipak/pybhpt/blob/main/CMakeLists.txt

Generates redshift_wrap.cpp from Cython .pyx files using the Python executable. Depends on multiple .pyx files and outputs to the build directory.

```cmake
file(MAKE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/cython")

add_custom_command(
  OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/cython/redshift_wrap.cpp"
  COMMAND ${Python_EXECUTABLE} -m cython
          "${CMAKE_CURRENT_SOURCE_DIR}/cython/redshift_wrap.pyx"
          --output-file "${CMAKE_CURRENT_BINARY_DIR}/cython/redshift_wrap.cpp"
          -3 -+ --module-name "cybhpt_full"
          -I "${CMAKE_CURRENT_SOURCE_DIR}/cpp/include"
  DEPENDS "${CMAKE_CURRENT_SOURCE_DIR}/cython/redshift_wrap.pyx"
          "${CMAKE_CURRENT_SOURCE_DIR}/cython/flux_wrap.pyx"
          "${CMAKE_CURRENT_SOURCE_DIR}/cython/teukolsky_wrap.pyx"
          "${CMAKE_CURRENT_SOURCE_DIR}/cython/geo_wrap.pyx"
          "${CMAKE_CURRENT_SOURCE_DIR}/cython/radialsolver_wrap.pyx"
          "${CMAKE_CURRENT_SOURCE_DIR}/cython/swsh_wrap.pyx"
  COMMENT "Generating redshift_wrap.cpp from Cython"
  VERBATIM
)
```

--------------------------------

### Print total fluxes

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/fluxes.ipynb

Prints the total energy, angular momentum, and Carter constant fluxes after summing over all modes.

```python
print("Total fluxes (E, Lz, Q): ", total_fluxes)
```

--------------------------------

### Optional LAPACK detection and linking

Source: https://github.com/znasipak/pybhpt/blob/main/CMakeLists.txt

Optionally finds LAPACK. If found, links it and defines PYBHPT_HAS_LAPACK; otherwise falls back to slow gsl_eigen_symm. Set PYBHPT_REQUIRE_LAPACK=ON to fail the build if LAPACK is missing.

```cmake
# LAPACK — used in swsh.cpp for fast symmetric eigenvalue computation.
# gsl_eigen_symm is ~100x slower than LAPACK dsyev for small matrices.
# By default LAPACK is optional (graceful fallback). Set -DPYBHPT_REQUIRE_LAPACK=ON
# (done in the wheel CI and the Apptainer image) to fail the build instead of
# silently degrading to the slow path. Hint the implementation with -DBLA_VENDOR=...
# (e.g. OpenBLAS) and point CMAKE_PREFIX_PATH at the conda env if needed.
option(PYBHPT_REQUIRE_LAPACK "Fail the build if LAPACK is not found" OFF)
if(PYBHPT_REQUIRE_LAPACK)
  find_package(LAPACK REQUIRED)
else()
  find_package(LAPACK)
endif()
if(LAPACK_FOUND)
  target_link_libraries(cybhpt_full PRIVATE ${LAPACK_LIBRARIES})
  target_compile_definitions(cybhpt_full PRIVATE PYBHPT_HAS_LAPACK)
  message(STATUS "LAPACK found: ${LAPACK_LIBRARIES}")
else()
  message(STATUS "LAPACK not found — spectral solver will use (slow) gsl_eigen_symm")
endif()
```

--------------------------------

### Import NumPy

Source: https://github.com/znasipak/pybhpt/blob/main/notebooks/flux-example.ipynb

Imports the NumPy library for numerical operations.

```python
import numpy as np
```

--------------------------------

### Plot radial derivatives for In and Up

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/teuk.ipynb

Plots the first and second radial derivatives of the 'In' and 'Up' solutions on a log scale. The output is a figure with two subplots.

```python
dRin = np.array([Psis.radialderivative('In', i) for i in range(len(r_grid))])
dRup = np.array([Psis.radialderivative('Up', i) for i in range(len(r_grid))])
d2Rin = np.array([Psis.radialderivative2('In', i) for i in range(len(r_grid))])
d2Rup = np.array([Psis.radialderivative2('Up', i) for i in range(len(r_grid))])

fig, axs = plt.subplots(2, 1, figsize=(8, 12), sharex=True)
axs[0].plot(r_grid, np.abs(dRin.real), label='Re$[\partial_r R]$')
axs[0].plot(r_grid, np.abs(dRin.imag), label='Im$[\partial_r R]$')
axs[0].plot(r_grid, np.abs(d2Rin.real), '--', label='Re$[\partial_r^2 R]$')
axs[0].plot(r_grid, np.abs(d2Rin.imag), '--', label='Im$[\partial_r^2 R]$')
axs[0].set_ylabel(r'In solutions')
axs[0].set_yscale('log')
axs[0].legend()

axs[1].plot(r_grid, np.abs(dRup.real), label='Re$[\partial_r^2 R]$')
axs[1].plot(r_grid, np.abs(dRup.imag), label='Im$[\partial_r^2 R]$')
axs[1].plot(r_grid, np.abs(d2Rup.real), '--', label='Re$[\partial_r^2 R]$')
axs[1].plot(r_grid, np.abs(d2Rup.imag), '--', label='Im$[\partial_r^2 R]$')
axs[1].set_ylabel(r'Up solutions')
axs[1].set_xlabel(r'$r$')
axs[1].set_yscale('log')

axs[1].legend()
plt.tight_layout()
plt.show()
```

```text
Result:
<Figure size 800x1200 with 2 Axes>
```

--------------------------------

### Define convergence test and FluxSummation class

Source: https://github.com/znasipak/pybhpt/blob/main/notebooks/flux-example.ipynb

Defines a convergence test function and a FluxSummation class that sums over frequency modes to compute total, horizon, and infinity fluxes for a particle on a Kerr geodesic. The class initializes with spin, semi-latus rectum, eccentricity, and inclination, and tracks peak modes.

```python
def convergenceTestModes(fluxMode, fluxTot, fluxComp, fluxComp2, tol = 1e-8, flux_position=[0, 1, 2]):
    if np.all(np.abs(np.array(fluxMode))[flux_position]/np.array(fluxTot)[flux_position] < tol) and np.all(np.abs(np.array(fluxComp))[flux_position]/np.array(fluxTot)[flux_position] < tol) and np.all(np.abs(np.array(fluxMode)[flux_position]) < np.abs(np.array(fluxComp)[flux_position])) and np.all(np.abs(np.array(fluxMode)[flux_position]) < np.abs(np.array(fluxComp2)[flux_position])):
        return True
    else:
        return False

class FluxSummation:
    '''
    Class for summing over all of the frequency modes that contribute to the total gravitational wave (GW()) fluxes
    produced by a particle on a geodesic orbit parametrized by the parameters $(a, p, e, x)$

    :param a: Kerr spin parameter
    :param p: Semi-latus rectum
    :param e: Orbital eccentricity
    :param x: Projection of the orbital inclination
    :param nsamples0: Initial number of sampling points to use in orbital integration of fluxes

    :attribute totalfluxes: Three dimensional array of GW fluxes, [dE/dt, dL/dt, dQ/dt]
    :attribute horizonfluxes: Three dimensional array of GW fluxes radiated down the black hole horizon
    :attribute infinityfluxes: Three dimensional array of GW fluxes radiated away at infinity
    '''
    def __init__(self, a, p, e, x, nsamples0 = None, tol = 1.e-8):
        self.tol = tol
        self.a = a
        self.p = p
        self.e = e
        self.x = x
        if nsamples0 is None:
            nsamples0 = 2**7
        self.nsamples0 = nsamples0
        self.nsamples = nsamples0
        self.geo_dict = {self.nsamples: KerrGeodesic(self.a, self.p, self.e, self.x, self.nsamples) }
        self.geo = self.geo_dict[self.nsamples]
        self.nlimit = 300
        self.klimit = 100
        self.mlimit = 100
        self.llimit = 100
        self.fluxref = 0.5*np.array([32./5.*np.abs(self.geo.frequencies[-1])**(10./3.), 32./5.*np.abs(self.geo.frequencies[-1])**(7./3.), 64./5.*np.abs(self.geo.frequencies[-1])**(2.)])
        print(f"Flux reference: {self.fluxref}")

        # self.nfluxes_full = np.zeros((2*self.nlimit, 3))
        # self.kfluxes_full = np.zeros((2*self.klimit, 3))
        self.totalfluxes = np.zeros(3)
        self.horizonfluxes = np.zeros(3)
        self.infinityfluxes = np.zeros(3)
        self.nfluxes = np.zeros(3)
        self.kfluxes = np.zeros(3)
        self.mfluxes = np.zeros(3)
        self.nmax = 0
        self.nmin = 0
        self.kmax = 0
        self.kmin = 0
        self.kmaxmodenums = np.zeros((2*self.nlimit + 1, 2), dtype=int)
        self.mmaxmodenums = np.zeros((2*self.nlimit + 1, 2*self.klimit + 1, 2), dtype=int)
        self.lmaxmodenums = np.zeros((2*self.nlimit + 1, 2*self.klimit + 1, self.mlimit + 1, 2), dtype=int)
        
        self.kpeakmodenums = np.zeros((2*self.nlimit + 1), dtype=int)
        self.kpeakmodes = np.zeros((2*self.llimit + 1, 3))
        self.mpeakmodenums = np.zeros((2*self.klimit + 1), dtype=int)
        self.mpeakmodes = np.zeros((2*self.llimit + 1, 3))
        self.fluxamplitudes = np.zeros((self.llimit + 1, self.mlimit + 1, 2*self.klimit + 1, 2*self.nlimit + 1), dtype=np.complex128)

        self.n = 0
        self.k = 0
        self.m = 0
        self.l = 2
        self.flux_positions = [0, 1, 2]

    def check_for_peak_m_mode(self):
        current_peak_value = self.mpeakmodes[self.k]
        if np.all(np.abs(self.mfluxes) > current_peak_value):
            current_peak_value = np.abs(self.mfluxes)
            self.mpeakmodes[self.k] = current_peak_value
            self.mpeakmodenums[self.k] = self.m
            # print(f"New m-peak is {self.m} for {(self.k, self.n)} with flux values {self.mfluxes}")

    def check_for_peak_k_mode(self):
        current_peak_value = self.kpeakmodes[self.n]
        if np.all(np.abs(self.kfluxes) > current_peak_value):
            current_peak_value = np.abs(self.kfluxes)
            self.kpeakmodes[self.n] = current_peak_value
            self.kpeakmodenums[self.n] = self.k
            # print(f"New k-peak is {self.k} for {(self.n)} with flux values {self.kfluxes}")
```

--------------------------------

### Set up circular Kerr geodesic

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/fluxes.ipynb

Defines a circular geodesic with eccentricity e=0 and cos(inclination) x=1, using a small number of samples (nsamples=4) for the Fourier decomposition.

```python
a, p, e, x = 0.99, 5, 0., 1.
geo = KerrGeodesic(a, p, e, x, nsamples=2**2)
```

--------------------------------

### Import pybhpt.radial and plotting libraries

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/radial.ipynb

Import the RadialTeukolsky class from pybhpt along with numpy and matplotlib for array operations and plotting.

```python
from pybhpt.radial import RadialTeukolsky
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

```

--------------------------------

### Compute metric coefficients for IRG and ORG gauges

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/tutorial.ipynb

Creates MetricCoefficients objects for the IRG and ORG gauges at theta = pi/2, using the radial grid and geodesic parameters.

```python
from pybhpt.metric import MetricCoefficients

th = np.array([0.5*np.pi])
habIRG = MetricCoefficients("IRG", a, r, th)
habORG = MetricCoefficients("ORG", a, r, th)
```

--------------------------------

### Define a background Kerr geodesic

Source: https://github.com/znasipak/pybhpt/blob/main/docs/notebooks/waveform.ipynb

Creates a KerrGeodesic object with parameters a=0.99, p=5, e=0.6, x=0.4, and 2^8 samples.

```python
a, p, e, x = 0.99, 5, 0.6, 0.4
geo = KerrGeodesic(a, p, e, x, nsamples = 2**8)
```