Coverage for pesummary/gw/conversions/tidal.py: 49.3%
148 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-05-02 08:42 +0000
1# Licensed under an MIT style license -- see LICENSE.md
3import numpy as np
4from pesummary.utils.decorators import array_input
5from pesummary.utils.utils import logger, iterator
7__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"]
9try:
10 from lalsimulation import (
11 CreateSimNeutronStarFamily, SimNeutronStarRadius,
12 SimNeutronStarLoveNumberK2, SimNeutronStarEOS4ParameterPiecewisePolytrope,
13 SimNSBH_compactness_from_lambda, SimIMRPhenomNSBHProperties,
14 SimNeutronStarEOS4ParameterSpectralDecomposition,
15 SimIMRPhenomNSBH_baryonic_mass_from_C
16 )
17 from lal import MRSUN_SI, MSUN_SI
18except ImportError:
19 pass
22@array_input()
23def lambda1_plus_lambda2(lambda1, lambda2):
24 """Return the sum of the primary objects tidal deformability and the
25 secondary objects tidal deformability
26 """
27 return lambda1 + lambda2
30@array_input()
31def lambda1_minus_lambda2(lambda1, lambda2):
32 """Return the primary objects tidal deformability minus the secondary
33 objests tidal deformability
34 """
35 return lambda1 - lambda2
38@array_input()
39def lambda_tilde_from_lambda1_lambda2(lambda1, lambda2, mass1, mass2):
40 """Return the dominant tidal term given samples for lambda1 and lambda2
41 """
42 from pesummary.gw.conversions import eta_from_m1_m2
43 eta = eta_from_m1_m2(mass1, mass2)
44 plus = lambda1_plus_lambda2(lambda1, lambda2)
45 minus = lambda1_minus_lambda2(lambda1, lambda2)
46 lambda_tilde = 8 / 13 * (
47 (1 + 7 * eta - 31 * eta**2) * plus
48 + (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2) * minus)
49 return lambda_tilde
52@array_input()
53def delta_lambda_from_lambda1_lambda2(lambda1, lambda2, mass1, mass2):
54 """Return the second dominant tidal term given samples for lambda1 and
55 lambda 2
56 """
57 from pesummary.gw.conversions import eta_from_m1_m2
58 eta = eta_from_m1_m2(mass1, mass2)
59 plus = lambda1_plus_lambda2(lambda1, lambda2)
60 minus = lambda1_minus_lambda2(lambda1, lambda2)
61 delta_lambda = 1 / 2 * (
62 (1 - 4 * eta) ** 0.5 * (1 - 13272 / 1319 * eta + 8944 / 1319 * eta**2)
63 * plus + (1 - 15910 / 1319 * eta + 32850 / 1319 * eta**2
64 + 3380 / 1319 * eta**3) * minus)
65 return delta_lambda
68@array_input()
69def lambda1_from_lambda_tilde(lambda_tilde, mass1, mass2):
70 """Return the individual tidal parameter given samples for lambda_tilde
71 """
72 from pesummary.gw.conversions import eta_from_m1_m2, q_from_m1_m2
73 eta = eta_from_m1_m2(mass1, mass2)
74 q = q_from_m1_m2(mass1, mass2)
75 lambda1 = 13 / 8 * lambda_tilde / (
76 (1 + 7 * eta - 31 * eta**2) * (1 + q**-5)
77 + (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2) * (1 - q**-5))
78 return lambda1
81@array_input()
82def lambda2_from_lambda1(lambda1, mass1, mass2):
83 """Return the individual tidal parameter given samples for lambda1
84 """
85 from pesummary.gw.conversions import q_from_m1_m2
86 q = q_from_m1_m2(mass1, mass2)
87 lambda2 = lambda1 / q**5
88 return lambda2
91def _lambda1_lambda2_from_eos(eos, mass_1, mass_2):
92 """Return lambda_1 and lambda_2 assuming a given equation of state
93 """
94 fam = CreateSimNeutronStarFamily(eos)
95 _lambda = []
96 for mass in [mass_1, mass_2]:
97 r = SimNeutronStarRadius(mass * MSUN_SI, fam)
98 k = SimNeutronStarLoveNumberK2(mass * MSUN_SI, fam)
99 c = mass * MRSUN_SI / r
100 _lambda.append((2. / 3.) * k / c**5.0)
101 return _lambda
104def wrapper_for_lambda1_lambda2_polytrope_EOS(args):
105 """Wrapper function to calculate the tidal deformability parameters from the
106 4_parameter_piecewise_polytrope_equation_of_state parameters for a pool
107 of workers
108 """
109 return _lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state(*args)
112def _lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state(
113 log_pressure_si, gamma_1, gamma_2, gamma_3, mass_1, mass_2
114):
115 """Wrapper function to calculate the tidal deformability parameters from the
116 4_parameter_piecewise_polytrope_equation_of_state parameters for a pool
117 of workers
118 """
119 eos = SimNeutronStarEOS4ParameterPiecewisePolytrope(
120 log_pressure_si, gamma_1, gamma_2, gamma_3
121 )
122 return _lambda1_lambda2_from_eos(eos, mass_1, mass_2)
125def wrapper_for_lambda1_lambda2_from_spectral_decomposition(args):
126 """Wrapper function to calculate the tidal deformability parameters from
127 the spectral decomposition parameters for a pool of workers
128 """
129 return _lambda1_lambda2_from_spectral_decomposition(*args)
132def _lambda1_lambda2_from_spectral_decomposition(
133 spectral_decomposition_gamma_0, spectral_decomposition_gamma_1,
134 spectral_decomposition_gamma_2, spectral_decomposition_gamma_3,
135 mass_1, mass_2
136):
137 """Wrapper function to calculate the tidal deformability parameters from
138 the spectral decomposition parameters for a pool of workers
139 """
140 gammas = [
141 spectral_decomposition_gamma_0, spectral_decomposition_gamma_1,
142 spectral_decomposition_gamma_2, spectral_decomposition_gamma_3
143 ]
144 eos = SimNeutronStarEOS4ParameterSpectralDecomposition(*gammas)
145 return _lambda1_lambda2_from_eos(eos, mass_1, mass_2)
148def _lambda1_lambda2_from_eos_multiprocess(function, args, multi_process=1):
149 """
150 """
151 import multiprocessing
153 with multiprocessing.Pool(multi_process) as pool:
154 lambdas = np.array(
155 list(
156 iterator(
157 pool.imap(function, args), tqdm=True, logger=logger,
158 total=len(args), desc="Calculating tidal parameters"
159 )
160 )
161 )
162 lambdas = np.array(lambdas).T
163 return lambdas[0], lambdas[1]
166@array_input()
167def lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state(
168 log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2, multi_process=1
169):
170 """Convert 4 parameter piecewise polytrope EOS parameters to the tidal
171 deformability parameters lambda_1, lambda_2
172 """
173 logger.warning(
174 "Calculating the tidal deformability parameters based on the 4 "
175 "parameter piecewise polytrope equation of state parameters. This may "
176 "take some time"
177 )
178 log_pressure_si = log_pressure - 1.
179 args = np.array(
180 [log_pressure_si, gamma_1, gamma_2, gamma_3, mass_1, mass_2]
181 ).T
182 return _lambda1_lambda2_from_eos_multiprocess(
183 wrapper_for_lambda1_lambda2_polytrope_EOS, args,
184 multi_process=multi_process[0]
185 )
188@array_input()
189def lambda1_lambda2_from_spectral_decomposition(
190 spectral_decomposition_gamma_0, spectral_decomposition_gamma_1,
191 spectral_decomposition_gamma_2, spectral_decomposition_gamma_3,
192 mass_1, mass_2, multi_process=1
193):
194 """Convert spectral decomposition parameters to the tidal deformability
195 parameters lambda_1, lambda_2
196 """
197 logger.warning(
198 "Calculating the tidal deformability parameters from the spectral "
199 "decomposition equation of state parameters. This may take some time"
200 )
201 args = np.array(
202 [
203 spectral_decomposition_gamma_0, spectral_decomposition_gamma_1,
204 spectral_decomposition_gamma_2, spectral_decomposition_gamma_3,
205 mass_1, mass_2
206 ]
207 ).T
208 return _lambda1_lambda2_from_eos_multiprocess(
209 wrapper_for_lambda1_lambda2_from_spectral_decomposition, args,
210 multi_process=multi_process[0]
211 )
214@array_input()
215def lambda1_from_4_parameter_piecewise_polytrope_equation_of_state(
216 log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2
217):
218 """Convert 4 parameter piecewise polytrope EOS parameters to the tidal
219 deformability parameters lambda_1
220 """
221 return lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state(
222 log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2
223 )[0]
226@array_input()
227def lambda2_from_4_parameter_piecewise_polytrope_equation_of_state(
228 log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2
229):
230 """Convert 4 parameter piecewise polytrope EOS parameters to the tidal
231 deformability parameters lambda_2
232 """
233 return lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state(
234 log_pressure, gamma_1, gamma_2, gamma_3, mass_1, mass_2
235 )[1]
238@array_input()
239def NS_compactness_from_lambda(lambda_x):
240 """Calculate neutron star compactness from its tidal deformability
241 """
242 data = np.zeros(len(lambda_x))
243 for num, _lambda in enumerate(lambda_x):
244 data[num] = SimNSBH_compactness_from_lambda(float(_lambda))
245 return data
248@array_input()
249def NS_baryonic_mass(compactness, NS_mass):
250 """Calculate the neutron star baryonic mass from its compactness and
251 gravitational mass in solar masses
252 """
253 data = np.zeros(len(NS_mass))
254 for num in np.arange(len(NS_mass)):
255 data[num] = SimIMRPhenomNSBH_baryonic_mass_from_C(
256 compactness[num], NS_mass[num]
257 )
258 return data
261@array_input()
262def _IMRPhenomNSBH_properties(mass_1, mass_2, spin_1z, lambda_2):
263 """Calculate NSBH specific properties using the IMRPhenomNSBH waveform
264 model given samples for mass_1, mass_2, spin_1z and lambda_2. mass_1 and
265 mass_2 should be in solar mass units
266 """
267 data = np.zeros((len(mass_1), 6))
268 for num in range(len(mass_1)):
269 data[num] = SimIMRPhenomNSBHProperties(
270 float(mass_1[num]) * MSUN_SI, float(mass_2[num]) * MSUN_SI,
271 float(spin_1z[num]), float(lambda_2[num])
272 )
273 transpose_data = data.T
274 # convert final mass and torus mass to solar masses
275 transpose_data[2] /= MSUN_SI
276 transpose_data[4] /= MSUN_SI
277 return transpose_data
280def _check_NSBH_approximant(approximant, *args, _raise=True):
281 """Check that the supplied NSBH waveform model is allowed
282 """
283 if approximant.lower() == "imrphenomnsbh":
284 return _IMRPhenomNSBH_properties(*args)
285 msg = (
286 "You have supplied the waveform model: '{}'. Currently only the "
287 "IMRPhenomNSBH waveform model can be used. Unable to calculate "
288 "the NSBH conversion".format(approximant)
289 )
290 if not _raise:
291 logger.warning(msg)
292 else:
293 raise ValueError(msg)
296@array_input()
297def NSBH_merger_type(
298 mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH",
299 percentages=True, percentage_round=2, _ringdown=None, _disruption=None,
300 _torus=None
301):
302 """Determine the merger type based on the disruption frequency, ringdown
303 frequency and torus mass. If percentages = True, a dictionary is returned
304 showing the number of samples which fall in each category. If
305 percentages = False, an array of length mass_1 is returned with
306 elements indicating the merger type for each sample
307 """
308 _type = np.zeros(len(mass_1), dtype='U15')
309 _type[:] = "disruptive"
310 if not all(param is not None for param in [_ringdown, _disruption, _torus]):
311 ringdown, disruption, torus, _, _, _ = _check_NSBH_approximant(
312 approximant, mass_1, mass_2, spin_1z, lambda_2
313 )
314 else:
315 ringdown = _ringdown
316 disruption = _disruption
317 torus = _torus
318 freq_ratio = disruption / ringdown
319 non_disruptive_inds = np.where(freq_ratio > 1)
320 _type[non_disruptive_inds] = "non_disruptive"
321 mildly_disruptive_inds = np.where((freq_ratio < 1) & (torus == 0))
322 _type[mildly_disruptive_inds] = "mildly_disruptive"
323 if percentages:
324 _percentages = {
325 "non_disruptive": 100 * len(non_disruptive_inds[0]) / len(mass_1),
326 "mildly_disruptive": 100 * len(mildly_disruptive_inds[0]) / len(mass_1)
327 }
328 _percentages["disruptive"] = (
329 100 - _percentages["non_disruptive"] - _percentages["mildly_disruptive"]
330 )
331 for key, value in _percentages.items():
332 _percentages[key] = np.round(value, percentage_round)
333 return _percentages
334 return _type
337@array_input()
338def NSBH_ringdown_frequency(
339 mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH"
340):
341 """Calculate the ringdown frequency given samples for mass_1, mass_2,
342 spin_1z, lambda_2. mass_1 and mass_2 should be in solar mass units.
343 """
344 return _check_NSBH_approximant(
345 approximant, mass_1, mass_2, spin_1z, lambda_2
346 )[0]
349@array_input()
350def NSBH_tidal_disruption_frequency(
351 mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH"
352):
353 """Calculate the tidal disruption frequency given samples for mass_1,
354 mass_2, spin_1z, lambda_2. mass_1 and mass_2 should be in solar mass units.
355 """
356 return _check_NSBH_approximant(
357 approximant, mass_1, mass_2, spin_1z, lambda_2
358 )[1]
361@array_input()
362def NSBH_baryonic_torus_mass(
363 mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH"
364):
365 """Calculate the baryonic torus mass given samples for mass_1, mass_2,
366 spin_1z, lambda_2. mass_1 and mass_2 should be in solar mass units.
367 """
368 return _check_NSBH_approximant(
369 approximant, mass_1, mass_2, spin_1z, lambda_2
370 )[2]
373@array_input()
374def NS_compactness_from_NSBH(
375 mass_1, mass_2, spin_1z, lambda_2, approximant="IMRPhenomNSBH"
376):
377 """Calculate the neutron star compactness given samples for mass_1, mass_2,
378 spin_1z, lambda_2. mass_1 and mass_2 should be in solar mass units.
379 """
380 return _check_NSBH_approximant(
381 approximant, mass_1, mass_2, spin_1z, lambda_2
382 )[3]