Source code for pesummary.gw.conversions.tidal

# Licensed under an MIT style license -- see LICENSE.md

import numpy as np
from pesummary.utils.decorators import array_input
from pesummary.utils.utils import logger, iterator

__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"]

try:
    from lalsimulation import (
        CreateSimNeutronStarFamily, SimNeutronStarRadius,
        SimNeutronStarLoveNumberK2, SimNeutronStarEOS4ParameterPiecewisePolytrope,
        SimNSBH_compactness_from_lambda, SimIMRPhenomNSBHProperties,
        SimNeutronStarEOS4ParameterSpectralDecomposition,
        SimIMRPhenomNSBH_baryonic_mass_from_C
    )
    from lal import MRSUN_SI, MSUN_SI
except ImportError:
    pass


[docs]@array_input() def lambda1_plus_lambda2(lambda1, lambda2): """Return the sum of the primary objects tidal deformability and the secondary objects tidal deformability """ return lambda1 + lambda2
[docs]@array_input() def lambda1_minus_lambda2(lambda1, lambda2): """Return the primary objects tidal deformability minus the secondary objests tidal deformability """ return lambda1 - lambda2
[docs]@array_input() def lambda_tilde_from_lambda1_lambda2(lambda1, lambda2, mass1, mass2): """Return the dominant tidal term given samples for lambda1 and lambda2 """ from pesummary.gw.conversions import eta_from_m1_m2 eta = eta_from_m1_m2(mass1, mass2) plus = lambda1_plus_lambda2(lambda1, lambda2) minus = lambda1_minus_lambda2(lambda1, lambda2) lambda_tilde = 8 / 13 * ( (1 + 7 * eta - 31 * eta**2) * plus + (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2) * minus) return lambda_tilde
[docs]@array_input() def delta_lambda_from_lambda1_lambda2(lambda1, lambda2, mass1, mass2): """Return the second dominant tidal term given samples for lambda1 and lambda 2 """ from pesummary.gw.conversions import eta_from_m1_m2 eta = eta_from_m1_m2(mass1, mass2) plus = lambda1_plus_lambda2(lambda1, lambda2) minus = lambda1_minus_lambda2(lambda1, lambda2) delta_lambda = 1 / 2 * ( (1 - 4 * eta) ** 0.5 * (1 - 13272 / 1319 * eta + 8944 / 1319 * eta**2) * plus + (1 - 15910 / 1319 * eta + 32850 / 1319 * eta**2 + 3380 / 1319 * eta**3) * minus) return delta_lambda
[docs]@array_input() def lambda1_from_lambda_tilde(lambda_tilde, mass1, mass2): """Return the individual tidal parameter given samples for lambda_tilde """ from pesummary.gw.conversions import eta_from_m1_m2, q_from_m1_m2 eta = eta_from_m1_m2(mass1, mass2) q = q_from_m1_m2(mass1, mass2) lambda1 = 13 / 8 * lambda_tilde / ( (1 + 7 * eta - 31 * eta**2) * (1 + q**-5) + (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2) * (1 - q**-5)) return lambda1
[docs]@array_input() def lambda2_from_lambda1(lambda1, mass1, mass2): """Return the individual tidal parameter given samples for lambda1 """ from pesummary.gw.conversions import q_from_m1_m2 q = q_from_m1_m2(mass1, mass2) lambda2 = lambda1 / q**5 return lambda2
def _lambda1_lambda2_from_eos(eos, mass_1, mass_2): """Return lambda_1 and lambda_2 assuming a given equation of state """ fam = CreateSimNeutronStarFamily(eos) _lambda = [] for mass in [mass_1, mass_2]: r = SimNeutronStarRadius(mass * MSUN_SI, fam) k = SimNeutronStarLoveNumberK2(mass * MSUN_SI, fam) c = mass * MRSUN_SI / r _lambda.append((2. / 3.) * k / c**5.0) return _lambda
[docs]def wrapper_for_lambda1_lambda2_polytrope_EOS(args): """Wrapper function to calculate the tidal deformability parameters from the 4_parameter_piecewise_polytrope_equation_of_state parameters for a pool of workers """ return _lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state(*args)
def _lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state( log_pressure_si, gamma_1, gamma_2, gamma_3, mass_1, mass_2 ): """Wrapper function to calculate the tidal deformability parameters from the 4_parameter_piecewise_polytrope_equation_of_state parameters for a pool of workers """ eos = SimNeutronStarEOS4ParameterPiecewisePolytrope( log_pressure_si, gamma_1, gamma_2, gamma_3 ) return _lambda1_lambda2_from_eos(eos, mass_1, mass_2)
[docs]def wrapper_for_lambda1_lambda2_from_spectral_decomposition(args): """Wrapper function to calculate the tidal deformability parameters from the spectral decomposition parameters for a pool of workers """ return _lambda1_lambda2_from_spectral_decomposition(*args)
def _lambda1_lambda2_from_spectral_decomposition( spectral_decomposition_gamma_0, spectral_decomposition_gamma_1, spectral_decomposition_gamma_2, spectral_decomposition_gamma_3, mass_1, mass_2 ): """Wrapper function to calculate the tidal deformability parameters from the spectral decomposition parameters for a pool of workers """ gammas = [ spectral_decomposition_gamma_0, spectral_decomposition_gamma_1, spectral_decomposition_gamma_2, spectral_decomposition_gamma_3 ] eos = SimNeutronStarEOS4ParameterSpectralDecomposition(*gammas) return _lambda1_lambda2_from_eos(eos, mass_1, mass_2) def _lambda1_lambda2_from_eos_multiprocess(function, args, multi_process=1): """ """ import multiprocessing with multiprocessing.Pool(multi_process[0]) as pool: lambdas = np.array( list( iterator( pool.imap(function, args), tqdm=True, logger=logger, total=len(args), desc="Calculating tidal parameters" ) ) ) lambdas = np.array(lambdas).T return lambdas[0], lambdas[1]
[docs]@array_input() def lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state( log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2, multi_process=1 ): """Convert 4 parameter piecewise polytrope EOS parameters to the tidal deformability parameters lambda_1, lambda_2 """ logger.warning( "Calculating the tidal deformability parameters based on the 4 " "parameter piecewise polytrope equation of state parameters. This may " "take some time" ) log_pressure_si = log_pressure - 1. args = np.array( [log_pressure_si, gamma_1, gamma_2, gamma_3, mass_1, mass_2] ).T return _lambda1_lambda2_from_eos_multiprocess( wrapper_for_lambda1_lambda2_polytrope_EOS, args, multi_process=multi_process[0] )
[docs]@array_input() def lambda1_lambda2_from_spectral_decomposition( spectral_decomposition_gamma_0, spectral_decomposition_gamma_1, spectral_decomposition_gamma_2, spectral_decomposition_gamma_3, mass_1, mass_2, multi_process=1 ): """Convert spectral decomposition parameters to the tidal deformability parameters lambda_1, lambda_2 """ logger.warning( "Calculating the tidal deformability parameters from the spectral " "decomposition equation of state parameters. This may take some time" ) args = np.array( [ spectral_decomposition_gamma_0, spectral_decomposition_gamma_1, spectral_decomposition_gamma_2, spectral_decomposition_gamma_3, mass_1, mass_2 ] ).T return _lambda1_lambda2_from_eos_multiprocess( wrapper_for_lambda1_lambda2_from_spectral_decomposition, args, multi_process=multi_process[0] )
[docs]@array_input() def lambda1_from_4_parameter_piecewise_polytrope_equation_of_state( log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2 ): """Convert 4 parameter piecewise polytrope EOS parameters to the tidal deformability parameters lambda_1 """ return lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state( log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2 )[0]
[docs]@array_input() def lambda2_from_4_parameter_piecewise_polytrope_equation_of_state( log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2 ): """Convert 4 parameter piecewise polytrope EOS parameters to the tidal deformability parameters lambda_2 """ return lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state( log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2 )[1]
[docs]@array_input() def NS_compactness_from_lambda(lambda_x): """Calculate neutron star compactness from its tidal deformability """ data = np.zeros(len(lambda_x)) for num, _lambda in enumerate(lambda_x): data[num] = SimNSBH_compactness_from_lambda(float(_lambda)) return data
[docs]@array_input() def NS_baryonic_mass(compactness, NS_mass): """Calculate the neutron star baryonic mass from its compactness and gravitational mass in solar masses """ data = np.zeros(len(NS_mass)) for num in np.arange(len(NS_mass)): data[num] = SimIMRPhenomNSBH_baryonic_mass_from_C( compactness[num], NS_mass[num] ) return data
@array_input() def _IMRPhenomNSBH_properties(mass_1, mass_2, spin_1z, lambda_2): """Calculate NSBH specific properties using the IMRPhenomNSBH waveform model given samples for mass_1, mass_2, spin_1z and lambda_2. mass_1 and mass_2 should be in solar mass units """ data = np.zeros((len(mass_1), 6)) for num in range(len(mass_1)): data[num] = SimIMRPhenomNSBHProperties( float(mass_1[num]) * MSUN_SI, float(mass_2[num]) * MSUN_SI, float(spin_1z[num]), float(lambda_2[num]) ) transpose_data = data.T # convert final mass and torus mass to solar masses transpose_data[2] /= MSUN_SI transpose_data[4] /= MSUN_SI return transpose_data def _check_NSBH_approximant(approximant, *args, _raise=True): """Check that the supplied NSBH waveform model is allowed """ if approximant.lower() == "imrphenomnsbh": return _IMRPhenomNSBH_properties(*args) msg = ( "You have supplied the waveform model: '{}'. Currently only the " "IMRPhenomNSBH waveform model can be used. Unable to calculate " "the NSBH conversion".format(approximant) ) if not _raise: logger.warning(msg) else: raise ValueError(msg)
[docs]@array_input() def NSBH_merger_type( mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH", percentages=True, percentage_round=2, _ringdown=None, _disruption=None, _torus=None ): """Determine the merger type based on the disruption frequency, ringdown frequency and torus mass. If percentages = True, a dictionary is returned showing the number of samples which fall in each category. If percentages = False, an array of length mass_1 is returned with elements indicating the merger type for each sample """ _type = np.zeros(len(mass_1), dtype='U15') _type[:] = "disruptive" if not all(param is not None for param in [_ringdown, _disruption, _torus]): ringdown, disruption, torus, _, _, _ = _check_NSBH_approximant( approximant, mass_1, mass_2, spin_1z, lambda_2 ) else: ringdown = _ringdown disruption = _disruption torus = _torus freq_ratio = disruption / ringdown non_disruptive_inds = np.where(freq_ratio > 1) _type[non_disruptive_inds] = "non_disruptive" mildly_disruptive_inds = np.where((freq_ratio < 1) & (torus == 0)) _type[mildly_disruptive_inds] = "mildly_disruptive" if percentages: _percentages = { "non_disruptive": 100 * len(non_disruptive_inds[0]) / len(mass_1), "mildly_disruptive": 100 * len(mildly_disruptive_inds[0]) / len(mass_1) } _percentages["disruptive"] = ( 100 - _percentages["non_disruptive"] - _percentages["mildly_disruptive"] ) for key, value in _percentages.items(): _percentages[key] = np.round(value, percentage_round) return _percentages return _type
[docs]@array_input() def NSBH_ringdown_frequency( mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH" ): """Calculate the ringdown frequency given samples for mass_1, mass_2, spin_1z, lambda_2. mass_1 and mass_2 should be in solar mass units. """ return _check_NSBH_approximant( approximant, mass_1, mass_2, spin_1z, lambda_2 )[0]
[docs]@array_input() def NSBH_tidal_disruption_frequency( mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH" ): """Calculate the tidal disruption frequency given samples for mass_1, mass_2, spin_1z, lambda_2. mass_1 and mass_2 should be in solar mass units. """ return _check_NSBH_approximant( approximant, mass_1, mass_2, spin_1z, lambda_2 )[1]
[docs]@array_input() def NSBH_baryonic_torus_mass( mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH" ): """Calculate the baryonic torus mass given samples for mass_1, mass_2, spin_1z, lambda_2. mass_1 and mass_2 should be in solar mass units. """ return _check_NSBH_approximant( approximant, mass_1, mass_2, spin_1z, lambda_2 )[2]
[docs]@array_input() def NS_compactness_from_NSBH( mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH" ): """Calculate the neutron star compactness given samples for mass_1, mass_2, spin_1z, lambda_2. mass_1 and mass_2 should be in solar mass units. """ return _check_NSBH_approximant( approximant, mass_1, mass_2, spin_1z, lambda_2 )[3]