Skip to content
Open
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
111 changes: 111 additions & 0 deletions pycbc/constants.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
"""
Constants module for PyCBC.

This module provides constants from various packages (LAL, astropy, numpy)
based on a preferred order of availability. This helps minimize dependency bloat
by relying on smaller, more general packages (astropy/numpy) when the large
LALSuite dependency is not available.

The priority order for constant lookup is: LAL > Astropy / SciPy / NumPy
"""

import logging

# We define a global logger for this module
logger = logging.getLogger('pycbc.constants')

# Attempt to import LALSuite
try:
import lal
_LAL_AVAILABLE = True
logger.debug("LAL constants module found.")
except ImportError:
lal = None
_LAL_AVAILABLE = False
logger.debug("LAL constants module not found.")


import numpy as np
from astropy import (
constants as aconstants,
units as aunits
)
# first, do mappings for constants which the fallback is directly in astropy/numpy
# All are in SI units
_FALLBACK_MAPPING = {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@GarethCabournDavies @spxiwh

In fact, almost everything in this mapping is something we should consider just not using lal for. A few of these may require checking (earth_si, I could imagine some specific value was chosen that is not known to different precision). Most of these though don't need to come from lal in any situation I can imagine.

I'd even prefer we don't have fallbacks (at least not without understanding the potential use case a bit more), because that means there could be unexpected behavior if say lal is or is not the source.

'C_SI': aconstants.c.value, # Speed of light
'G_SI': aconstants.G.value, # Gravitational constant
'MSUN_SI': aconstants.M_sun.value, # Mass of the Sun
'PC_SI': aconstants.pc.value, # Parsec
'REARTH_SI': aconstants.R_earth.value, # Earth equatorial radius
'YRJUL_SI': aunits.year.to(aunits.s), # years in seconds
'PI': np.pi,
'TWOPI': 2 * np.pi,
'GAMMA': np.euler_gamma,
}

# We need to define some constants from astropy values
MSUN_SI = _FALLBACK_MAPPING['MSUN_SI']
C_SI = _FALLBACK_MAPPING['C_SI']
G_SI = _FALLBACK_MAPPING['G_SI']

MTSUN_SI = MSUN_SI * G_SI / (C_SI ** 3)
MRSUN_SI = MSUN_SI * G_SI / (C_SI * C_SI)

# Add in these hybrid constants
_FALLBACK_MAPPING.update({
'MTSUN_SI': MTSUN_SI,
'MRSUN_SI': MRSUN_SI
})

# Define the Constant Lookup Function
def get_constant(name):
"""
Retrieves a constant value by name from the preferred order of packages.

The order of preference is: LAL > Astropy > NumPy

Parameters
----------
name : str
The name of the constant to retrieve (e.g., 'C_SI' for the speed of light).

Returns
-------
float or object
The value of the constant.

Raises
------
NotImplementedError
If the constant is not found in any of the available packages.
"""


if _LAL_AVAILABLE:
return getattr(lal, name)

elif name in _FALLBACK_MAPPING:
return _FALLBACK_MAPPING[name]

raise NotImplementedError(
'Contact the PyCBC team, this should never happen. You are '
'trying to use a LAL constant which doesnt have an astropy/numpy '
'fallback defined.'
)

# Expose the constants as attributes of the module

C_SI = get_constant('C_SI')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@GarethCabournDavies All constants should be gotten from or derived from astropy where it it exists.

Copy link
Member

@ahnitz ahnitz Nov 3, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or thinks like PI should just be gotten from say numpy / scipy as appropriate.

Copy link
Contributor Author

@GarethCabournDavies GarethCabournDavies Nov 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's what this will do 'under the hood' - if I remember right there is some reason (I can't remember exactly what) that lal defines its own PI (and TWOPI) that needs to exact match something on their end. I will see if I can find out why

Or are you saying we should use astropy/numpy/scipy constants as preference over LAL?

Copy link
Contributor

@spxiwh spxiwh Nov 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ahnitz While I agree that it is better in general to use astropy for constants to be consistent, both of the main GW collaborations (IGWN+3G (lumping these together)) and LISA, define constants separately from astropy, and sometimes this can matter. (Cosmology choices tend to be the most obvious gotcha, but LISA positioning might also be quite sensitive as well).

I think it would be good to consider a mechanism where astropy is the default, but you should be able to say "I want constants from LAL" or "I want constants from lisa_constants" ... We could also try to persuade both those communities to be more consistent, but I'm not sure that's going to happen ..

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will have a look at how to implement that. Maybe the pip install pycbc[igwn]-style mechanism may be useful there?

Note: I found that LISA constants for all the bits I could find that were defined by astropy calls straight through anyway

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ahnitz While I agree that it is better in general to use astropy for constants to be consistent, both of the main GW collaborations (IGWN+3G (lumping these together)) and LISA, define constants separately from astropy, and sometimes this can matter. (Cosmology choices tend to be the most obvious gotcha, but LISA positioning might also be quite sensitive as well).

I think it would be good to consider a mechanism where astropy is the default, but you should be able to say "I want constants from LAL" or "I want constants from lisa_constants" ... We could also try to persuade both those communities to be more consistent, but I'm not sure that's going to happen ..

I'm also thinking about using the latest lisa_constants for Mojito Light in PyCBC.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@spxiwh I think I largely agree. That makes sense only when the constants are unique / specialized / specific variations (as you say cosmology is a good example or things specific to an actual detector).

There should be a preference here though to get constants from more standard sources when they do exist. There are a lot of cases where the constants aren't even different numerically and in those cases we should just switch to keep the set of nonstandard constants limited.

G_SI = get_constant('G_SI')

MTSUN_SI = get_constant('MTSUN_SI')
MRSUN_SI = get_constant('MTSUN_SI')
MSUN_SI = get_constant('MSUN_SI')

PC_SI = get_constant('PC_SI')
REARTH_SI = get_constant('REARTH_SI')
YRJUL_SI = get_constant('YRJUL_SI')

PI = get_constant('PI')
GAMMA = get_constant('PI')
36 changes: 18 additions & 18 deletions pycbc/conversions.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,11 +32,11 @@
import numpy
import logging

import lal

from pycbc.detector import Detector
import pycbc.cosmology
from pycbc import neutron_stars as ns
from pycbc.constants import YRJUL_SI, MSUN_SI, MTSUN_SI, C_SI, G_SI, PI

from .coordinates import (
spherical_to_cartesian as _spherical_to_cartesian,
Expand Down Expand Up @@ -99,7 +99,7 @@ def formatreturn(arg, input_is_array=False):

def sec_to_year(sec):
""" Converts number of seconds to number of years """
return sec / lal.YRJUL_SI
return sec / YRJUL_SI


def hypertriangle(*params, bounds=(0, 1)):
Expand Down Expand Up @@ -398,7 +398,7 @@ def tau0_from_mtotal_eta(mtotal, eta, f_lower):
the given frequency.
"""
# convert to seconds
mtotal = mtotal * lal.MTSUN_SI
mtotal = mtotal * MTSUN_SI
# formulae from arxiv.org:0706.4437
return _a0(f_lower) / (mtotal**(5./3.) * eta)

Expand All @@ -407,7 +407,7 @@ def tau0_from_mchirp(mchirp, f_lower):
r"""Returns :math:`\tau_0` from the chirp mass and the given frequency.
"""
# convert to seconds
mchirp = mchirp * lal.MTSUN_SI
mchirp = mchirp * MTSUN_SI
# formulae from arxiv.org:0706.4437
return _a0(f_lower) / mchirp ** (5./3.)

Expand All @@ -417,7 +417,7 @@ def tau3_from_mtotal_eta(mtotal, eta, f_lower):
the given frequency.
"""
# convert to seconds
mtotal = mtotal * lal.MTSUN_SI
mtotal = mtotal * MTSUN_SI
# formulae from arxiv.org:0706.4437
return _a3(f_lower) / (mtotal**(2./3.) * eta)

Expand All @@ -443,7 +443,7 @@ def mchirp_from_tau0(tau0, f_lower):
"""
mchirp = (_a0(f_lower) / tau0) ** (3./5.) # in seconds
# convert back to solar mass units
return mchirp / lal.MTSUN_SI
return mchirp / MTSUN_SI


def mtotal_from_tau0_tau3(tau0, tau3, f_lower,
Expand All @@ -452,7 +452,7 @@ def mtotal_from_tau0_tau3(tau0, tau3, f_lower,
mtotal = (tau3 / _a3(f_lower)) / (tau0 / _a0(f_lower))
if not in_seconds:
# convert back to solar mass units
mtotal /= lal.MTSUN_SI
mtotal /= MTSUN_SI
return mtotal


Expand Down Expand Up @@ -1026,8 +1026,8 @@ def spin_from_pulsar_freq(mass, radius, freq):
The spin frequency of the pulsar, in Hz.
"""
omega = 2 * numpy.pi * freq
mt = mass * lal.MTSUN_SI
mominert = (2/5.) * mt * (radius * 1000 / lal.C_SI)**2
mt = mass * MTSUN_SI
mominert = (2/5.) * mt * (radius * 1000 / C_SI)**2
return mominert * omega / mt**2


Expand Down Expand Up @@ -1403,7 +1403,7 @@ def final_mass_from_f0_tau(f0, tau, l=2, m=2):
# from Berti et al. 2006
spin = final_spin_from_f0_tau(f0, tau, l=l, m=m)
a, b, c = _berti_mass_constants[l,m]
return (a + b*(1-spin)**c)/(2*numpy.pi*f0*lal.MTSUN_SI)
return (a + b*(1-spin)**c)/(2*numpy.pi*f0*MTSUN_SI)

def freqlmn_from_other_lmn(f0, tau, current_l, current_m, new_l, new_m):
"""Returns the QNM frequency (in Hz) of a chosen new (l,m) mode from the
Expand Down Expand Up @@ -1544,14 +1544,14 @@ def get_final_from_initial(mass1, mass2, spin1x=0., spin1y=0., spin1z=0.,
if approximant == 'NRSur7dq4':
from lalsimulation import nrfits
try:
res = nrfits.eval_nrfit(m1*lal.MSUN_SI,
m2*lal.MSUN_SI,
res = nrfits.eval_nrfit(m1*MSUN_SI,
m2*MSUN_SI,
spin1, spin2, 'NRSur7dq4Remnant',
['FinalMass', 'FinalSpin'],
f_ref=f_ref)
except RuntimeError:
continue
final_mass[ii] = res['FinalMass'][0] / lal.MSUN_SI
final_mass[ii] = res['FinalMass'][0] / MSUN_SI
sf = res['FinalSpin']
final_spin[ii] = (sf**2).sum()**0.5
if sf[-1] < 0:
Expand Down Expand Up @@ -1696,7 +1696,7 @@ def velocity_to_frequency(v, M):
f : float
Gravitational-wave frequency
"""
return v**(3.0) / (M * lal.MTSUN_SI * lal.PI)
return v**(3.0) / (M * MTSUN_SI * PI)

def frequency_to_velocity(f, M):
""" Calculate the invariant velocity from the total
Expand All @@ -1714,7 +1714,7 @@ def frequency_to_velocity(f, M):
v : float or numpy.array
Invariant velocity
"""
return (lal.PI * M * lal.MTSUN_SI * f)**(1.0/3.0)
return (PI * M * MTSUN_SI * f)**(1.0/3.0)


def f_schwarzchild_isco(M):
Expand Down Expand Up @@ -1773,13 +1773,13 @@ def nltides_coefs(amplitude, n, m1, m2):

# Calculate chirp mass
mc = mchirp_from_mass1_mass2(m1, m2)
mc *= lal.lal.MSUN_SI
mc *= MSUN_SI

# Calculate constants in phasing
a = (96./5.) * \
(lal.lal.G_SI * lal.lal.PI * mc * f_ref / lal.lal.C_SI**3.)**(5./3.)
(G_SI * PI * mc * f_ref / C_SI**3.)**(5./3.)
b = 6. * amplitude
t_of_f_factor = -1./(lal.lal.PI*f_ref) * b/(a*a * (n-4.))
t_of_f_factor = -1./(PI*f_ref) * b/(a*a * (n-4.))
phi_of_f_factor = -2.*b / (a*a * (n-3.))

return f_ref, t_of_f_factor, phi_of_f_factor
Expand Down
7 changes: 4 additions & 3 deletions pycbc/inject/inject.py
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
from pycbc.filter import resample_to_delta_t
import pycbc.io
from pycbc.io.ligolw import LIGOLWContentHandler
from pycbc.constants import REARTH_SI, C_SI

logger = logging.getLogger('pycbc.inject.inject')

Expand Down Expand Up @@ -226,7 +227,7 @@ def apply(
+ str(strain.dtype))

lalstrain = strain.lal()
earth_travel_time = lal.REARTH_SI / lal.C_SI
earth_travel_time = REARTH_SI / C_SI
t0 = float(strain.start_time) - earth_travel_time
t1 = float(strain.end_time) + earth_travel_time

Expand Down Expand Up @@ -597,7 +598,7 @@ def apply(
t0 = float(strain.start_time)
t1 = float(strain.end_time)
else:
earth_travel_time = lal.REARTH_SI / lal.C_SI
earth_travel_time = REARTH_SI / C_SI
t0 = float(strain.start_time) - earth_travel_time
t1 = float(strain.end_time) + earth_travel_time

Expand Down Expand Up @@ -1338,7 +1339,7 @@ def apply(self, strain, detector_name, f_lower=None, distance_scale=1):

lalstrain = strain.lal()
#detector = Detector(detector_name)
earth_travel_time = lal.REARTH_SI / lal.C_SI
earth_travel_time = REARTH_SI / C_SI
t0 = float(strain.start_time) - earth_travel_time
t1 = float(strain.end_time) + earth_travel_time

Expand Down
4 changes: 2 additions & 2 deletions pycbc/io/gracedb.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@

import pycbc
from pycbc import version as pycbc_version
from pycbc import pnutils
from pycbc import pnutils, constants
from pycbc.io.ligolw import (
return_empty_sngl,
create_process_table,
Expand Down Expand Up @@ -213,7 +213,7 @@ def __init__(self, coinc_ifos, ifos, coinc_results, **kwargs):
coinc_inspiral_row.end_time = sngl_populated.end_time
coinc_inspiral_row.end_time_ns = sngl_populated.end_time_ns
coinc_inspiral_row.snr = network_snrsq ** 0.5
far = 1.0 / (lal.YRJUL_SI * coinc_results['foreground/ifar'])
far = 1.0 / (constants.YRJUL_SI * coinc_results['foreground/ifar'])
coinc_inspiral_row.combined_far = far
coinc_inspiral_table.append(coinc_inspiral_row)
outdoc.childNodes[0].appendChild(coinc_inspiral_table)
Expand Down
Loading
Loading