From c825a43e288ecedeb55143e8e930233fbc2cbc3c Mon Sep 17 00:00:00 2001 From: SamuelH-97 <78215902+SamuelH-97@users.noreply.github.com> Date: Thu, 13 Oct 2022 15:59:05 +0100 Subject: [PATCH] Moving em_progenitors.py and ns_sequence data from tmpltbank to new neutron_stars package (#3937) * Remnant mass function moved from em_progenitors to conversions. Relevant NS functions moved to neutron_stars. Calls to these functions fixed to reflect new locations * no change, pushing to resolve issue with checkout * importing pycbc.tmpltbank rather than pycbc in attempt to resolve import issues (currently unresolved) * moved neutron stars code and eos data to new package and created __init__ file * Renamed em_progenitors as ns_functions and moved to neutron_stars package. Edited how said functions are imported elsewhere where necessary. Added neutron_stars package path to setup file * First suite of changes suggested in PR. Removes some cases of double imports and old code comments from during development * Removed the logging and sys calls. Changed the exception to a ValueError with clear error messages. Co-authored-by: Andrew R. Williamson * retrigger checks * EM-bright work * Constraint string substitutions for static args * Speed up remnant mass calculation (1st attempt) * Adding derivatives in ns_functions for root solver * Add to comment * Missing bracket * Fix bugs * Enable spherical or cartesian spins for m_rem calc * Reorganising neutron_stars module * Some fixes to isso solver * Adding 2D interpolation for ISSO radii * Refined bivariate interpolator for ISSO calc * Change in setup.py, needed for tox * Fix docstrings * Fixing doc strings * Remove old pycbc_dark_vs_bright_injections code * Added newline at end of file to conform to codeclimate check * removed trailing whitespace * removed trailing whitespaces * Removed imports outside of toplevel * removed trailing whitespaces * Empty-Commit, retrigger checks * correctly importing pickle and os.path at top of pg_isso_solver.py * Creating 3 variables in which previously repeated lines of calculation are performed (in need of new names, potenatially better solution available) * fixed issue with unnecessary indents * changed variables to functions after negecting to realise they relied on called variables. Still in need of more informative names, potentially room for more elegant solution * removed trailing whitespaces * Fixing error in trig functions (missing return) * removed functions aimed at avoiding repeated code. Deemed unnecessary * rectified mislabeled parameter chi_lims to bounds as required by function * introducing latest version of remnant mass function in coord_utils get_randm_mass function. Minor rewrite to earlier portion sof function to acommodate this change * added short description of concat_grid function * Moved ns_sequence_file initialisation above it's first use to avoid error * fixed issue with methods used to call variables in remnant mass function * introduced eos as a variable in the get_random_mass function in order to work with updated remnant_mass funtion * adding new lines between function definitions * re-introduced eta_nsbh in get_random_mass as it is required later in the function * moved module level imports to top of init file * moved module level imports back to bottom of init file as moving to top seemed to break code * removed trailing whitespace * Initialised input_is_array variable and added functionality to check if input variable is an array. If it is, input_is_array is set to True * resolved issue with undefined variable erroneously added in code which checks if an input is an array * addrd r prefix to doctring in an attempt to resolve issues related to use of backslashes * Replace code which checks if input is an array with more reliable method Co-authored-by: Tito Dal Canton * Updated docstrings to more accurately represent the forms of inputs and outputs in a number of functions * Fixed typo in docstring Co-authored-by: Tito Dal Canton * Reformatted code to agree with codeclimate Co-authored-by: Tito Dal Canton * Removed hardcoded output path for generate_isso_bivariate_interp() * Fix wrong import name * Updated comments in constraint function Co-authored-by: Francesco Pannarale * Remove ISSO interpolant, will bring back a better version later * Forgot a few things * Update pycbc/tmpltbank/coord_utils.py Co-authored-by: Tito Dal Canton Co-authored-by: Samuel Higginbotham Co-authored-by: Andrew R. Williamson Co-authored-by: Andrew Williamson Co-authored-by: Tito Dal Canton Co-authored-by: Francesco Pannarale Co-authored-by: Tito Dal Canton --- bin/pycbc_create_injections | 2 +- bin/pycbc_dark_vs_bright_injections | 201 --------- pycbc/conversions.py | 97 ++++- pycbc/distributions/constraints.py | 11 +- pycbc/io/record.py | 10 +- pycbc/neutron_stars/__init__.py | 9 + pycbc/neutron_stars/eos_utils.py | 197 +++++++++ .../ns_data}/equil_2H.dat | 0 pycbc/neutron_stars/pg_isso_solver.py | 383 ++++++++++++++++++ pycbc/tmpltbank/__init__.py | 6 - pycbc/tmpltbank/coord_utils.py | 20 +- pycbc/tmpltbank/em_progenitors.py | 360 ---------------- pycbc/tmpltbank/option_utils.py | 2 +- pycbc/waveform/parameters.py | 5 + setup.py | 2 +- 15 files changed, 724 insertions(+), 581 deletions(-) delete mode 100755 bin/pycbc_dark_vs_bright_injections create mode 100644 pycbc/neutron_stars/__init__.py create mode 100644 pycbc/neutron_stars/eos_utils.py rename pycbc/{tmpltbank/ns_sequences => neutron_stars/ns_data}/equil_2H.dat (100%) create mode 100644 pycbc/neutron_stars/pg_isso_solver.py delete mode 100644 pycbc/tmpltbank/em_progenitors.py diff --git a/bin/pycbc_create_injections b/bin/pycbc_create_injections index 0c51a3b0a69..a9616808a94 100644 --- a/bin/pycbc_create_injections +++ b/bin/pycbc_create_injections @@ -206,7 +206,7 @@ dists = distributions.read_distributions_from_config(cp, opts.dist_section) # construct class that will draw the samples randomsampler = JointDistribution(variable_params, *dists, - **{"constraints" : constraints}) + **{"constraints": constraints}) if opts.ninjections: draw_size = opts.ninjections diff --git a/bin/pycbc_dark_vs_bright_injections b/bin/pycbc_dark_vs_bright_injections deleted file mode 100755 index fb310965843..00000000000 --- a/bin/pycbc_dark_vs_bright_injections +++ /dev/null @@ -1,201 +0,0 @@ -#!/usr/bin/env python - -# Copyright (C) 2015 Francesco Pannarale -# -# This program is free software; you can redistribute it and/or modify it -# under the terms of the GNU General Public License as published by the -# Free Software Foundation; either version 3 of the License, or (at your -# option) any later version. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General -# Public License for more details. -# -# You should have received a copy of the GNU General Public License along -# with this program; if not, write to the Free Software Foundation, Inc., -# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - -""" -Given an xml(.gz) table with CBC injections, this script separates them into: -(1) (potentially) EM bright sources, i.e. BNS + (potentially) bright NSBH; -(2) EM dim sources, i.e. BBH + dark NSBH. -The two sets are stored into two separate output files. -""" - -__author__ = "Francesco Pannarale" -__email__ = "francesco.pannarale@ligo.org" -__version__ = "1.1" -__date__ = "28.09.2015" - -import argparse -import logging -import numpy -import pycbc -import pycbc.inject -import pycbc.tmpltbank.em_progenitors -from pycbc.conversions import eta_from_mass1_mass2 -from ligo.lw import utils as ligolw_utils -from ligo.lw import lsctables -from random import sample -from copy import deepcopy - -parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument('-i', dest='inj_xml', required=True, help='Input LIGOLW injections file.') -parser.add_argument('--output-bright', dest='output_bright', required=True, - help="Output LIGOLW file containing potentially EM bright injections.") -parser.add_argument('--output-dim', dest='output_dim', required=True, - help="Output LIGOLW file containing EM dim injections.") -parser.add_argument('--eos', dest='eos', required=True, - help="Select the EOS to be used for the NS when calculating" - "the remnant disk mass. Only 2H is currently supported.") -parser.add_argument('--ns-bh-boundary', type=float, required=True, - help="Mass boundary between neutron stars and black holes. " - "Components below this mass are considered neutron " - "stars. Components at/above are considered black holes." - "UNITS=Solar mass") -parser.add_argument('--remnant-mass-threshold', type=float, required=True, - help="Setting this filters EM dim NS-BH binaries: if the" - "remnant disk mass does not exceed this value, the NS-BH" - "binary is dropped from the bank. UNITS=Solar mass") -parser.add_argument("-z", "--write-compress", action="store_true", help="Write compressed xml.gz files.") -parser.add_argument("--frame-axis-view", action="store_true", - help="The x, y, z spin components are wrt the line of site.") -parser.add_argument("-v", "--verbose", action="store_true", default=False, - help="Extended standard output.") -parser.add_argument("-k", "--max-keep", type=int, default=None, - help="Upper limit on the number of injections in the " - "bright output file.") -opts = parser.parse_args() - -log_fmt = '%(asctime)s %(message)s' -log_date_fmt = '%Y-%m-%d %H:%M:%S' -logging.basicConfig(level=logging.INFO, format=log_fmt, datefmt=log_date_fmt) - -if opts.verbose: - logging.info("Loading injections") -injections = pycbc.inject.InjectionSet(opts.inj_xml) - -if opts.verbose: - logging.info("Loading neutron star equilibrium configurations") -ns_sequence, max_ns_g_mass = pycbc.tmpltbank.em_progenitors.load_ns_sequence(opts.eos) - -# (Potentially) EM bright sources: table and name of the file that will store it -output_bright = opts.output_bright -if opts.write_compress: - if not output_bright.endswith('gz'): - output_bright = output_bright+'.gz' -# make a list of columns that are present in the input table. -# The : split is needed for columns like `process:process_id`, -# which must be listed as `process:process_id` in `lsctables.New()`, -# but are listed as just `process_id` in the `columnnames` attribute -used_columns = [] -for col in injections.table.validcolumns: - att = col.split(':')[-1] - if att in injections.table.columnnames: - used_columns.append(col) - -out_sim_inspiral_bright = lsctables.New(lsctables.SimInspiralTable, - columns=used_columns) - -if opts.max_keep is not None: - out_sim_inspiral_bright_trimmed = deepcopy(out_sim_inspiral_bright) -# EM dim sources: table and name of the file that will store it -output_dim = opts.output_dim -out_sim_inspiral_dim = deepcopy(out_sim_inspiral_bright) -if opts.write_compress: - if not output_dim.endswith('gz'): - output_dim = output_dim+'.gz' - -for i, inj in enumerate(injections.table): - if opts.verbose: - logging.info('%d/%d', i, len(injections.table)) - m1 = inj.mass1 - m2 = inj.mass2 - # BNS are all considered potentially EM bright - if numpy.logical_and(m1 < opts.ns_bh_boundary, m2 < opts.ns_bh_boundary): - out_sim_inspiral_bright.append(inj) - # NSBH systems - elif numpy.logical_not(numpy.logical_and(m1 >= opts.ns_bh_boundary, m2 >= opts.ns_bh_boundary)): - eta = inj.eta if inj.eta is not None else eta_from_mass1_mass2(m1, m2) - s1x = inj.spin1x - s1y = inj.spin1y - s1z = inj.spin1z - s2x = inj.spin2x - s2y = inj.spin2y - s2z = inj.spin2z - incl = inj.inclination - # 1 = NS, 2 = BH - if m1 < m2: - ns_mass = m1 - bh_spin_magnitude = numpy.sqrt(s2x*s2x + s2y*s2y + s2z*s2z) - if opts.frame_axis_view: - # Scalar product between spin and orbital angular momentum - # (the x, y, z spin components are wrt the line of site) - bh_spin_para = s2x*numpy.sin(incl)+s2z*numpy.cos(incl) - else: - # The x, y, z spin components are wrt the orbital angular momentum - bh_spin_para = s2z - bh_spin_inclination = numpy.arccos(bh_spin_para/bh_spin_magnitude) - # 2 = NS, 1 = BH - else: - ns_mass = m2 - bh_spin_magnitude = numpy.sqrt(s1x*s1x + s1y*s1y + s1z*s1z) - if opts.frame_axis_view: - # Scalar product between spin and orbital angular momentum - # (the x, y, z spin components are wrt the line of site) - bh_spin_para = s1x*numpy.sin(incl)+s1z*numpy.cos(incl) - else: - # The x, y, z spin components are wrt the orbital angular momentum - bh_spin_para = s1z - bh_spin_inclination = numpy.arccos(bh_spin_para/bh_spin_magnitude) - # remnant_mass is the remnant disk mass, it is compared to the - # threshold set by the user to discriminate dim and bright - remnant_mass = pycbc.tmpltbank.em_progenitors.remnant_mass( \ - eta, ns_mass, ns_sequence, bh_spin_magnitude, bh_spin_inclination) - if remnant_mass > opts.remnant_mass_threshold: - out_sim_inspiral_bright.append(inj) - else: - out_sim_inspiral_dim.append(inj) - # BBH are all considered EM dark - else: - out_sim_inspiral_dim.append(inj) - -if opts.max_keep is not None and \ - len(out_sim_inspiral_bright) > int(opts.max_keep): - out_sim_inspiral_bright_trimmed.extend(sample(out_sim_inspiral_bright, - int(opts.max_keep))) - logging.info("Found %d/%d (potentially) EM bright injections. Trimming to " - "keep only %d. Storing them to %s." - % (len(out_sim_inspiral_bright), len(injections.table), - len(out_sim_inspiral_bright_trimmed), - output_bright)) -else: - logging.info("Found %d/%d (potentially) EM bright injections. Storing them" - " to %s." % (len(out_sim_inspiral_bright), - len(injections.table), output_bright)) -logging.info("Found %d/%d EM dim injections. Storing them to %s." - % (len(out_sim_inspiral_dim), len(injections.table), output_dim)) - -if opts.verbose: - logging.info('Writing output') -llw_doc = injections.indoc -llw_root = llw_doc.childNodes[0] -llw_root.removeChild(injections.table) -if opts.max_keep is not None and \ - len(out_sim_inspiral_bright) > int(opts.max_keep): - llw_root.appendChild(out_sim_inspiral_bright_trimmed) - ligolw_utils.write_filename(llw_doc, output_bright, - gz=output_bright.endswith('gz')) - llw_root.removeChild(out_sim_inspiral_bright_trimmed) - output_bright = output_bright.replace("POTENTIALLY_BRIGHT", - "POTENTIALLY_BRIGHT_UNTRIMMED") - logging.info(output_bright) -llw_root.appendChild(out_sim_inspiral_bright) -ligolw_utils.write_filename(llw_doc, output_bright, - gz=output_bright.endswith('gz')) -llw_root.removeChild(out_sim_inspiral_bright) -llw_root.appendChild(out_sim_inspiral_dim) -ligolw_utils.write_filename(llw_doc, output_dim, - gz=output_dim.endswith('gz')) -logging.info('Done') diff --git a/pycbc/conversions.py b/pycbc/conversions.py index ebf28142a25..b77abb5a457 100644 --- a/pycbc/conversions.py +++ b/pycbc/conversions.py @@ -23,7 +23,7 @@ # ============================================================================= # """ -This modules provides a library of functions that calculate waveform parameters +This module provides a library of functions that calculate waveform parameters from other parameters. All exposed functions in this module's namespace return one parameter given a set of inputs. """ @@ -33,7 +33,10 @@ import lal from pycbc.detector import Detector import pycbc.cosmology -from .coordinates import spherical_to_cartesian as _spherical_to_cartesian +from .coordinates import ( + spherical_to_cartesian as _spherical_to_cartesian, + cartesian_to_spherical as _cartesian_to_spherical) +from pycbc import neutron_stars as ns pykerr = pycbc.libutils.import_optional('pykerr') lalsim = pycbc.libutils.import_optional('lalsimulation') @@ -417,6 +420,92 @@ def lambda_from_mass_tov_file(mass, tov_file, distance=0.): lambdav = numpy.interp(mass_src, mass_from_file, lambda_from_file) return lambdav + +def remnant_mass_from_mass1_mass2_spherical_spin_eos( + mass1, mass2, spin1a=0.0, spin1pol=0.0, eos='2H'): + """ + Function that determines the remnant disk mass of an NS-BH system + using the fit to numerical-relativity results discussed in + Foucart, Hinderer & Nissanke, PRD 98, 081501(R) (2018). + The BH spin may be misaligned with the orbital angular momentum. + In such cases the ISSO is approximated following the approach of + Stone, Loeb & Berger, PRD 87, 084053 (2013), which was originally + devised for a previous NS-BH remnant mass fit of + Foucart, PRD 86, 124007 (2012). + Note: NS spin is assumed to be 0! + + Parameters + ----------- + mass1 : float + The mass of the black hole, in solar masses. + mass2 : float + The mass of the neutron star, in solar masses. + spin1a : float, optional + The dimensionless magnitude of the spin of mass1. Default = 0. + spin1pol : float, optional + The tilt angle of the spin of mass1. Default = 0 (aligned w L). + eos : str, optional + Name of the equation of state being adopted. Default is '2H'. + + Returns + ---------- + remnant_mass: float + The remnant mass in solar masses + """ + mass1, mass2, spin1a, spin1pol, input_is_array = ensurearray( + mass1, mass2, spin1a, spin1pol) + # mass1 must be greater than mass2 + try: + if any(mass2 > mass1) and input_is_array: + raise ValueError(f'Require mass1 >= mass2') + except TypeError: + if mass2 > mass1 and not input_is_array: + raise ValueError(f'Require mass1 >= mass2. {mass1} < {mass2}') + ns_compactness, ns_b_mass = ns.initialize_eos(mass2, eos) + eta = eta_from_mass1_mass2(mass1, mass2) + remnant_mass = ns.foucart18( + eta, ns_compactness, ns_b_mass, spin1a, spin1pol) + return formatreturn(remnant_mass, input_is_array) + + +def remnant_mass_from_mass1_mass2_cartesian_spin_eos( + mass1, mass2, spin1x=0.0, spin1y=0.0, spin1z=0.0, eos='2H'): + """ + Function that determines the remnant disk mass of an NS-BH system + using the fit to numerical-relativity results discussed in + Foucart, Hinderer & Nissanke, PRD 98, 081501(R) (2018). + The BH spin may be misaligned with the orbital angular momentum. + In such cases the ISSO is approximated following the approach of + Stone, Loeb & Berger, PRD 87, 084053 (2013), which was originally + devised for a previous NS-BH remnant mass fit of + Foucart, PRD 86, 124007 (2012). + Note: NS spin is assumed to be 0! + + Parameters + ----------- + mass1 : float + The mass of the black hole, in solar masses. + mass2 : float + The mass of the neutron star, in solar masses. + spin1x : float, optional + The dimensionless x-component of the spin of mass1. Default = 0. + spin1y : float, optional + The dimensionless y-component of the spin of mass1. Default = 0. + spin1z : float, optional + The dimensionless z-component of the spin of mass1. Default = 0. + eos: str, optional + Name of the equation of state being adopted. Default is '2H'. + + Returns + ---------- + remnant_mass: float + The remnant mass in solar masses + """ + spin1a, _, spin1pol = _cartesian_to_spherical(spin1x, spin1y, spin1z) + return remnant_mass_from_mass1_mass2_spherical_spin_eos( + mass1, mass2, spin1a, spin1pol, eos=eos) + + # # ============================================================================= # @@ -1537,5 +1626,7 @@ def nltides_gw_phase_diff_isco(f_low, f0, amplitude, n, m1, m2): 'optimal_dec_from_detector', 'optimal_ra_from_detector', 'chi_eff_from_spherical', 'chi_p_from_spherical', 'nltides_gw_phase_diff_isco', 'spin_from_pulsar_freq', - 'freqlmn_from_other_lmn', 'taulmn_from_other_lmn' + 'freqlmn_from_other_lmn', 'taulmn_from_other_lmn', + 'remnant_mass_from_mass1_mass2_spherical_spin_eos', + 'remnant_mass_from_mass1_mass2_cartesian_spin_eos' ] diff --git a/pycbc/distributions/constraints.py b/pycbc/distributions/constraints.py index 5272e0178c2..2da1f3478f9 100644 --- a/pycbc/distributions/constraints.py +++ b/pycbc/distributions/constraints.py @@ -16,6 +16,7 @@ This modules provides classes for evaluating multi-dimensional constraints. """ +import re import scipy.spatial import numpy import h5py @@ -37,7 +38,15 @@ def __init__(self, constraint_arg, static_args=None, transforms=None, static_args.items(), key=lambda x: len(x[0]), reverse=True)) ) for arg, val in static_args.items(): - constraint_arg = constraint_arg.replace(arg, str(val)) + swp = f"'{val}'" if isinstance(val, str) else str(val) + # Substitute static arg name for value if it appears in the + # constraint_arg string at the beginning of a word and is not + # followed by an underscore or equals sign. + # This ensures that static_args that are also kwargs in function calls are + # handled correctly, i.e., the kwarg is not touched while its value is replaced + # with the static_arg value. + constraint_arg = re.sub( + r'\b{}(?!\_|\=)'.format(arg), swp, constraint_arg) self.constraint_arg = constraint_arg self.transforms = transforms for kwarg in kwargs.keys(): diff --git a/pycbc/io/record.py b/pycbc/io/record.py index f0353246c0f..27398e56848 100644 --- a/pycbc/io/record.py +++ b/pycbc/io/record.py @@ -1802,7 +1802,8 @@ class WaveformArray(_FieldArrayWithDefaults): parameters.spin_px, parameters.spin_py, parameters.spin_pz, parameters.spin_sx, parameters.spin_sy, parameters.spin_sz, parameters.spin1_a, parameters.spin1_azimuthal, parameters.spin1_polar, - parameters.spin2_a, parameters.spin2_azimuthal, parameters.spin2_polar] + parameters.spin2_a, parameters.spin2_azimuthal, parameters.spin2_polar, + parameters.remnant_mass] @property def primary_mass(self): @@ -1912,5 +1913,12 @@ def spin2_polar(self): return coordinates.cartesian_to_spherical_polar( self.spin2x, self.spin2y, self.spin2z) + @property + def remnant_mass(self): + """Returns the remnant mass for an NS-BH binary.""" + return conversions.remnant_mass_from_mass1_mass2_cartesian_spin_eos( + self.mass1, self.mass2, self.spin1x, + self.spin1y, self.spin1z) + __all__ = ['FieldArray', 'WaveformArray'] diff --git a/pycbc/neutron_stars/__init__.py b/pycbc/neutron_stars/__init__.py new file mode 100644 index 00000000000..33c02d9529d --- /dev/null +++ b/pycbc/neutron_stars/__init__.py @@ -0,0 +1,9 @@ +import os.path +# Setup the directory with the NS equilibrium sequence(s) +NS_DATA_DIRECTORY = os.path.join( + os.path.dirname(__file__), 'ns_data') +NS_SEQUENCES = [ + f.replace('equil_', '').replace('.dat', '') + for f in os.listdir(NS_DATA_DIRECTORY) if f.endswith('.dat')] +from pycbc.neutron_stars.eos_utils import * +from pycbc.neutron_stars.pg_isso_solver import * diff --git a/pycbc/neutron_stars/eos_utils.py b/pycbc/neutron_stars/eos_utils.py new file mode 100644 index 00000000000..a26c8441ca8 --- /dev/null +++ b/pycbc/neutron_stars/eos_utils.py @@ -0,0 +1,197 @@ +# Copyright (C) 2022 Francesco Pannarale, Andrew Williamson, +# Samuel Higginbotham +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +""" +Utility functions for handling NS equations of state +""" +import os.path +import numpy as np +from scipy.interpolate import interp1d +import lalsimulation as lalsim +from . import NS_SEQUENCES, NS_DATA_DIRECTORY +from .pg_isso_solver import PG_ISSO_solver + + +def load_ns_sequence(eos_name): + """ + Load the data of an NS non-rotating equilibrium sequence generated + using the equation of state (EOS) chosen by the user. + File format is: grav mass (Msun), baryonic mass (Msun), compactness + + Parameters + ----------- + eos_name: string + NS equation of state label ('2H' is the only supported + choice at the moment) + + Returns + ---------- + ns_sequence: numpy.array + contains the sequence data in the form NS gravitational + mass (in solar masses), NS baryonic mass (in solar + masses), NS compactness (dimensionless) + max_ns_g_mass: float + the maximum NS gravitational mass (in solar masses) in + the sequence (this is the mass of the most massive stable + NS) + """ + ns_sequence_file = os.path.join( + NS_DATA_DIRECTORY, 'equil_{}.dat'.format(eos_name)) + if eos_name not in NS_SEQUENCES: + raise NotImplementedError( + f'{eos_name} does not have an implemented NS sequence file! ' + f'To implement, the file {ns_sequence_file} must exist and ' + 'contain: NS gravitational mass (in solar masses), NS baryonic ' + 'mass (in solar masses), NS compactness (dimensionless)') + ns_sequence = np.loadtxt(ns_sequence_file) + max_ns_g_mass = max(ns_sequence[:, 0]) + return (ns_sequence, max_ns_g_mass) + + +def interp_grav_mass_to_baryon_mass(ns_g_mass, ns_sequence): + """ + Determines the baryonic mass of an NS given its gravitational + mass and an NS equilibrium sequence (in solar masses). + + Parameters + ----------- + ns_g_mass: float + NS gravitational mass (in solar masses) + ns_sequence: numpy.array + contains the sequence data in the form NS gravitational + mass (in solar masses), NS baryonic mass (in solar + masses), NS compactness (dimensionless) + + Returns + ---------- + float + """ + x = ns_sequence[:, 0] + y = ns_sequence[:, 1] + f = interp1d(x, y) + + return f(ns_g_mass) + + +def interp_grav_mass_to_compactness(ns_g_mass, ns_sequence): + """ + Determines the dimensionless compactness parameter of an NS given + its gravitational mass and an NS equilibrium sequence. + + Parameters + ----------- + ns_g_mass: float + NS gravitational mass (in solar masses) + ns_sequence: numpy.array + contains the sequence data in the form NS gravitational + mass (in solar masses), NS baryonic mass (in solar + masses), NS compactness (dimensionless) + + Returns + ---------- + float + """ + x = ns_sequence[:, 0] + y = ns_sequence[:, 2] + f = interp1d(x, y) + + return f(ns_g_mass) + + +def initialize_eos(ns_mass, eos): + """Load an equation of state and return the compactness and baryonic + mass for a given neutron star mass + + Parameters + ---------- + ns_mass : {float, array} + The mass of the neutron star, in solar masses. + eos : str + Name of the equation of state. + + Returns + ------- + ns_compactness : float + Compactness parameter of the neutron star. + ns_b_mass : float + Baryonic mass of the neutron star. + """ + if isinstance(ns_mass, np.ndarray): + input_is_array = True + if eos in NS_SEQUENCES: + ns_seq, ns_max = load_ns_sequence(eos) + try: + if any(ns_mass > ns_max) and input_is_array: + raise ValueError( + f'Maximum NS mass for {eos} is {ns_max}, received masses ' + f'up to {max(ns_mass[ns_mass > ns_max])}') + except TypeError: + if ns_mass > ns_max and not input_is_array: + raise ValueError( + f'Maximum NS mass for {eos} is {ns_max}, received ' + f'{ns_mass}') + # Interpolate NS compactness and rest mass + ns_compactness = interp_grav_mass_to_compactness(ns_mass, ns_seq) + ns_b_mass = interp_grav_mass_to_baryon_mass(ns_mass, ns_seq) + elif eos in lalsim.SimNeutronStarEOSNames: + #eos_obj = lalsim.SimNeutronStarEOSByName(eos) + #eos_fam = lalsim.CreateSimNeutronStarFamily(eos_obj) + #r_ns = lalsim.SimNeutronStarRadius(ns_mass * lal.MSUN_SI, eos_obj) + #ns_compactness = lal.G_SI * ns_mass * lal.MSUN_SI / (r_ns * lal.C_SI**2) + raise NotImplementedError( + 'LALSimulation EOS interface not yet implemented!') + else: + raise NotImplementedError( + f'{eos} is not implemented! Available are: ' + f'{NS_SEQUENCES + list(lalsim.SimNeutronStarEOSNames)}') + return (ns_compactness, ns_b_mass) + + +def foucart18( + eta, ns_compactness, ns_b_mass, bh_spin_mag, bh_spin_pol): + """Function that determines the remnant disk mass of an NS-BH system + using the fit to numerical-relativity results discussed in + `Foucart, Hinderer & Nissanke, PRD 98, 081501(R) (2018)`_. + + .. _Foucart, Hinderer & Nissanke, PRD 98, 081501(R) (2018): + https://doi.org/10.1103/PhysRevD.98.081501 + + Parameters + ---------- + eta : {float, array} + The symmetric mass ratio of the system + (note: primary is assumed to be the BH). + ns_compactness : {float, array} + NS compactness parameter. + ns_b_mass : {float, array} + Baryonic mass of the NS. + bh_spin_mag: {float, array} + Dimensionless spin magnitude of the BH. + bh_spin_pol : {float, array} + The tilt angle of the BH spin. + """ + isso = PG_ISSO_solver(bh_spin_mag, bh_spin_pol) + # Fit parameters and tidal correction + alpha = 0.406 + beta = 0.139 + gamma = 0.255 + delta = 1.761 + fit = ( + alpha / eta ** (1/3) * (1 - 2 * ns_compactness) + - beta * ns_compactness / eta * isso + + gamma + ) + return ns_b_mass * np.where(fit > 0.0, fit, 0.0) ** delta diff --git a/pycbc/tmpltbank/ns_sequences/equil_2H.dat b/pycbc/neutron_stars/ns_data/equil_2H.dat similarity index 100% rename from pycbc/tmpltbank/ns_sequences/equil_2H.dat rename to pycbc/neutron_stars/ns_data/equil_2H.dat diff --git a/pycbc/neutron_stars/pg_isso_solver.py b/pycbc/neutron_stars/pg_isso_solver.py new file mode 100644 index 00000000000..ea0aab33f09 --- /dev/null +++ b/pycbc/neutron_stars/pg_isso_solver.py @@ -0,0 +1,383 @@ +# Copyright (C) 2022 Francesco Pannarale, Andrew Williamson, +# Samuel Higginbotham +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +""" +Innermost Stable Spherical Orbit (ISSO) solver in the Perez-Giz (PG) +formalism. See `Stone, Loeb, Berger, PRD 87, 084053 (2013)`_. + +.. _Stone, Loeb, Berger, PRD 87, 084053 (2013): + http://dx.doi.org/10.1103/PhysRevD.87.084053 +""" + +import numpy as np +from scipy.optimize import root_scalar + + +def ISCO_eq(r, chi): + r"""Polynomial that enables the calculation of the Kerr innermost + stable circular orbit (ISCO) radius via its roots, + + .. math:: Z(r) = [r (r-6)]^2 - \chi^2 [2r (3r+14) - 9 \chi^2]\,. + + Parameters + ----------- + r: float + the radial coordinate in BH mass units + chi: float + the BH dimensionless spin parameter + + Returns + ---------- + float + """ + return (r * (r - 6))**2 - chi**2 * (2 * r * (3 * r + 14) - 9 * chi**2) + + +def ISCO_eq_dr(r, chi): + """Partial derivative of :func:`ISCO_eq` with respect to r. + + Parameters + ---------- + r: float + the radial coordinate in BH mass units + chi: float + the BH dimensionless spin parameter + + Returns + ------- + float + """ + return 4 * r**3 - 36 * r**2 + 12 * (6 - chi**2) * r - 28 * chi**2 + + +def ISCO_eq_dr2(r, chi): + """Double partial derivative of :func:`ISCO_eq` with respect to r. + + Parameters + ---------- + r: float + the radial coordinate in BH mass units + chi: float + the BH dimensionless spin parameter + + Returns + ------- + float + """ + return 12 * r**2 - 72 * r + 12 * (6 - chi**2) + + +def ISSO_eq_at_pole(r, chi): + r"""Polynomial that enables the calculation of the Kerr polar + (:math:`\iota = \pm \pi / 2`) innermost stable spherical orbit + (ISSO) radius via the roots of + + .. math:: + + P(r) &= r^3 [r^2 (r - 6) + \chi^2 (3 r + 4)] \\ + &\quad + \chi^4 [3 r (r - 2) + \chi^2] \, , + + where :math:`\chi` is the BH dimensionless spin parameter. Physical + solutions are between 6 and + :math:`1 + \sqrt{3} + \sqrt{3 + 2 \sqrt{3}}`. + + Parameters + ---------- + r: float + the radial coordinate in BH mass units + chi: float + the BH dimensionless spin parameter + + Returns + ------- + float + """ + chi2 = chi * chi + return ( + r**3 * (r**2 * (r - 6) + chi2 * (3 * r + 4)) + + chi2 * chi2 * (3 * r * (r - 2) + chi2)) + + +def ISSO_eq_at_pole_dr(r, chi): + """Partial derivative of :func:`ISSO_eq_at_pole` with respect to r. + + Parameters + ---------- + r: float + the radial coordinate in BH mass units + chi: float + the BH dimensionless spin parameter + + Returns + ------- + float + """ + chi2 = chi * chi + twlvchi2 = 12 * chi2 + sxchi4 = 6 * chi2 * chi2 + return ( + 6 * r**5 - 30 * r**4 + twlvchi2 * r**3 + twlvchi2 * r**2 + sxchi4 * r + + sxchi4) + + +def ISSO_eq_at_pole_dr2(r, chi): + """Double partial derivative of :func:`ISSO_eq_at_pole` with + respect to r. + + Parameters + ---------- + r: float + the radial coordinate in BH mass units + chi: float + the BH dimensionless spin parameter + + Returns + ------- + float + """ + chi2 = chi * chi + return ( + 30 * r**4 - 120 * r**3 + 36 * chi2 * r**2 + 24 * chi2 * r + + 6 * chi2 * chi2) + + +def PG_ISSO_eq(r, chi, incl): + r"""Polynomial that enables the calculation of a generic innermost + stable spherical orbit (ISSO) radius via the roots in :math:`r` of + + .. math:: + + S(r) &= r^8 Z(r) + \chi^2 (1 - \cos(\iota)^2) \\ + &\quad * [\chi^2 (1 - \cos(\iota)^2) Y(r) - 2 r^4 X(r)]\,, + + where + + .. math:: + + X(r) &= \chi^2 (\chi^2 (3 \chi^2 + 4 r (2 r - 3)) \\ + &\quad + r^2 (15 r (r - 4) + 28)) - 6 r^4 (r^2 - 4) \, , + + .. math:: + + Y(r) &= \chi^4 (\chi^4 + r^2 [7 r (3 r - 4) + 36]) \\ + &\quad + 6 r (r - 2) \\ + &\qquad * (\chi^6 + 2 r^3 + [\chi^2 (3 r + 2) + 3 r^2 (r - 2)]) \, , + + and :math:`Z(r) =` :func:`ISCO_eq`. Physical solutions are between + the equatorial ISSO (i.e. the ISCO) radius (:func:`ISCO_eq`) and + the polar ISSO radius (:func:`ISSO_eq_at_pole`). + See `Stone, Loeb, Berger, PRD 87, 084053 (2013)`_. + + .. _Stone, Loeb, Berger, PRD 87, 084053 (2013): + http://dx.doi.org/10.1103/PhysRevD.87.084053 + + Parameters + ---------- + r: float + the radial coordinate in BH mass units + chi: float + the BH dimensionless spin parameter + incl: float + inclination angle between the BH spin and the orbital angular + momentum in radians + + Returns + ------- + float + """ + chi2 = chi * chi + chi4 = chi2 * chi2 + r2 = r * r + r4 = r2 * r2 + three_r = 3 * r + r_minus_2 = r - 2 + sin_incl2 = (np.sin(incl))**2 + + X = ( + chi2 * ( + chi2 * (3 * chi2 + 4 * r * (2 * r - 3)) + + r2 * (15 * r * (r - 4) + 28)) + - 6 * r4 * (r2 - 4)) + Y = ( + chi4 * (chi4 + r2 * (7 * r * (three_r - 4) + 36)) + + 6 * r * r_minus_2 * ( + chi4 * chi2 + 2 * r2 * r * ( + chi2 * (three_r + 2) + 3 * r2 * r_minus_2))) + Z = ISCO_eq(r, chi) + + return r4 * r4 * Z + chi2 * sin_incl2 * (chi2 * sin_incl2 * Y - 2 * r4 * X) + + +def PG_ISSO_eq_dr(r, chi, incl): + """Partial derivative of :func:`PG_ISSO_eq` with respect to r. + + Parameters + ---------- + r: float + the radial coordinate in BH mass units + chi: float + the BH dimensionless spin parameter + incl: float + inclination angle between the BH spin and the orbital angular + momentum in radians + + Returns + ------- + float + """ + sini = np.sin(incl) + sin2i = sini * sini + sin4i = sin2i * sin2i + chi2 = chi * chi + chi4 = chi2 * chi2 + chi6 = chi4 * chi2 + chi8 = chi4 * chi4 + chi10 = chi6 * chi4 + return ( + 12 * r**11 - 132 * r**10 + + r**9 * (120 * chi2 * sin2i - 60 * chi2 + 360) - r**8 * 252 * chi2 + + 8 * r**7 * ( + 36 * chi4 * sin4i - 30 * chi4 * sin2i + 9 * chi4 + - 48 * chi2 * sin2i) + + 7 * r**6 * (120 * chi4 * sin2i - 144 * chi4 * sin4i) + + 6 * r**5 * ( + 36 * chi6 * sin4i - 16 * chi6 * sin2i + 144 * chi4 * sin4i + - 56 * chi4 * sin2i) + + r**4 * (120 * chi6 * sin2i - 240 * chi6 * sin4i) + + r**3 * (84 * chi8 * sin4i - 24 * chi8 * sin2i - 192 * chi6 * sin4i) + - 84 * r**2 * chi8 * sin4i + + r * (12 * chi10 * sin4i + 72 * chi8 * sin4i) - 12 * chi10 * sin4i) + + +def PG_ISSO_eq_dr2(r, chi, incl): + """Second partial derivative of :func:`PG_ISSO_eq` with respect to + r. + + Parameters + ---------- + r: float + the radial coordinate in BH mass units + chi: float + the BH dimensionless spin parameter + incl: float + inclination angle between the BH spin and the orbital angular + momentum in radians + + Returns + ------- + float + """ + sini = np.sin(incl) + sin2i = sini * sini + sin4i = sin2i * sin2i + chi2 = chi * chi + chi4 = chi2 * chi2 + chi6 = chi4 * chi2 + chi8 = chi4 * chi4 + return ( + 132 * r**10 - 1320 * r**9 + + 90 * r**8 * (12 * chi2 * sin2i - 6 * chi2 + 36) - 2016 * chi2 * r**7 + + 56 * r**6 * ( + 36 * chi4 * sin4i - 30 * chi4 * sin2i + 9 * chi4 + - 48 * chi2 * sin2i) + + 42 * r**5 * (120 * chi4 * sin2i - 144 * chi4 * sin4i) + + 30 * r**4 * ( + 36 * chi6 * sin4i - 16 * chi6 * sin2i + 144 * chi4 * sin4i + - 56 * chi4 * sin2i) + + r**3 * (480 * chi6 * sin2i - 960 * chi6 * sin4i) + + r**2 * ( + 252 * chi8 * sin4i - 72 * chi8 * sin2i - 576 * chi6 * sin4i) + - r * 168 * chi8 * sin4i + + 12 * chi8 * chi2 * sin4i + 72 * chi8 * sin4i) + + +def PG_ISSO_solver(chi, incl): + """Function that determines the radius of the innermost stable + spherical orbit (ISSO) for a Kerr BH and a generic inclination + angle between the BH spin and the orbital angular momentum. + This function finds the appropriate root of :func:`PG_ISSO_eq`. + + Parameters + ---------- + chi: {float, array} + the BH dimensionless spin parameter + incl: {float, array} + the inclination angle between the BH spin and the orbital + angular momentum in radians + + Returns + ------- + solution: array + the radius of the orbit in BH mass units + """ + # Auxiliary variables + cos_incl = np.cos(incl) + sgnchi = np.sign(cos_incl)*chi + if np.isscalar(sgnchi): + sgnchi = np.array(sgnchi, copy=False, ndmin=1) + chi = np.array(chi, copy=False, ndmin=1) + incl = np.array(incl, copy=False, ndmin=1) + + # ISCO radius for the given spin magnitude + initial_guess = [2 if s > 0.99 else (9 if s < 0 else 5) for s in sgnchi] + rISCO_limit = np.array([ + root_scalar( + ISCO_eq, x0=g0, fprime=ISCO_eq_dr, fprime2=ISCO_eq_dr2, + args=(sc)).root + for g0, sc in zip(initial_guess, sgnchi)]) + # If the inclination is 0 or pi, just output the ISCO radius + equatorial = np.isclose(incl, 0) | np.isclose(incl, np.pi) + if all(equatorial): + return rISCO_limit + + # ISSO radius for an inclination of pi/2 + initial_guess = [9 if c < 0 else 6 for c in chi] + rISSO_at_pole_limit = np.array([ + root_scalar( + ISSO_eq_at_pole, x0=g0, fprime=ISSO_eq_at_pole_dr, + fprime2=ISSO_eq_at_pole_dr2, args=(c)).root + for g0, c in zip(initial_guess, chi)]) + # If the inclination is pi/2, just output the ISSO radius at the pole(s) + polar = np.isclose(incl, np.pi / 2) + if all(polar): + return rISSO_at_pole_limit + + # Otherwise, find the ISSO radius for a generic inclination + initial_hi = np.maximum(rISCO_limit, rISSO_at_pole_limit) + initial_lo = np.minimum(rISCO_limit, rISSO_at_pole_limit) + brackets = [ + (bl, bh) if c == 1 else None + for bl, bh, c in zip(initial_lo, initial_hi, chi)] + solution = np.array([ + root_scalar( + PG_ISSO_eq, x0=g0, fprime=PG_ISSO_eq_dr, bracket=bracket, + fprime2=PG_ISSO_eq_dr2, args=(c, inc), xtol=1e-12).root + for g0, bracket, c, inc in zip(initial_hi, brackets, chi, incl)]) + oob = (solution < 1) | (solution > 9) + if any(oob): + solution = np.array([ + root_scalar( + PG_ISSO_eq, x0=g0, fprime=PG_ISSO_eq_dr, bracket=bracket, + fprime2=PG_ISSO_eq_dr2, args=(c, inc)).root + if ob else sol for g0, bracket, c, inc, ob, sol + in zip(initial_lo, brackets, chi, incl, oob, solution) + ]) + oob = (solution < 1) | (solution > 9) + if any(oob): + raise RuntimeError('Unable to obtain some solutions!') + return solution diff --git a/pycbc/tmpltbank/__init__.py b/pycbc/tmpltbank/__init__.py index 57c1894b2c1..952cdfccb30 100644 --- a/pycbc/tmpltbank/__init__.py +++ b/pycbc/tmpltbank/__init__.py @@ -1,5 +1,3 @@ -import os.path - from pycbc.tmpltbank.calc_moments import * from pycbc.tmpltbank.lambda_mapping import * from pycbc.tmpltbank.coord_utils import * @@ -7,8 +5,4 @@ from pycbc.tmpltbank.brute_force_methods import * from pycbc.tmpltbank.option_utils import * from pycbc.tmpltbank.partitioned_bank import * -from pycbc.tmpltbank.em_progenitors import * from pycbc.tmpltbank.bank_conversions import * - -# Setup the directory with the NS equilibrium sequence(s) -NS_SEQUENCE_FILE_DIRECTORY = os.path.join(os.path.dirname(__file__), 'ns_sequences') diff --git a/pycbc/tmpltbank/coord_utils.py b/pycbc/tmpltbank/coord_utils.py index 06c033ad3e4..c1bb38a7f62 100644 --- a/pycbc/tmpltbank/coord_utils.py +++ b/pycbc/tmpltbank/coord_utils.py @@ -19,7 +19,7 @@ from pycbc.tmpltbank.lambda_mapping import get_chirp_params from pycbc import conversions from pycbc import pnutils -from pycbc.tmpltbank.em_progenitors import load_ns_sequence, remnant_masses +from pycbc.neutron_stars import load_ns_sequence def estimate_mass_range(numPoints, massRangeParams, metricParams, fUpper,\ covary=True): @@ -181,7 +181,7 @@ def get_random_mass_point_particles(numPoints, massRangeParams): return mass1, mass2, spin1z, spin2z -def get_random_mass(numPoints, massRangeParams): +def get_random_mass(numPoints, massRangeParams, eos='2H'): """ This function will generate a large set of points within the chosen mass and spin space, and with the desired minimum remnant disk mass (this applies @@ -195,7 +195,9 @@ def get_random_mass(numPoints, massRangeParams): Number of systems to simulate massRangeParams : massRangeParameters instance Instance holding all the details of mass ranges and spin ranges. - + eos : string + Name of equation of state of neutron star. + Returns -------- mass1 : float @@ -223,7 +225,7 @@ def get_random_mass(numPoints, massRangeParams): # only systems that can yield (at least) the desired remnant # disk mass and that pass the mass and spin range cuts. else: - ns_sequence, max_ns_g_mass = load_ns_sequence(massRangeParams.ns_eos) + max_ns_g_mass = load_ns_sequence(massRangeParams.ns_eos)[1] boundary_mass = massRangeParams.ns_bh_boundary_mass if max_ns_g_mass < boundary_mass: warn_msg = "WARNING: " @@ -291,8 +293,14 @@ def get_random_mass(numPoints, massRangeParams): # threshold mass] mask_bright_nsbh = numpy.zeros(len(mass1_nsbh), dtype=bool) if eta_nsbh.size != 0: - remnant = remnant_masses(eta_nsbh, mass2_nsbh, ns_sequence, - spin1z_nsbh, 0.) + remnant = conversions.remnant_mass_from_mass1_mass2_cartesian_spin_eos( + mass1_nsbh, + mass2_nsbh, + spin1x=0.0, + spin1y=0.0, + spin1z=spin1z_nsbh, + eos=eos + ) mask_bright_nsbh[remnant > massRangeParams.remnant_mass_threshold] = True diff --git a/pycbc/tmpltbank/em_progenitors.py b/pycbc/tmpltbank/em_progenitors.py deleted file mode 100644 index 610b3dbda5b..00000000000 --- a/pycbc/tmpltbank/em_progenitors.py +++ /dev/null @@ -1,360 +0,0 @@ -# Copyright (C) 2015 Francesco Pannarale -# -# This program is free software; you can redistribute it and/or modify it -# under the terms of the GNU General Public License as published by the -# Free Software Foundation; either version 3 of the License, or (at your -# option) any later version. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General -# Public License for more details. -# -# You should have received a copy of the GNU General Public License along -# with this program; if not, write to the Free Software Foundation, Inc., -# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - -import math -import numpy as np -import pycbc -import os.path -import scipy.optimize -import scipy.interpolate -import logging -import sys - -############################################################################## -# Innermost Stable Spherical Orbit (ISSO) solver in the Perez-Giz (PG) # -# formalism [see Stone, Loeb, Berger, PRD 87, 084053 (2013)]. # -############################################################################## - -# Equation that determines the ISCO radius (in BH mass units) -def ISCO_eq(r, chi): - """ - Polynomial that enables the calculation of the Kerr - inntermost stable circular orbit (ISCO) radius via its - roots. - - Parameters - ----------- - r: float - the radial coordinate in BH mass units - chi: float - the BH dimensionless spin parameter - - Returns - ---------- - float - (r*(r-6))**2-chi**2*(2*r*(3*r+14)-9*chi**2) - """ - return (r*(r-6))**2-chi**2*(2*r*(3*r+14)-9*chi**2) - -# Equation that determines the ISSO radius (in BH mass units) at one of the -# poles -def ISSO_eq_at_pole(r, chi): - """ - Polynomial that enables the calculation of the Kerr polar - (inclination = +/- pi/2) innermost stable spherical orbit - (ISSO) radius via its roots. Physical solutions are - between 6 and 1+sqrt[3]+sqrt[3+2sqrt[3]]. - - Parameters - ---------- - r: float - the radial coordinate in BH mass units - chi: float - the BH dimensionless spin parameter - - Returns - ------- - float - r**3*(r**2*(r-6)+chi**2*(3*r+4))+chi**4*(3*r*(r-2)+chi**2) - """ - return r**3*(r**2*(r-6)+chi**2*(3*r+4))+chi**4*(3*r*(r-2)+chi**2) - -# Equation that determines the ISSO radius (in BH mass units) for a generic -# orbital inclination -def PG_ISSO_eq(r, chi, incl): - """Polynomial that enables the calculation of a generic innermost - stable spherical orbit (ISSO) radius via its roots. Physical - solutions are between the equatorial ISSO (aka the ISCO) radius - and the polar ISSO radius. See Stone, Loeb, Berger, PRD 87, 084053 (2013). - - Parameters - ---------- - r: float - the radial coordinate in BH mass units - chi: float - the BH dimensionless spin parameter - incl: float - inclination angle between the BH spin and the orbital angular - momentum in radians - - Returns - ------- - float - ``r**8*Z+chi**2*(1-cos_incl**2)*(chi**2*(1-cos_incl**2)*Y-2*r**4*X)`` - where - ``X=chi**2*(chi**2*(3*chi**2+4*r*(2*r-3))+r**2*(15*r*(r-4)+28))-6*r**4*(r**2-4)`` - ``Y=chi**4*(chi**4+r**2*(7*r*(3*r-4)+36))+6*r*(r-2)*(chi**6+2*r**3*(chi**2*(3*r+2)+3*r**2*(r-2)))`` - ``Z=ISCO_eq(r,chi)`` - """ - chi2 = chi*chi - chi4 = chi2*chi2 - r2 = r*r - r4 = r2*r2 - three_r = 3*r - r_minus_2 = r - 2 - sin_incl2 = (math.sin(incl))**2 - - X = chi2*(chi2*(3*chi2+4*r*(2*r-3))+r2*(15*r*(r-4)+28))-6*r4*(r2-4) - Y = chi4 * (chi4+r2*(7 * r * (three_r-4)+36))+6 * r * r_minus_2 * \ - (chi4*chi2+2*r2*r*(chi2*(three_r+2)+3*r2*r_minus_2)) - Z = ISCO_eq(r, chi) - - return r4*r4*Z+chi2*sin_incl2*(chi2*sin_incl2*Y-2*r4*X) - -# ISSO radius solver -def PG_ISSO_solver(chi,incl): - """Function that determines the radius of the innermost stable - spherical orbit (ISSO) for a Kerr BH and a generic inclination - angle between the BH spin and the orbital angular momentum. - This function finds the appropriat root of PG_ISSO_eq. - - Parameters - ---------- - chi: float - the BH dimensionless spin parameter - incl: float - the inclination angle between the BH spin and the orbital - angular momentum in radians - - Returns - ------- - solution: float - the radius of the orbit in BH mass units - """ - # Auxiliary variables - cos_incl=math.cos(incl) - sgnchi = np.sign(cos_incl)*chi - - # ISCO radius for the given spin magnitude - if sgnchi > 0.99: - initial_guess = 2 # Works better than 6 for chi=1 - elif sgnchi < 0: - initial_guess = 9 - else: - initial_guess = 5 # Works better than 6 for chi=0.5 - rISCO_limit = scipy.optimize.fsolve(ISCO_eq, initial_guess, args=sgnchi) - - # ISSO radius for an inclination of pi/2 - if chi < 0: - initial_guess = 9 - else: - initial_guess = 6 - rISSO_at_pole_limit = scipy.optimize.fsolve(ISSO_eq_at_pole, initial_guess, - args=chi) - - # If the inclination is 0 or pi, just output the ISCO radius - if incl in [0, math.pi]: - solution = rISCO_limit - # If the inclination is pi/2, just output the ISSO radius at the pole(s) - elif incl == math.pi/2: - solution = rISSO_at_pole_limit - # Otherwise, find the ISSO radius for a generic inclination - else: - initial_guess = max(rISCO_limit,rISSO_at_pole_limit) - solution = scipy.optimize.fsolve(PG_ISSO_eq, initial_guess, - args=(chi, incl)) - if solution < 1 or solution > 9: - initial_guess = min(rISCO_limit,rISSO_at_pole_limit) - solution = scipy.optimize.fsolve(PG_ISSO_eq, initial_guess, - args=(chi, incl)) - - return solution - -############################################################################## -# 2H 2-piecewise polytropic EOS, NS non-rotating equilibrium sequence # -# File format is: grav mass (Msun) baryonic mass (Msun) compactness # -# # -# Eventually, the function should call an NS sequence generator within LAL # -# the EOS prescribed by the user and store it. # -############################################################################## -def load_ns_sequence(eos_name): - """ - Load the data of an NS non-rotating equilibrium sequence - generated using the equation of state (EOS) chosen by the - user. [Only the 2H 2-piecewise polytropic EOS is currently - supported. This yields NSs with large radiss (15-16km).] - - Parameters - ----------- - eos_name: string - NS equation of state label ('2H' is the only supported - choice at the moment) - - Returns - ---------- - ns_sequence: 3D-array - contains the sequence data in the form NS gravitational - mass (in solar masses), NS baryonic mass (in solar - masses), NS compactness (dimensionless) - max_ns_g_mass: float - the maximum NS gravitational mass (in solar masses) in - the sequence (this is the mass of the most massive stable - NS) - """ - ns_sequence = [] - - if eos_name == '2H': - ns_sequence_path = os.path.join(pycbc.tmpltbank.NS_SEQUENCE_FILE_DIRECTORY, 'equil_2H.dat') - ns_sequence = np.loadtxt(ns_sequence_path) - else: - err_msg = "Only the 2H EOS is currently supported." - err_msg += "If you plan to use a different NS EOS," - err_msg += "be sure not to filter too many templates!\n" - logging.error(err_msg) - sys.exit(1) - raise Exception('Unsupported EOS!') - - max_ns_g_mass = max(ns_sequence[:,0]) - - return (ns_sequence, max_ns_g_mass) - -############################################################################## -# Given an NS equilibrium sequence and gravitational mass (in Msun), return # -# the NS baryonic mass (in Msun). # -############################################################################## -def ns_g_mass_to_ns_b_mass(ns_g_mass, ns_sequence): - """ - Determines the baryonic mass of an NS given its gravitational - mass and an NS equilibrium sequence. - - Parameters - ----------- - ns_g_mass: float - NS gravitational mass (in solar masses) - ns_sequence: 3D-array - contains the sequence data in the form NS gravitational - mass (in solar masses), NS baryonic mass (in solar - masses), NS compactness (dimensionless) - - Returns - ---------- - float - The NS baryonic mass (in solar massesr**3*(r**2*(r-6)+chi**2*(3*r+4))+ - chi**4*(3*r*(r-2)+chi**2)) - """ - x = ns_sequence[:,0] - y = ns_sequence[:,1] - f = scipy.interpolate.interp1d(x, y) - - return f(ns_g_mass) - -############################################################################## -# Given an NS equilibrium sequence and gravitational mass (in Msun), return # -# the NS compactness. # -############################################################################## -def ns_g_mass_to_ns_compactness(ns_g_mass, ns_sequence): - """ - Determines the compactness of an NS given its - gravitational mass and an NS equilibrium sequence. - - Parameters - ----------- - ns_g_mass: float - NS gravitational mass (in solar masses) - ns_sequence: 3D-array - contains the sequence data in the form NS gravitational - mass (in solar masses), NS baryonic mass (in solar - masses), NS compactness (dimensionless) - - Returns - ---------- - float - The NS compactness (dimensionless) - """ - x = ns_sequence[:,0] - y = ns_sequence[:,2] - f = scipy.interpolate.interp1d(x, y) - - return f(ns_g_mass) - -############################################################################## -# NS-BH merger remnant mass [Foucart, Hinderer, Nissanke PRD 98, 081501(R) # -# (2018)]. # -# # -# * THIS ASSUMES THE NS SPIN IS 0 (the user is warned about this). # -# * Physical parameters passed to remnant_mass (see sanity checks below) # -# must not be unphysical. # -# * Tilted BH spin cases are handled by using rISSO(chi,incl) where rISCO # -# appears in the formula. This approximation was suggested in Stone, Loeb, # -# Berger, PRD 87, 084053 (2013) for the previous NS-BH remnant mass fit of # -# Foucart, PRD 86, 124007 (2012). # -############################################################################## -def remnant_mass(eta, ns_g_mass, ns_sequence, chi, incl): - """ - Function that determines the remnant disk mass of - an NS-BH system using the fit to numerical-relativity - results discussed in Foucart, Hinderer, Nissanke PRD 98, 081501(R) (2018). - - Parameters - ----------- - eta: float - the symmetric mass ratio of the binary - ns_g_mass: float - NS gravitational mass (in solar masses) - ns_sequence: 3D-array - contains the sequence data in the form NS gravitational - mass (in solar masses), NS baryonic mass (in solar - masses), NS compactness (dimensionless) - chi: float - the BH dimensionless spin parameter - incl: float - the inclination angle between the BH spin and the orbital - angular momentum in radians - - Returns - ---------- - remnant_mass: float - The remnant mass in solar masses - """ - - # Sanity checks on eta and chi - if not (eta>0. and eta<=0.25 and abs(chi)<=1 and ns_g_mass>0): - err_msg = "The BH spin magnitude must be <= 1," - err_msg += "eta must be between 0 and 0.25," - err_msg += "and the NS mass must be positive." - logging.error(err_msg) - info_msg = "The function remnant_mass was launched with" - info_msg += "ns_mass={0}, eta={1}, chi={2}," - info_msg += "inclination={3}\n'.format(ns_g_mass, eta, chi, incl)" - logging.info(info_msg) - sys.exit(1) - raise Exception('Unphysical parameters!') - - # NS compactness and rest mass - ns_compactness = ns_g_mass_to_ns_compactness(ns_g_mass, ns_sequence) - ns_b_mass = ns_g_mass_to_ns_b_mass(ns_g_mass, ns_sequence) - - # Fit parameters and tidal correction - alpha = 0.406 - beta = 0.139 - gamma = 0.255 - delta = 1.761 - # The remnant mass over the NS rest mass - remnant_mass = (max(alpha/eta**(1./3.)*(1 - 2*ns_compactness) - - beta*ns_compactness / eta*PG_ISSO_solver(chi, incl) - + gamma, 0.)) ** delta - # Convert to solar masses - remnant_mass = remnant_mass*ns_b_mass - - return remnant_mass - - -############################################################################# -# Vectorized version of remnant_mass. The third argument (NS equilibrium # -# sequence) is excluded from vectorisation. # -############################################################################# -remnant_masses = np.vectorize(remnant_mass) -remnant_masses.excluded.add(2) diff --git a/pycbc/tmpltbank/option_utils.py b/pycbc/tmpltbank/option_utils.py index 39cc0ca085e..8e4a53bf3f8 100644 --- a/pycbc/tmpltbank/option_utils.py +++ b/pycbc/tmpltbank/option_utils.py @@ -21,7 +21,7 @@ import os from pycbc.tmpltbank.lambda_mapping import get_ethinca_orders, pycbcValidOrdersHelpDescriptions from pycbc import pnutils -from pycbc.tmpltbank.em_progenitors import load_ns_sequence +from pycbc.neutron_stars import load_ns_sequence from pycbc.types import positive_float, nonnegative_float class IndentedHelpFormatterWithNL(argparse.ArgumentDefaultsHelpFormatter): diff --git a/pycbc/waveform/parameters.py b/pycbc/waveform/parameters.py index 925fff39a93..47cb9b43b0f 100644 --- a/pycbc/waveform/parameters.py +++ b/pycbc/waveform/parameters.py @@ -383,6 +383,11 @@ def docstr(self, prefix='', include_label=True): numrel_data = Parameter("numrel_data", dtype=str, default="", label=None, description="Sets the NR flags; only needed for NR waveforms.") +remnant_mass = Parameter("remnant_mass", + dtype=float, label=r"$m_{\mathrm{rem}}$", + description="Remnant mass of NS-BH merger. See " + "conversions.remnant_mass_" + "from_mass1_mass2_spin1x_spin1y_spin1z_eos") # # General location parameters diff --git a/setup.py b/setup.py index ed375cb853f..c3cfcb313e2 100755 --- a/setup.py +++ b/setup.py @@ -304,7 +304,7 @@ def run(self): packages = find_packages(), package_data = {'pycbc.workflow': find_files('pycbc/workflow'), 'pycbc.results': find_files('pycbc/results'), - 'pycbc.tmpltbank': find_files('pycbc/tmpltbank')}, + 'pycbc.neutron_stars': find_files('pycbc/neutron_stars')}, ext_modules = ext, python_requires='>=3.7', classifiers=[