pySEOBNR

https://waveforms.docs.ligo.org/software/pyseobnr
Admin
pySEOBNR is a Python package for generating effective-one-body gravitational waveform models for...
Tokens:259,850
Snippets:1,088
Trust Score:10
Update:1 month ago
Show doc for...
Tokens:
Raw
### Install documentation dependencies

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/installation.rst.txt

Install pandoc and the required documentation build packages.

```console
(.conda_venv) $ conda install --yes --channel conda-forge pandoc
```

```console
(.venv) $ apt-get update && apt-get install -y pandoc
```

```console
(.venv) $ brew install pandoc
```

```console
(.venv) $ pip install ".[docs]"
(.venv) $ cd docs
(.venv) $ make html
```

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

### Install and configure pre-commit

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/how_to_contribute.rst.txt

Commands to install the pre-commit package and initialize hooks within the repository.

```console
# install the pre-commit package
pip install pre-commit
cd pyseobnr
pre-commit install
```

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

### Install pyseobnr from source

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/installation.rst.txt

Clone the repository and install the package using pip in a clean environment.

```console
(.venv) $ git clone https://git.ligo.org/waveforms/software/pyseobnr.git
(.venv) $ cd pyseobnr
(.venv) $ pip install -U pip wheel setuptools numpy
(.venv) $ pip install .
```

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

### Importing necessary libraries

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/notebooks/example_precession.md.txt

Initial setup for plotting and numerical operations.

```python
import warnings

import matplotlib.pyplot as plt
import numpy as np
```

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

### Generate documentation locally

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/how_to_contribute.rst.txt

Use this command to build the project documentation and automatically install necessary dependencies.

```console
tox -e docs
```

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

### bbh_UIBfits_setup

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/utils/nr_utils.html

A common setup function that prepares input parameters for UIB final-state and luminosity fit functions.

```APIDOC
## bbh_UIBfits_setup

### Description
Prepares and validates input parameters (masses and spins) for UIB fit functions. It calculates derived quantities such as symmetric mass ratio, total spin, and effective spin.

### Parameters
- **m1** (float/array) - Component mass 1
- **m2** (float/array) - Component mass 2
- **chi1** (float/array) - Dimensionless spin of BH 1
- **chi2** (float/array) - Dimensionless spin of BH 2
```

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

### Install optional dependencies

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/installation.rst.txt

Install extra packages for running sanity checks or external waveform tools.

```console
(.venv) $ pip install ".[checks]"
```

```console
(.venv) $ pip install git+https://bitbucket.org/sergei_ossokine/waveform_tools
```

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

### Setup UIB fit parameters

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/utils/nr_utils.html

Prepares common variables for UIB final-state and luminosity fit functions, including mass ratios, symmetric mass ratios, and spin combinations.

```python
def bbh_UIBfits_setup(m1, m2, chi1, chi2):
    """
    common setup function for UIB final-state and luminosity fit functions
    """

    # Vectorize the function if arrays are provided as input
    m1 = np.vectorize(float)(np.array(m1))
    m2 = np.vectorize(float)(np.array(m2))
    chi1 = np.vectorize(float)(np.array(chi1))
    chi2 = np.vectorize(float)(np.array(chi2))

    if np.any(m1 < 0):
        raise ValueError("m1 must not be negative")
    if np.any(m2 < 0):
        raise ValueError("m2 must not be negative")

    if np.any(abs(chi1) > 1):
        raise ValueError("chi1 has to be in [-1, 1]")
    if np.any(abs(chi2) > 1):
        raise ValueError("chi2 has to be in [-1, 1]")

    # binary masses
    m = m1 + m2
    if np.any(m <= 0):
        raise ValueError("m1+m2 must be positive")
    msq = m * m
    m1sq = m1 * m1
    m2sq = m2 * m2

    # symmetric mass ratio
    eta = m1 * m2 / msq
    if np.any(eta > 0.25):
        print(
            "Truncating eta from above to 0.25. This should only be necessary in some rounding corner cases, "
            "but better check your m1 and m2 inputs..."
        )
        eta = np.minimum(eta, 0.25)
    if np.any(eta < 0.0):
        print(
            "Truncating negative eta to 0.0. This should only be necessary in some rounding corner cases, "
            "but better check your m1 and m2 inputs..."
        )
        eta = np.maximum(eta, 0.0)
    eta2 = eta * eta
    eta3 = eta2 * eta
    eta4 = eta2 * eta2

    # spin variables (in m = 1 units)
    S1 = chi1 * m1sq / msq  # spin angular momentum 1
    S2 = chi2 * m2sq / msq  # spin angular momentum 2
    Stot = S1 + S2  # total spin
    Shat = (chi1 * m1sq + chi2 * m2sq) / (
        m1sq + m2sq
    )  # effective spin, = msq*Stot/(m1sq+m2sq)
    Shat2 = Shat * Shat
    Shat3 = Shat2 * Shat
    Shat4 = Shat2 * Shat2

    # asymmetric spin combination (spin difference), where the paper assumes m1>m2
    # to make our implementation symmetric under simultaneous exchange of m1<->m2 and chi1<->chi2,
    # we flip the sign here when m2>m1
    chidiff = chi1 - chi2
    if np.any(m2 > m1):
        chidiff = np.sign(m1 - m2) * chidiff
    chidiff2 = chidiff * chidiff

    # typical squareroots and functions of eta
    sqrt2 = 2.0**0.5
    sqrt3 = 3.0**0.5
    sqrt1m4eta = (1.0 - 4.0 * eta) ** 0.5

    return (
        m,
        eta,
        eta2,
        eta3,
        eta4,
        Stot,
        Shat,
        Shat2,
        Shat3,
        Shat4,
        chidiff,
        chidiff2,
        sqrt2,
        sqrt3,
        sqrt1m4eta,
    )
```

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

### Modifying minimum starting separation

Source: https://waveforms.docs.ligo.org/software/pyseobnr/source/notebooks/example_SEOBNRv5EHM.html

Configures the minimum starting separation to 7M, preventing the activation of secular backward integration.

```python
settings = dict(
    secular_bwd_int=True,
    r_start_min=7.0,
)

t, modes = generate_modes_opt(
    q,
    chi_1,
    chi_2,
    omega_start,
    eccentricity=0.5,
    rel_anomaly=0.0,
    approximant="SEOBNRv5EHM",
    settings=settings,
)
```

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

### compute_starting_values(r_stop, dynamics, threshold, params)

Source: https://waveforms.docs.ligo.org/software/pyseobnr/api/pyseobnr.eob.utils.utils_eccentric.compute_starting_values.html

Computes the exact values of starting time, eccentricity, relativistic anomaly, and dimensionless starting orbit-averaged orbital frequency that correspond to a specific separation r_stop.

```APIDOC
## compute_starting_values(r_stop, dynamics, threshold, params)

### Description
Computes the exact values of starting time, eccentricity, relativistic anomaly, and dimensionless starting orbit-averaged orbital frequency that correspond to r = r_stop.

### Parameters
- **r_stop** (float) - Required - Value of the separation at which to compute the starting values.
- **dynamics** (ndarray) - Required - Dynamical variables (t, eccentricity, rel_anomaly, omega_avg^{2/3}).
- **threshold** (bool) - Required - If False, then perform a root solving procedure to find the time at which r = r_stop.
- **params** (EOBParams) - Required - EOB parameters of the system.

### Returns
- **tuple** - (t_start, eccentricity, rel_anomaly, omega_avg) representing the starting values at r = r_stop.
```

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

### Setup BBH HBR2016 Calculations

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/utils/nr_utils.html

Prepares mass and spin inputs for HBR final spin fits by vectorizing inputs and calculating the mass ratio.

```python
def _bbh_HBR2016_setup(m1, m2, chi1, chi2):
    """
    Setup function for the Hofmann, Barausse, and Rezzolla final spin fits to vectorize the masses
    and spins and calculate the mass ratio.
    """

    # Vectorize if arrays are provided as input
    m1 = np.vectorize(float)(np.array(m1))
    m2 = np.vectorize(float)(np.array(m2))
    chi1 = np.vectorize(float)(np.array(chi1))
    chi2 = np.vectorize(float)(np.array(chi2))

    return m1, m2, chi1, chi2, m2 / m1
```

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

### Disabling secular backwards integration

Source: https://waveforms.docs.ligo.org/software/pyseobnr/source/notebooks/example_SEOBNRv5EHM.html

Demonstrates how to deactivate secular backwards integration by setting the 'secular_bwd_int' flag to False, which will cause an error if the starting separation is below 10M.

```python
try:
    settings = dict(
        secular_bwd_int=False,
    )

    t, modes = generate_modes_opt(
        q,
        chi_1,
        chi_2,
        omega_start,
        eccentricity=0.5,
        rel_anomaly=0.0,
        approximant="SEOBNRv5EHM",
        settings=settings,
    )
except Exception as e:
    print(e)
```

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

### bbh_UIBfits_setup

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/api/pyseobnr.eob.utils.nr_utils.rst.txt

Initializes the UIB fits for binary black hole systems.

```APIDOC
## bbh_UIBfits_setup

### Description
Initializes the UIB fits required for binary black hole waveform modeling.

### Signature
`bbh_UIBfits_setup()`
```

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

### Compute Starting Values for EOB Dynamics

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/utils/utils_eccentric.html

Determines the starting time, eccentricity, relativistic anomaly, and frequency at a specific radial separation.

```python
def compute_starting_values(
    r_stop: float,
    dynamics: np.ndarray,
    threshold: bool,
    params: EOBParams,
) -> tuple[float, float, float, float]:
    """
    Compute the exact values of starting time, eccentricity,
    relativistic anomaly, and dimensionless starting orbit-averaged
    orbital frequency that correspond to r = r_stop.

    Args:
        r_stop: Value of the separation at which to compute the
            starting values
        dynamics: Dynamical variables
            (t, eccentricity, rel_anomaly, omega_avg^{2/3})
        threshold: If False, then perform a root solving procedure
            to find the time at which r = r_stop
        params: EOB parameters of the system

    Returns:
        tuple: (t_start, eccentricity, rel_anomaly, omega_avg).
            Starting values at r = r_stop.
    """

    if threshold:
        t_start, eccentricity, rel_anomaly, omega_avg = (
            dynamics[0, 0],
            dynamics[0, 1],
            float(dynamics[0, 2]) % (2 * np.pi),
            dynamics[0, 3] ** 1.5,
        )
```

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

### SEOBNRv5RRForce.__init__

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/api/pyseobnr.eob.waveform.waveform.SEOBNRv5RRForce.rst.txt

Initializes the SEOBNRv5RRForce class instance.

```APIDOC
## __init__(self)

### Description
Initializes a new instance of the SEOBNRv5RRForce class.

### Usage
```python
from pyseobnr.eob.waveform.waveform import SEOBNRv5RRForce
force = SEOBNRv5RRForce()
```
```

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

### compute_starting_values

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/api/pyseobnr.eob.utils.utils_eccentric.compute_starting_values.rst.txt

Calculates the starting values for eccentric EOB waveform generation.

```APIDOC
## compute_starting_values

### Description
Calculates the starting values required for the evolution of eccentric EOB waveforms.

### Signature
`pyseobnr.eob.utils.utils_eccentric.compute_starting_values(...)`
```

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

### compute_quasiprecessing_PNdynamics_opt

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/dynamics/pn_evolution_opt.html

Computes the dynamics starting from omega_start, with spins defined at omega_ref.

```APIDOC
## compute_quasiprecessing_PNdynamics_opt(omega_ref, omega_start, m_1, m_2, chi_1, chi_2, rtol=1e-10, atol=1e-12, backend='DOP853')

### Description
Compute the dynamics starting from omega_start, with spins defined at omega_ref. PN evolution equations are integrated to get spin and orbital angular momentum, then aligned-spin EOB dynamics are integrated.

### Parameters
- **omega_ref** (float) - Reference frequency
- **omega_start** (float) - Starting frequency
- **m_1** (float) - Mass of primary
- **m_2** (float) - Mass of secondary
- **chi_1** (np.ndarray) - Dimensionless spin of the primary
- **chi_2** (np.ndarray) - Dimensionless spin of the secondary
- **rtol** (float, optional) - Relative tolerance for EOB integration. Defaults to 1e-12.
- **atol** (float, optional) - Absolute tolerance for EOB integration. Defaults to 1e-12.
- **backend** (str, optional) - The backend to use for ODE integration. Defaults to 'DOP853'.

### Returns
- **tuple** - Aligned-spin EOB dynamics, PN time, PN dynamics, PN splines
```

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

### Initialize SEOBNRv5EHM Parameters

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/models/SEOBNRv5EHM.html

Sets up eccentricity, radiation reaction, and integration flags during model initialization.

```python
# Parameters employed in the eccentric model
        self.eccentricity = eccentricity
        self.rel_anomaly = rel_anomaly
        self.omega_start = omega_start
        self.f_start = omega_start / (self.M * lal.MTSUN_SI * np.pi)
        self._radiation_reaction: RadiationReactionForceEcc = RR
        self._radiation_reaction_is_initialized = False
        self.EccIC = self.settings["EccIC"]
        self.flags_ecc = (
            dict(
                flagPN12=1,
                flagPN1=1,
                flagPN32=1,
                flagPN2=1,
                flagPN52=1,
                flagPN3=1,
                flagPA=1,
                flagPA_modes=1,
                flagTail=1,
                flagMemory=1,
            )
            | self.settings["flags_ecc"]
        )
        self.atol = self.settings["atol"]
        self.dissipative_ICs = self.settings["dissipative_ICs"]
        self.e_stop = self.settings["e_stop"]
        self.h_0 = self.settings["h_0"]
        self.IC_messages = self.settings["IC_messages"]
        self.inspiral_modes = self.settings["inspiral_modes"]
        self.nqc_method = self.settings["nqc_method"]
        self.r_min = self.settings["r_min"]
        self.r_start_min = self.settings["r_start_min"]
        self.rtol = self.settings["rtol"]
        self.secular_bwd_int = self.settings["secular_bwd_int"]
        self.t_backwards = self.settings["t_backwards"]
        self.t_bwd_secular = self.settings["t_bwd_secular"]
        self.t_max = self.settings["t_max"]
        self.warning_bwd_int = self.settings["warning_bwd_int"]
        self.warning_secular_bwd_int = self.settings["warning_secular_bwd_int"]
        self.y_init = self.settings["y_init"]
        assert self.settings["integrator"] in get_args(IntegratorType)
        self.integrator: IntegratorType = cast(
            IntegratorType, self.settings["integrator"]
        )
```

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

### Configuring parameters and initializing the generator

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/notebooks/example_SEOBNRv5EHM.md.txt

Sets up physical and extrinsic parameters in a dictionary and retrieves the waveform generator for SEOBNRv5EHM.

```python
# start with the usual parameter definitions

# Masses
m1 = 30.0
m2 = 30.0

# spings
s1x = 0.0
s1y = 0.0
s1z = 0.5
s2x = 0.0
s2y = 0.0
s2z = -0.2

# Extrinsic parameters
deltaT = 1.0 / 4096.0 / 2.0
f_min = 10.0
f_ref = 10.0
luminosity_distance = 1000.0

iota = np.pi / 3.0
phase = 0.0

# Eccentric parameters
eccentricity = 0.2
longitude_ascending_nodes = 0.0
meanPerAno = 0.5

# Conditioning
condition = 1
deltaF = 1.0 / 8.0  # /20.
# Create dict for gwsignal generator
gwsignal_dict = {
    "mass1": m1 * u.solMass,
    "mass2": m2 * u.solMass,
    "spin1x": s1x * u.dimensionless_unscaled,
    "spin1y": s1y * u.dimensionless_unscaled,
    "spin1z": s1z * u.dimensionless_unscaled,
    "spin2x": s2x * u.dimensionless_unscaled,
    "spin2y": s2y * u.dimensionless_unscaled,
    "spin2z": s2z * u.dimensionless_unscaled,
    #'deltaF' : delta_frequency * u.Hz,
    "deltaT": deltaT * u.s,
    #'deltaF': deltaF*u.Hz,
    "f22_start": f_min * u.Hz,
    "f_max": f_max * u.Hz,
    "f22_ref": f_ref * u.Hz,
    "phi_ref": phase * u.rad,
    "distance": luminosity_distance * u.Mpc,
    "inclination": iota * u.rad,
    "eccentricity": eccentricity * u.dimensionless_unscaled,
    "longAscNodes": longitude_ascending_nodes * u.rad,
    "meanPerAno": meanPerAno * u.rad,
    # 'ModeArray': mode_array,
    "condition": 0,
    "lmax": 4,
    "lmax_nyquist": 4,
    # Some additional specific SEOBNRv5EHM parameters
    "secular_bwd_int": True,
    "warning_secular_bwd_int": False,
    "t_backwards": 0,
    "warning_bwd_int": False,
}

waveform_approximant = "SEOBNRv5EHM"

try:
    wf_gen = gwsignal_get_waveform_generator(waveform_approximant)
except ValueError as e:
    if str(e) != "Approximant not implemented in GWSignal!":
        raise

    wf_gen = None
    print("SEOBNRv5EHM not supported by this version of lal")
```

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

### Configure GSL for pygsl-lite

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/installation.rst.txt

Force the installation of pygsl-lite to link against a specific system GSL library.

```console
(.venv) $ GSL_HOME=/path/to/your/GSL/HOME pip install --force-reinstall pygsl-lite
```

```console
(.venv) $ GSL_HOME=$(brew --prefix gsl) pip install --force-reinstall pygsl-lite
```

```console
(.venv) $ pip install \
    --use-pep517 \
    --config-setting="--global-option=build_ext" \
    --config-setting="--build-option=-I$(brew --prefix gsl)/include/" \
    --config-setting="--build-option=-L$(brew --prefix gsl)/lib/" \
    --force-reinstall \
    --no-binary pygsl_lite \
    --no-cache-dir \
    pygsl_lite
```

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

### Conditioned Time-Domain Polarizations (Tapered)

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/generate_waveform.html

Generates time-domain polarizations with tapering applied at the start of the waveform.

```python
    def generate_td_polarizations_conditioned_1(self):
        """
        Generate time-domain polarizations, with tappering at the beginning
        of the waveform, returned as LAL REAL8TimeSeries.
        """

        hp_lal, hc_lal = self.generate_td_polarizations()
        lalsim.SimInspiralREAL8WaveTaper(hp_lal.data, 1)
        lalsim.SimInspiralREAL8WaveTaper(hc_lal.data, 1)

        return hp_lal, hc_lal
```

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

### Importing PySEOBNR and dependencies

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/notebooks/getting_started.md.txt

Initializes the environment by silencing LAL warnings and importing necessary libraries including GenerateWaveform.

```python
import warnings

# silence warnings coming from LAL
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")

import numpy as np
from matplotlib import pyplot as plt

# import the library
from pyseobnr.generate_waveform import GenerateWaveform
```

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

### Initialize Environment and Imports

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/notebooks/pseob_example.md.txt

Imports necessary libraries and suppresses LAL warnings.

```python
import warnings
from copy import deepcopy

import numpy as np
from matplotlib import pyplot as plt

# silence warnings coming from LAL
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
```

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

### Initialize SEOBNRv5HM model parameters

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/models/SEOBNRv5HM.html

Configures mode selection, Hamiltonian initialization, and validates post-adiabatic and initial condition settings.

```python
# All the modes we will need to output
        self.return_modes = self.settings.get("return_modes", None)

        # Check that the modes are valid, i.e. something we
        # can return
        self.max_ell_returned = self._validate_modes(user_settings=settings)
        self.lmax_nyquist = self.settings.get("lmax_nyquist", self.max_ell_returned)
        # Now deal with which mixed modes the user wants, if any

        self.mixed_modes = [(3, 2), (4, 3)]
        self.mixed_modes = [x for x in self.mixed_modes if x in self.return_modes]

        # All the modes we need to compute. This can be a larger list
        # than the returned modes, e.g. when we need certain modes to
        # do mode mixing
        self.computed_modes = deepcopy(self.return_modes)
        # Make sure the array contains what we need
        self._ensure_consistency()
        self._validate_antisymmetric_parameters()
        self._validate_convention_parameters()
        self._initialize_params(phys_pars=self.phys_pars)

        # Initialize the Hamiltonian
        self.H = H(self.eob_pars)

        # Uncomment and comment the line above to make the python Hamiltonian work
        # self.H = H()

        if self.settings["postadiabatic_type"] not in ["numeric", "analytic"]:
            raise ValueError("Incorrect value for postadiabatic_type")

        if self.settings["initial_conditions"] not in get_args(
            InitialConditionTypesPrecessing
        ):
            raise ValueError("Incorrect value for initial_conditions")

        if self.settings["initial_conditions_postadiabatic_type"] not in get_args(
            InitialConditionPostadiabaticTypesPrecessing
        ):
            raise ValueError(
                "Incorrect value for initial_conditions_postadiabatic_type"
            )

        # sign of the final spin
        self._sign_final_spin: int | None = None

        # Shift for setting phase of (2,2)-coprecessing mode to 0 at t_ref
        # applied to the polarizations
        self.alpha_diff = 0.0  # Initialize to 0 - default it has no effect
```

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

### Compute starting orbital values

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/models/SEOBNRv5EHM.html

Calculates initial orbital frequencies and separation based on the chosen EccIC method.

```python
def _compute_starting_values(self):
    """
    Compute the starting values of the orbit-averaged orbital frequency,
    the instantaneous orbital frequency, and the starting separation of
    the system.
    """

    flags = self.flags_ecc.copy()
    flags["flagPA"] = 0  # This is employed in conservative equations
    if self.EccIC == 0:
        self.omega_inst = self.omega_start
        self.RR.evolution_equations.compute(
            e=self.eccentricity,
            omega=self.omega_start,
            z=self.rel_anomaly,
        )
        self.x_avg = self.RR.evolution_equations.get("xavg_omegainst")
        self.omega_avg = self.x_avg**1.5
    elif self.EccIC == 1:
        self.omega_avg = self.omega_start
        self.x_avg = self.omega_avg ** (2.0 / 3.0)
        self.omega_inst = dot_phi_omega_avg_e_z(
            self.omega_avg,
            self.eccentricity,
            self.rel_anomaly,
            self.nu,
            self.m_1 - self.m_2,
            self.chi_A,
            self.chi_S,
            flags,
        )
    else:
        raise NotImplementedError(
            "Select a supported value of 'EccIC': "
            "0 (for starting instantaneous orbital frequency) "
            "or 1 (for starting orbit-averaged orbital frequency)."
        )
    flags["flagPA"] = 1
    self.r_start = r_omega_avg_e_z(
        self.omega_avg,
        self.eccentricity,
        self.rel_anomaly,
        self.nu,
        self.m_1 - self.m_2,
        self.chi_A,
        self.chi_S,
        flags,
    )
    self.eob_pars.ecc_params.omega_avg = self.omega_avg
    self.eob_pars.ecc_params.x_avg = self.x_avg
    self.eob_pars.ecc_params.omega_inst = self.omega_inst
    self.eccentricity_ref = self.eccentricity
    self.rel_anomaly_ref = self.rel_anomaly
    self.omega_avg_ref = self.omega_avg
```

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

### SEOBNRv5HM_opt.__init__

Source: https://waveforms.docs.ligo.org/software/pyseobnr/api/pyseobnr.models.SEOBNRv5HM.SEOBNRv5HM_opt.html

Initializes the SEOBNRv5HM_opt model with mass ratio, spin components, initial frequency, Hamiltonian, and radiation reaction force.

```APIDOC
## class SEOBNRv5HM_opt

### Description
Initializes the SEOBNRv5 approximant model.

### Parameters
- **q** (float) - Required - Mass ratio m1/m2 >= 1
- **chi_1** (float) - Required - z component of the dimensionless spin of primary
- **chi_2** (float) - Required - z component of the dimensionless spin of secondary
- **omega0** (float) - Required - Initial orbital frequency, in geometric units
- **H** (Hamiltonian) - Required - Hamiltonian class
- **RR** (Callable) - Required - RR force
- **settings** (Dict[Any, Any]) - Optional - The settings. Defaults to None.
```

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

### compute_dynamics_ecc_secular_opt

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/dynamics/integrate_ode_ecc.html

Integrates the secular eccentric dynamics backwards in time to find starting conditions corresponding to a specific separation.

```APIDOC
## compute_dynamics_ecc_secular_opt

### Description
Main function to integrate the secular eccentric dynamics backwards in time. The RHS of the employed ODE system is given by Eqs. (18a)-(18c) in [Gamboa2024a].

### Parameters
- **eccentricity** (float) - Required - Eccentricity of the system
- **rel_anomaly** (float) - Required - Relativistic anomaly
- **omega_avg** (float) - Required - Orbit-averaged orbital frequency
- **RR** (SEOBNRv5RRForceEcc) - Required - Instance of the radiation reaction force containing an initialized instance of the secular evolution equations
- **params** (EOBParams) - Required - EOB parameters of the system
- **integrator** (str) - Optional - Either 'rkf45' or 'rk8pd' (default). 
- **atol** (float) - Optional - Absolute tolerance of the ODE solver (default 1e-12)
- **rtol** (float) - Optional - Relative tolerance of the ODE solver (default 1e-11)
- **h_0** (float | None) - Optional - Initial time step of the ODE integrator (default -1.0)
- **r_stop** (float | None) - Optional - Separation at which the ODE solver is stopped

### Returns
- **tuple[float, float, float, float]** - (t_start, eccentricity, rel_anomaly, omega_avg), the values of starting time, eccentricity, relativistic anomaly, and orbit-averaged orbital frequency that correspond to r = r_stop
```

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

### MergerRingdownFits.__init__

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/fits/MR_fits.html

Initializes the MergerRingdownFits class with binary system parameters.

```APIDOC
## MergerRingdownFits(m1, m2, chi1, chi2)

### Description
Initializes the class with the masses and dimensionless spin components of the binary system to prepare for coefficient calculations.

### Parameters
- **m1** (float) - Required - Mass of the primary.
- **m2** (float) - Required - Mass of the secondary.
- **chi1** (list) - Required - Dimensionless spin components of the primary.
- **chi2** (list) - Required - Dimensionless spin components of the secondary.
```

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

### Update Spin Parameters and Calibration

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/dynamics/integrate_ode_prec.html

Updates spin parameters and calibration coefficients at the start of the EOB integration process.

```python
X1 = params.p_params.X_1
X2 = params.p_params.X_2

tmp = splines["everything"](omega_start)
chi1_LN_start = tmp[0]
chi2_LN_start = tmp[1]
chi1_L_start = tmp[2]
chi2_L_start = tmp[3]
chi1_v_start = tmp[4:7]
chi2_v_start = tmp[7:10]

params.p_params.omega = omega_start
params.p_params.chi1_v = tuple(chi1_v_start)
params.p_params.chi2_v = tuple(chi2_v_start)

params.p_params.chi_1, params.p_params.chi_2 = chi1_LN_start, chi2_LN_start
params.p_params.chi1_L, params.p_params.chi2_L = chi1_L_start, chi2_L_start

params.p_params.update_spins(chi1_LN_start, chi2_LN_start)

ap_start = chi1_LN_start * X1 + chi2_LN_start * X2
am_start = chi1_LN_start * X1 - chi2_LN_start * X2

dSO_start = dSO_poly_fit(params.p_params.nu, ap_start, am_start)
H.calibration_coeffs.dSO = dSO_start
```

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

### Initialize ODE Integration

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/dynamics/integrate_ode_prec.html

Sets up the GSL ODE integrator and initializes the state variables for the evolution.

```python
sys = odeiv2.system(rhs_wrapper, None, 4, [H, RR, m_1, m_2, params])

T = odeiv2.step_rkf45
s = step(T, 4)
c = control_y_new(atol, rtol)
e = evolve(4)

# Use an agnostic and small initial step for the integrator
h = 0.1

# initial values
t = 0.0
y = y0

# Convert to numpy array to get the correct value
ts = [t]
res_gsl = [y]
omegas = [params.p_params.omega]
augm_dyn = [
    [
        params.p_params.H_val,
        params.p_params.omega,
        params.p_params.omega_circ,
        params.p_params.chi_1,
        params.p_params.chi_2,
    ]
]

# omega_previous = omega_start
r_previous = r0
```

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

### Initialize SEOBNRv5EHM Model Parameters

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/models/SEOBNRv5EHM.html

Initializes model parameters, computes reference shifts, sets Hamiltonian coefficients, and processes GSF contributions.

```python
# Initialize the parameters of the model
self._initialize_params(phys_pars=self.phys_pars)
assert id(self.H.eob_params) == id(self.eob_pars)
self._compute_starting_values()

# Compute the shift from reference point to peak of (2,2) mode
NR_deltaT_fit = NR_deltaT_NS(self.nu) + NR_deltaT(self.nu, self.ap, self.am)
self.NR_deltaT = NR_deltaT_fit

# Set the Hamiltonian coefficients
self._set_H_coeffs()

# Set the GSF contributions to the waveform
gsf_coeffs = GSF_amplitude_fits(self.nu)
keys = gsf_coeffs.keys()

# The following is just a fancy way of passing the coeffs
for key in keys:
    tmp = re.findall(r"h(\d)(\d)_v(\d+)", key)
    if tmp:
        l, m, v = [int(x) for x in tmp[0]]
        self.eob_pars.flux_params.extra_coeffs[l, m, v] = gsf_coeffs[key]
    else:
        tmp = re.findall(r"h(\d)(\d)_vlog(\d+)", key)
        if tmp:
            l, m, v = [int(x) for x in tmp[0]]
            self.eob_pars.flux_params.extra_coeffs_log[l, m, v] = gsf_coeffs[
                key
            ]

self._evaluate_model()
```

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

### Configure secular backwards integration

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_sources/source/notebooks/example_SEOBNRv5EHM.md.txt

Manage secular backwards integration settings, which are active by default for starting separations below 10M.

```python
settings = dict(
    # Setting this to False will avoid the warning message
    # warning_secular_bwd_int=False,
)

t, modes = generate_modes_opt(
    q,
    chi_1,
    chi_2,
    omega_start,
    eccentricity=0.5,
    rel_anomaly=0.0,
    approximant="SEOBNRv5EHM",
    settings=settings,
)
```

```python
plt.figure()
plt.plot(t, np.abs(modes["2,2"]))
plt.xlabel("Time (M)")
plt.ylabel(r"$|h_{2}|$")
plt.axvline(x=0, ls="--", color="k")
plt.grid(True)
```

```python
try:
    settings = dict(
        secular_bwd_int=False,
    )

    t, modes = generate_modes_opt(
        q,
        chi_1,
        chi_2,
        omega_start,
        eccentricity=0.5,
        rel_anomaly=0.0,
        approximant="SEOBNRv5EHM",
        settings=settings,
    )
except Exception as e:
    print(e)
```

```python
settings = dict(
    secular_bwd_int=True,
    r_start_min=7.0,
)

t, modes = generate_modes_opt(
    q,
    chi_1,
    chi_2,
    omega_start,
    eccentricity=0.5,
    rel_anomaly=0.0,
    approximant="SEOBNRv5EHM",
    settings=settings,
)
```

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

### Get all dynamics for anti-symmetric modes

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/waveform/compute_antisymmetric.html

Computes in-plane basis vectors and spin projections required for anti-symmetric mode evaluation.

```python
def get_all_dynamics(
    dyn: dict[str, np.ndarray],
    t: np.ndarray,
    mA: float,
    mB: float,
) -> dict[str, np.ndarray]:
    """Computes the in-plane basis vectors :math:`\boldsymbol{n}` and
    :math:`\boldsymbol{\lambda}` as well
    as the needed projections of the individual spins onto these vectors (and L_N)
    for the evaluation of the anti-symmetric modes.

    :param dyn: dictionary of dynamic quantities on the angles. In particular, the quantities
        ``orbphase``, ``L_N``, ``q_copr``, ``chiA`` and ``chiB`` should be provided and
        should be arrays of time size as ``t``.
    :param t: time array for the anti-symmetric modes, see :py:func:`.compute_time_for_asym`
    :param mA: mass of the primary
    :param mB: mass of the secondary
    :return: dictionary with the needed quantities to evaluate the PN anti-symmetric modes. Contains
        a copy of the keys passed in ``dyn``.

    .. seealso::

        Section §3.2 (and notation section) of [Estelles2025]_.
    """

    orbital_phase = dyn["orbphase"]
    intrp = CubicSpline(t, orbital_phase)
    omega = intrp.derivative()(t)
    v = omega ** (1.0 / 3)

    L_N = dyn["L_N"]
    q = dyn["q_copr"]
    chiA = dyn["chiA"]
    chiB = dyn["chiB"]

    q2 = quaternion.from_float_array(
        np.column_stack(
            (np.cos(orbital_phase / 2), (L_N.T * np.sin(orbital_phase / 2)).T)
        )
    )

    q2_times_q = q2 * q
    x = quaternion.as_float_array(q2_times_q * quaternion.x * q2_times_q.conjugate())[
        :, 1:
    ]
```

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

### Initialize pSEOB support parameters

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/models/SEOBNRv5Base.html

Extends parameter initialization to include deviation dictionaries for plunge-merger and QNM modeling.

```python
def _initialize_params(
    self, *, phys_pars: dict | None, eob_pars: EOBParams | None = None
):
    super()._initialize_params(phys_pars=phys_pars, eob_pars=eob_pars)
    default_deviation_dict = {
        f"{ell},{emm}": 0.0 for ell, emm in self.model_valid_modes
    }

    # Plunge-merger deviations
    self.dA_dict = default_deviation_dict | self.settings.get("dA_dict", {})
    self.dw_dict = default_deviation_dict | self.settings.get("dw_dict", {})

    self.dTpeak = self.settings.get("dTpeak", 0.0)

    # EOB Hamiltonian deviation
    self.da6 = self.settings.get("da6", 0.0)
    self.ddSO = self.settings.get("ddSO", 0.0)

    # QNM deviations
    self.domega_dict = default_deviation_dict | self.settings.get("domega_dict", {})
    self.dtau_dict = default_deviation_dict | self.settings.get("dtau_dict", {})
    for value in self.dtau_dict.values():
        if value <= -1:
            raise ValueError(
                "dtau must be larger than -1, otherwise the remnant rings up instead "
                "of ringing down."
            )

    # Whether one is sampling over the deltaT parameter that determines the merger-ringdown attachment.
    # This does not allow attaching the merger-ringdown at the last point of the dynamics.
    self.deltaT_sampling = self.settings.get("deltaT_sampling", False)
```

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

### Evolve PN Dynamics

Source: https://waveforms.docs.ligo.org/software/pyseobnr/_modules/pyseobnr/eob/dynamics/pn_evolution_opt.html

Integrates the PN equations of motion using solve_ivp, handling both forward and backward integration to reach the desired starting frequency.

```python
        else:
            yout = 1

        return yout

    spinTaylor_stoppingTest.omegas = [omega_ref]
    spinTaylor_stoppingTest.omegaPeaks = 0

    spinTaylor_stoppingTest.terminal = True
    spinTaylor_stoppingTest.direction = -1

    res_forward = solve_ivp(
        prec_eqns_20102022,
        (0, 2e9),
        z0,
        events=spinTaylor_stoppingTest,
        args=(nu, m_1, m_2, X1, X2),
        rtol=rtol,
        atol=atol,
        method=backend,
    )

    if omega_start < omega_ref:
        # We want to start at a lower frequency than omega_ref
        # Thus we first integrate backwards in time
        def term_backwards(t, y, *args):
            return y[-1] - omega_start

        term_backwards.terminal = True

        res_back = solve_ivp(
            prec_eqns_20102022,
            (0, -2e9),
            z0,
            events=term_backwards,
            args=(nu, m_1, m_2, X1, X2),
            rtol=rtol,
            atol=atol,
            method=backend,
        )

        # Combine the forwards and backward pieces
        combined_t = np.concatenate((res_back.t[::-1][:-1], res_forward.t[1:]))
        combined_y = np.vstack((res_back.y.T[::-1][:-1], res_forward.y.T[1:, :]))

        # Set t=0 at the lowest frequency
        combined_t -= combined_t[0]
    else:
        combined_t = res_forward.t[:]
        combined_y = res_forward.y.T[:]

    # Use this to avoid interpolation errors due to repeating the last point
    tdiff_abs = abs(combined_t[-1] - combined_t[-2])
    if tdiff_abs < 1e-10:
        combined_t = combined_t[:-1]
        combined_y = combined_y[:-1]

    omega_diff = combined_y[-1, -1] - combined_y[-2, -1]
    if omega_diff < 0:
        combined_t = combined_t[:-1]
        combined_y = combined_y[:-1]

    return combined_t, combined_y
```
Tokens:
Raw