Commit 0c59df40 authored by Maria Haney's avatar Maria Haney

Merge branch 'nrtidalv2-final' into 'master'

Adding NRTidalv2 approximants

See merge request lscsoft/lalsuite!1200
parents 72da3dac 810d429d
This diff is collapsed.
......@@ -1206,9 +1206,10 @@ static gsl_matrix *XLALSimIMRPhenomBProjectExtrinsicParam(
* (should be used instead of IMRPhenomP,
* unless there are good reasons not to).
* * IMRPhenomPv2_NRTidal models precessing binaries,
* adds NRTides (https://arxiv.org/pdf/1706.02969.pdf)
* adds NR-tuned tidal effects
* to IMRPhenomD phasing and twists the waveform
* to generate the corresponding precessing waveform.
* Two flavors of NRTidal models are available: original (_NRTidal, based on https://arxiv.org/pdf/1706.02969.pdf) and an improved version 2 (_NRTidalv2, based on https://arxiv.org/pdf/1905.06011.pdf).
* * IMRPhenomHM models spinning, non-precessing binaries,
* based on IMRPhenomD that also includes higher order modes.
*
......@@ -1230,6 +1231,10 @@ static gsl_matrix *XLALSimIMRPhenomBProjectExtrinsicParam(
* Riemenschneider, Setyawati, Hinderer. Concluded with git hash
* f15615215a7e70488d32137a827d63192cbe3ef6 (February 2019).
*
* @review IMRPhneomPv2_NRTidalv2 reviewed by Haney, Ossokine, Pürrer. Concluded
* January 2020, for details and final git hash see review wiki:
* https://git.ligo.org/waveforms/reviews/nrtidal_v2/-/wikis/home
*
* @review IMRPhenomHM review wiki page can be found here
* https://git.ligo.org/waveforms/reviews/phenomhm/wikis/home
* @{
......
......@@ -49,7 +49,8 @@ static int IMRPhenomDGenerateFD(
const REAL8 chi1_in, /**< aligned-spin of companion 1 */
const REAL8 chi2_in, /**< aligned-spin of companion 2 */
const REAL8 distance, /**< distance to source (m) */
LALDict *extraParams /**< linked list containing the extra testing GR parameters */
LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
NRTidal_version_type NRTidal_version /**< NRTidal version; either NRTidal_V or NRTidalv2_V or NoNRT_V in case of BBH baseline */
);
/**
......@@ -112,7 +113,8 @@ int XLALSimIMRPhenomDGenerateFD(
const REAL8 f_min, /**< Starting GW frequency (Hz) */
const REAL8 f_max, /**< End frequency; 0 defaults to Mf = \ref f_CUT */
const REAL8 distance, /**< Distance of source (m) */
LALDict *extraParams /**< linked list containing the extra testing GR parameters */
LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
NRTidal_version_type NRTidal_version /**< Version of NRTides; can be one of NRTidal versions or NoNRT_V for the BBH baseline */
) {
/* external: SI; internal: solar masses */
const REAL8 m1 = m1_SI / LAL_MSUN_SI;
......@@ -162,7 +164,7 @@ int XLALSimIMRPhenomDGenerateFD(
freqs->data[1] = f_max_prime;
int status = IMRPhenomDGenerateFD(htilde, freqs, deltaF, phi0, fRef,
m1, m2, chi1, chi2,
distance, extraParams);
distance, extraParams, NRTidal_version);
XLAL_CHECK(XLAL_SUCCESS == status, status, "Failed to generate IMRPhenomD waveform.");
XLALDestroyREAL8Sequence(freqs);
......@@ -204,7 +206,8 @@ int XLALSimIMRPhenomDFrequencySequence(
const REAL8 chi1, /**< Aligned-spin parameter of companion 1 */
const REAL8 chi2, /**< Aligned-spin parameter of companion 2 */
const REAL8 distance, /**< Distance of source (m) */
LALDict *extraParams /**< linked list containing the extra testing GR parameters */
LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
NRTidal_version_type NRTidal_version /**< NRTidal version; either NRTidal_V or NRTidalv2_V or NoNRT_V in case of BBH baseline */
) {
/* external: SI; internal: solar masses */
const REAL8 m1 = m1_SI / LAL_MSUN_SI;
......@@ -232,7 +235,7 @@ int XLALSimIMRPhenomDFrequencySequence(
int status = IMRPhenomDGenerateFD(htilde, freqs, 0, phi0, fRef,
m1, m2, chi1, chi2,
distance, extraParams);
distance, extraParams, NRTidal_version);
XLAL_CHECK(XLAL_SUCCESS == status, status, "Failed to generate IMRPhenomD waveform.");
return XLAL_SUCCESS;
......@@ -261,7 +264,8 @@ static int IMRPhenomDGenerateFD(
const REAL8 chi1_in, /**< aligned-spin of companion 1 */
const REAL8 chi2_in, /**< aligned-spin of companion 2 */
const REAL8 distance, /**< distance to source (m) */
LALDict *extraParams /**< linked list containing the extra testing GR parameters */
LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
NRTidal_version_type NRTidal_version /**< NRTidal version; either NRTidal_V or NRTidalv2_V or NoNRT_V in case of BBH baseline */
) {
LIGOTimeGPS ligotimegps_zero = LIGOTIMEGPSZERO; // = {0, 0}
......@@ -269,18 +273,38 @@ static int IMRPhenomDGenerateFD(
// At the end we will check if we created a LALDict in extraParams
// and destroy it if we did.
LALDict *extraParams_in = extraParams;
REAL8Sequence *amp_tidal = NULL; /* Tidal amplitude series; required only for IMRPhenomD_NRTidalv2 */
REAL8 dquadmon1_in = 0., dquadmon2_in = 0., lambda1_in = 0, lambda2_in = 0.;
if (NRTidal_version == NRTidalv2_V) {
dquadmon1_in = XLALSimInspiralWaveformParamsLookupdQuadMon1(extraParams);
dquadmon2_in = XLALSimInspiralWaveformParamsLookupdQuadMon2(extraParams);
lambda1_in = XLALSimInspiralWaveformParamsLookupTidalLambda1(extraParams);
lambda2_in = XLALSimInspiralWaveformParamsLookupTidalLambda2(extraParams);
}
REAL8 chi1, chi2, m1, m2;
if (m1_in>m2_in) {
REAL8 chi1, chi2, m1, m2, dquadmon1, dquadmon2, lambda1, lambda2;
if (m1_in>=m2_in) {
chi1 = chi1_in;
chi2 = chi2_in;
m1 = m1_in;
m2 = m2_in;
dquadmon1 = dquadmon1_in;
dquadmon2 = dquadmon2_in;
lambda1 = lambda1_in;
lambda2 = lambda2_in;
} else { // swap spins and masses
chi1 = chi2_in;
chi2 = chi1_in;
m1 = m2_in;
m2 = m1_in;
dquadmon1 = dquadmon2_in;
dquadmon2 = dquadmon1_in;
lambda1 = lambda2_in;
lambda2 = lambda1_in;
if (NRTidal_version == NRTidalv2_V) {
XLALSimInspiralWaveformParamsInsertdQuadMon1(extraParams, dquadmon1);
XLALSimInspiralWaveformParamsInsertdQuadMon2(extraParams, dquadmon2);
}
}
int status = init_useful_powers(&powers_of_pi, LAL_PI);
......@@ -398,25 +422,55 @@ static int IMRPhenomDGenerateFD(
const REAL8 phi_precalc = 2.*phi0 + phifRef;
int status_in_for = XLAL_SUCCESS;
int ret = XLAL_SUCCESS;
/* Now generate the waveform */
#pragma omp parallel for
for (UINT4 i=0; i<freqs->length; i++) { // loop over frequency points in sequence
double Mf = M_sec * freqs->data[i];
int j = i + offset; // shift index for frequency series if needed
UsefulPowers powers_of_f;
status_in_for = init_useful_powers(&powers_of_f, Mf);
if (XLAL_SUCCESS != status_in_for)
{
XLALPrintError("init_useful_powers failed for Mf, status_in_for=%d", status_in_for);
status = status_in_for;
if (NRTidal_version == NRTidalv2_V) {
/* Generate the tidal amplitude (Eq. 24 of arxiv: 1905.06011) to add to BBH baseline; only for IMRPhenomD_NRTidalv2 */
amp_tidal = XLALCreateREAL8Sequence(freqs->length);
ret = XLALSimNRTunedTidesFDTidalAmplitudeFrequencySeries(amp_tidal, freqs, m1, m2, lambda1, lambda2);
XLAL_CHECK(XLAL_SUCCESS == ret, ret, "Failed to generate tidal amplitude series to construct IMRPhenomD_NRTidalv2 waveform.");
/* Generated tidal amplitude corrections */
#pragma omp parallel for
for (UINT4 i=0; i<freqs->length; i++) { // loop over frequency points in sequence
double Mf = M_sec * freqs->data[i];
double ampT = amp_tidal->data[i];
int j = i + offset; // shift index for frequency series if needed
UsefulPowers powers_of_f;
status_in_for = init_useful_powers(&powers_of_f, Mf);
if (XLAL_SUCCESS != status_in_for)
{
XLALPrintError("init_useful_powers failed for Mf, status_in_for=%d", status_in_for);
status = status_in_for;
}
else {
REAL8 amp = IMRPhenDAmplitude(Mf, pAmp, &powers_of_f, &amp_prefactors);
REAL8 phi = IMRPhenDPhase(Mf, pPhi, pn, &powers_of_f, &phi_prefactors, 1.0, 1.0);
phi -= t0*(Mf-MfRef) + phi_precalc;
((*htilde)->data->data)[j] = amp0 * (amp+2*sqrt(LAL_PI/5.)*ampT) * cexp(-I * phi);
}
}
else {
REAL8 amp = IMRPhenDAmplitude(Mf, pAmp, &powers_of_f, &amp_prefactors);
REAL8 phi = IMRPhenDPhase(Mf, pPhi, pn, &powers_of_f, &phi_prefactors, 1.0, 1.0);
phi -= t0*(Mf-MfRef) + phi_precalc;
((*htilde)->data->data)[j] = amp0 * amp * cexp(-I * phi);
} else {
#pragma omp parallel for
for (UINT4 i=0; i<freqs->length; i++) { // loop over frequency points in sequence
double Mf = M_sec * freqs->data[i];
int j = i + offset; // shift index for frequency series if needed
UsefulPowers powers_of_f;
status_in_for = init_useful_powers(&powers_of_f, Mf);
if (XLAL_SUCCESS != status_in_for)
{
XLALPrintError("init_useful_powers failed for Mf, status_in_for=%d", status_in_for);
status = status_in_for;
}
else {
REAL8 amp = IMRPhenDAmplitude(Mf, pAmp, &powers_of_f, &amp_prefactors);
REAL8 phi = IMRPhenDPhase(Mf, pPhi, pn, &powers_of_f, &phi_prefactors, 1.0, 1.0);
phi -= t0*(Mf-MfRef) + phi_precalc;
((*htilde)->data->data)[j] = amp0 * amp * cexp(-I * phi);
}
}
}
......@@ -424,6 +478,7 @@ static int IMRPhenomDGenerateFD(
LALFree(pPhi);
LALFree(pn);
XLALDestroyREAL8Sequence(freqs);
XLALDestroyREAL8Sequence(amp_tidal);
/* If extraParams was allocated in this function and not passed in
......
......@@ -44,7 +44,8 @@ static int IMRPhenomD_NRTidal_Core(
REAL8 lambda2, /**< Dimensionless tidal deformability of NS 2 */
LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
const REAL8Sequence *freqs_in, /**< Frequency points at which to evaluate the waveform (Hz) */
REAL8 deltaF /**< Sampling frequency (Hz) */
REAL8 deltaF, /**< Sampling frequency (Hz) */
NRTidal_version_type NRTidal_version /**< Version of NRTides; can be any one of NRTidal_V (arXiv:1706.02969), NRTidalv2_V (arXiv:1905.06011) or NRTidalv2NoAmpCorr_V (arXiv:1905.06011, without amplitude corrections) */
);
// Implementation //////////////////////////////////////////////////////////////
......@@ -62,7 +63,8 @@ int IMRPhenomD_NRTidal_Core(
REAL8 lambda2, /**< Dimensionless tidal deformability of NS 2 */
LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
const REAL8Sequence *freqs_in, /**< Frequency points at which to evaluate the waveform (Hz) */
REAL8 deltaF) /**< Sampling frequency (Hz) */
REAL8 deltaF, /**< Sampling frequency (Hz) */
NRTidal_version_type NRTidal_version /**< Version of NRTides; can be one of NRTidal or NRTidalv2NoAmpCorr */ )
{
/* Check output arrays */
if(!htilde) XLAL_ERROR(XLAL_EFAULT);
......@@ -77,20 +79,36 @@ int IMRPhenomD_NRTidal_Core(
if(fRef == 0.0)
fRef = fLow;
REAL8 dquadmon1 = XLALSimInspiralWaveformParamsLookupdQuadMon1(extraParams);
REAL8 dquadmon2 = XLALSimInspiralWaveformParamsLookupdQuadMon2(extraParams);
/* Internally we need m1 > m2, so change around if this is not the case */
if (m1_SI < m2_SI) {
// Swap m1 and m2
double m1temp = m1_SI;
double chi1temp = chi1;
double lambda1temp = lambda1;
double dquadmon1temp = dquadmon1;
m1_SI = m2_SI;
chi1 = chi2;
lambda1 = lambda2;
m2_SI = m1temp;
chi2 = chi1temp;
lambda2 = lambda1temp;
if (lambda1 != lambda2){
lambda1 = lambda2;
XLALSimInspiralWaveformParamsInsertTidalLambda1(extraParams, lambda1);
lambda2 = lambda1temp;
XLALSimInspiralWaveformParamsInsertTidalLambda2(extraParams, lambda2);
}
if (dquadmon1 != dquadmon2) {
dquadmon1 = dquadmon2;
XLALSimInspiralWaveformParamsInsertdQuadMon1(extraParams, dquadmon1);
dquadmon2 = dquadmon1temp;
XLALSimInspiralWaveformParamsInsertdQuadMon2(extraParams, dquadmon2);
}
}
// Call IMRPhenomD. We call either the FrequencySequence version
// or the regular LAL version depending on how we've been called.
......@@ -122,7 +140,7 @@ int IMRPhenomD_NRTidal_Core(
chi1, chi2,
fLow, f_max_nr_tidal,
distance,
extraParams);
extraParams, NRTidal_version);
// if uniform sampling and fHigh > NRTIDAL_FMAX then resize htilde
// so that it goes up to the user fHigh but is filled with zeros
......@@ -144,13 +162,30 @@ int IMRPhenomD_NRTidal_Core(
m1_SI, m2_SI,
chi1, chi2,
distance,
extraParams);
extraParams, NRTidal_version);
}
XLAL_CHECK(XLAL_SUCCESS == ret, ret, "Failed to generate IMRPhenomD waveform.");
UINT4 offset;
REAL8Sequence *freqs = NULL;
REAL8Sequence *phi_tidal = NULL;
REAL8Sequence *amp_tidal = NULL;
REAL8Sequence *planck_taper = NULL;
/* Initialising terms for HO spin terms added in PhenomD_NRTidalv2 */
REAL8 SS_3p5PN = 0., SSS_3p5PN = 0.;
const REAL8 m1 = m1_SI / LAL_MSUN_SI;
const REAL8 m2 = m2_SI / LAL_MSUN_SI;
const REAL8 M = m1 + m2;
const REAL8 m_sec = M * LAL_MTSUN_SI; /* Total mass in seconds */
REAL8 eta = m1 * m2 / (M*M); /* Symmetric mass-ratio */
const REAL8 piM = LAL_PI * m_sec;
REAL8 X_A = m1/M;
REAL8 X_B = m2/M;
REAL8 pn_fac = 3.*pow(piM,2./3.)/(128.*eta);
/* End of initialising new parameters */
if (deltaF > 0) { // uniform frequencies
// Recreate freqs using only the lower and upper bounds
UINT4 iStart = (UINT4) (fLow / deltaF);
......@@ -172,25 +207,31 @@ int IMRPhenomD_NRTidal_Core(
COMPLEX16 *data=(*htilde)->data->data;
// Get FD tidal phase correction and amplitude factor from arXiv:1706.02969
REAL8Sequence *phi_tidal = XLALCreateREAL8Sequence(freqs->length);
REAL8Sequence *amp_tidal = XLALCreateREAL8Sequence(freqs->length);
ret = XLALSimNRTunedTidesFDTidalPhaseFrequencySeries(
phi_tidal, amp_tidal, freqs,
m1_SI, m2_SI, lambda1, lambda2
);
XLAL_CHECK(XLAL_SUCCESS == ret, ret, "XLALSimNRTunedTidesFDTidalPhaseFrequencySeries Failed.");
phi_tidal = XLALCreateREAL8Sequence(freqs->length);
planck_taper = XLALCreateREAL8Sequence(freqs->length);
if (NRTidal_version == NRTidalv2_V) {
ret = XLALSimNRTunedTidesFDTidalPhaseFrequencySeries(phi_tidal, amp_tidal, planck_taper, freqs, m1_SI, m2_SI, lambda1, lambda2, NRTidalv2NoAmpCorr_V);
XLAL_CHECK(XLAL_SUCCESS == ret, ret, "XLALSimNRTunedTidesFDTidalPhaseFrequencySeries Failed.");
XLALSimInspiralGetHOSpinTerms(&SS_3p5PN, &SSS_3p5PN, X_A, X_B, chi1, chi2, dquadmon1+1., dquadmon2+1.);
}
else {
XLALSimNRTunedTidesFDTidalPhaseFrequencySeries(phi_tidal, amp_tidal, planck_taper, freqs, m1_SI, m2_SI, lambda1, lambda2, NRTidal_version);
XLAL_CHECK(XLAL_SUCCESS == ret, ret, "XLALSimNRTunedTidesFDTidalPhaseFrequencySeries Failed.");
}
// Assemble waveform from amplitude and phase
for (size_t i=0; i<freqs->length; i++) { // loop over frequency points in sequence
int j = i + offset; // shift index for frequency series if needed
// Apply tidal phase correction and amplitude taper
COMPLEX16 Corr = amp_tidal->data[i] * cexp(-I*phi_tidal->data[i]);
double f = freqs->data[i];
COMPLEX16 Corr = planck_taper->data[i] * cexp(-I*phi_tidal->data[i] - I*pn_fac*(SS_3p5PN + SSS_3p5PN)*pow(f,2./3.));
data[j] *= Corr;
}
XLALDestroyREAL8Sequence(freqs);
XLALDestroyREAL8Sequence(phi_tidal);
XLALDestroyREAL8Sequence(amp_tidal);
XLALDestroyREAL8Sequence(planck_taper);
return XLAL_SUCCESS;
}
......@@ -254,14 +295,15 @@ int XLALSimIMRPhenomDNRTidalFrequencySequence(
REAL8 chi2, /**< Dimensionless aligned component spin of NS 2 */
REAL8 lambda1, /**< Dimensionless tidal deformability of NS 1 */
REAL8 lambda2, /**< Dimensionless tidal deformability of NS 2 */
LALDict *extraParams /**< linked list containing the extra testing GR parameters */
LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
NRTidal_version_type NRTidal_version /**< Version of NRTides; can be any one of NRTidal_V (arXiv:1706.02969), NRTidalv2_V (arXiv:1905.06011) or NRTidalv2NoAmpCorr_V (arXiv:1905.06011, without amplitude corrections) */
) {
if (!freqs) XLAL_ERROR(XLAL_EFAULT);
// Call the internal core function with deltaF = 0 to indicate that freqs is non-uniformly
// spaced and we want the strain only at these frequencies
int retcode = IMRPhenomD_NRTidal_Core(htilde,
phiRef, fRef, distance, m1_SI, m2_SI, chi1, chi2, lambda1, lambda2, extraParams, freqs, 0);
phiRef, fRef, distance, m1_SI, m2_SI, chi1, chi2, lambda1, lambda2, extraParams, freqs, 0, NRTidal_version);
return(retcode);
}
......@@ -289,7 +331,8 @@ int XLALSimIMRPhenomDNRTidal(
REAL8 chi2, /**< Dimensionless aligned component spin of NS 2 */
REAL8 lambda1, /**< Dimensionless tidal deformability of NS 1 */
REAL8 lambda2, /**< Dimensionless tidal deformability of NS 2 */
LALDict *extraParams /**< linked list containing the extra testing GR parameters */
LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
NRTidal_version_type NRTidal_version /**< Version of NRTides; can be any one of NRTidal_V (arXiv:1706.02969), NRTidalv2_V (arXiv:1905.06011) or NRTidalv2NoAmpCorr_V (arXiv:1905.06011, without amplitude corrections) */
) {
// Use fLow, fHigh, deltaF to compute freqs sequence
// Instead of building a full sequence we only transfer the boundaries and let
......@@ -299,7 +342,7 @@ int XLALSimIMRPhenomDNRTidal(
freqs->data[1] = fHigh;
int retcode = IMRPhenomD_NRTidal_Core(htilde,
phiRef, fRef, distance, m1_SI, m2_SI, chi1, chi2, lambda1, lambda2, extraParams, freqs, deltaF);
phiRef, fRef, distance, m1_SI, m2_SI, chi1, chi2, lambda1, lambda2, extraParams, freqs, deltaF, NRTidal_version);
XLALDestroyREAL8Sequence(freqs);
......
This diff is collapsed.
......@@ -98,6 +98,7 @@ static int PhenomPCore(
* spacing deltaF. Otherwise, the frequency points are spaced non-uniformly.
* Then we will use deltaF = 0 to create the frequency series we return. */
IMRPhenomP_version_type IMRPhenomP_version, /**< IMRPhenomPv1 uses IMRPhenomC, IMRPhenomPv2 uses IMRPhenomD, IMRPhenomPv2_NRTidal is a tidal version of IMRPhenomPv2 */
NRTidal_version_type NRTidal_version, /**< either NRTidal or NRTidalv2 for BNS waveform; NoNRT_V for BBH waveform */
LALDict *extraParams /**< linked list containing the extra testing GR parameters */
);
......@@ -201,8 +202,9 @@ static void nudge(REAL8 *x, REAL8 X, REAL8 epsilon);
static int PhenomPCoreOneFrequency_withTides(
const REAL8 fHz, /**< Frequency (Hz) */
const REAL8 ampTidal, /**< tidal amplitude at a frequency sample; planck window */
COMPLEX16 phaseTidal, /**< tidal phasing at a frequency sample from NRTidal infrastructure*/
const REAL8 window, /**< planck window */
const REAL8 ampTidal, /**< tidal amplitude at a frequency sample */
const REAL8 phaseTidal, /**< tidal phasing at a frequency sample from NRTidal infrastructure*/
const REAL8 distance, /**< Distance of source (m) */
const REAL8 M, /**< Total mass (Solar masses) */
const REAL8 phic, /**< Orbital phase at the peak of the underlying non precessing model (rad) */
......
This diff is collapsed.
#ifndef _LALSIM_IMR_SEOBNRv4_ROM_NRTidal_H
#define _LALSIM_IMR_SEOBNRv4_ROM_NRTidal_H
static void Self_spin_phase_contributions(
const REAL8 m1_SI, /**< Mass of neutron star 1 (kg) */
const REAL8 m2_SI, /**< Mass of neutron star 2 (kg) */
const REAL8 chi1L, /**< Dimensionless aligned component spin of NS 1 */
const REAL8 chi2L, /**< Dimensionless aligned component spin of NS 2 */
const REAL8 qm_def1, /**< Quadrupole deformation parameter of body 1 (dimensionless) */
const REAL8 qm_def2, /**< Quadrupole deformation parameter of body 2 (dimensionless) */
/**< qm_def1,2 = 1 for BH */
REAL8 *pfa_v4_contrib, /**< self-spin contribution to v^4 */
REAL8 *pfa_v6_contrib /**< self-spin contribution to v^6 */
);
static int SEOBNRv4ROM_NRTidal_Core(
struct tagCOMPLEX16FrequencySeries **hptilde, /**< Output: Frequency-domain waveform h+ */
struct tagCOMPLEX16FrequencySeries **hctilde, /**< Output: Frequency-domain waveform hx */
REAL8 phiRef, /**< Phase at reference time */
REAL8 fRef, /**< Reference frequency (Hz); 0 defaults to fLow */
REAL8 distance, /**< Distance of source (m) */
REAL8 inclination, /**< Inclination of source (rad) */
REAL8 m1_SI, /**< Mass of neutron star 1 (kg) */
REAL8 m2_SI, /**< Mass of neutron star 2 (kg) */
REAL8 chi1, /**< Dimensionless aligned component spin of NS 1 */
REAL8 chi2, /**< Dimensionless aligned component spin of NS 2 */
REAL8 Lambda1, /**< Dimensionless tidal deformability of NS 1 */
REAL8 Lambda2, /**< Dimensionless tidal deformability of NS 2 */
const REAL8Sequence *freqs_in, /**< Frequency points at which to evaluate the waveform (Hz) */
REAL8 deltaF /**< Sampling frequency (Hz) */
);
struct tagCOMPLEX16FrequencySeries **hptilde, /**< Output: Frequency-domain waveform h+ */
struct tagCOMPLEX16FrequencySeries **hctilde, /**< Output: Frequency-domain waveform hx */
REAL8 phiRef, /**< Phase at reference time */
REAL8 fRef, /**< Reference frequency (Hz); 0 defaults to fLow */
REAL8 distance, /**< Distance of source (m) */
REAL8 inclination, /**< Inclination of source (rad) */
REAL8 m1_SI, /**< Mass of neutron star 1 (kg) */
REAL8 m2_SI, /**< Mass of neutron star 2 (kg) */
REAL8 chi1, /**< Dimensionless aligned component spin of NS 1 */
REAL8 chi2, /**< Dimensionless aligned component spin of NS 2 */
REAL8 Lambda1, /**< Dimensionless tidal deformability of NS 1 */
REAL8 Lambda2, /**< Dimensionless tidal deformability of NS 2 */
const REAL8Sequence *freqs_in, /**< Frequency points at which to evaluate the waveform (Hz) */
REAL8 deltaF, /**< Sampling frequency (Hz) */
LALDict *LALparams, /**< linked list containing the extra testing GR parameters */
NRTidal_version_type NRTidal_version /** < either NRTidal or NRTidalv2NoAmpCorr */
);
#endif
This diff is collapsed.
......@@ -361,6 +361,8 @@ typedef enum tagApproximant {
* @remarks Implemented in lalsimulation (frequency domain). */
SEOBNRv4_ROM_NRTidal, /**< Low-mass double-spin frequency domain reduced order model of spin-aligned EOBNR model SEOBNRv4 [Bohe et al, arXiv:1611.03703] with tidal phase corrections [Dietrich et al, arXiv:1706.02969]
* @remarks Implemented in lalsimulation (frequency domain). */
SEOBNRv4_ROM_NRTidalv2, /**< based on NRTidalv2; https://arxiv.org/abs/1905.06011.
* @remarks Implemented in lalsimulation (time domain and frequency domain). */
SEOBNRv4T_surrogate, /**< Double-spin frequency domain surrogate model of spin-aligned tidal EOBNR model SEOBNRv4T
* @remarks Implemented in lalsimulation (frequency domain). */
HGimri, /**< Time domain inspiral-merger-ringdown waveform for quasi-circular intermediate mass-ratio inspirals [Huerta & Gair arXiv:1009.1985]
......@@ -379,6 +381,8 @@ typedef enum tagApproximant {
* @remarks Implemented in lalsimulation (frequency domain). */
IMRPhenomD_NRTidal, /**< Uses arxiv:1706.02969 to upgrad IMRPhenomD to a tidal approximant
* @remarks Implemented in lalsimulation (frequency domain). */
IMRPhenomD_NRTidalv2, /**< NRTidalv2; https://arxiv.org/abs/1905.06011
* @remarks Implemented in lalsimulation (time domain and frequency domain).*/
IMRPhenomHM, /**< Frequency domain with higher modes (non-precessing spins) inspiral-merger-ringdown templates, based on IMRPhenomD.
* @remarks Implemented in lalsimulation (frequency domain). Ref London et al, arXiv:1708.00404 */
IMRPhenomP, /**< Frequency domain (generic spins) inspiral-merger-ringdown templates of Hannam et al., arXiv:1308.3271 [gr-qc]. Based on IMRPhenomC.
......@@ -386,6 +390,8 @@ typedef enum tagApproximant {
IMRPhenomPv2, /**< Frequency domain (generic spins) inspiral-merger-ringdown templates of Hannam et al., arXiv:1308.3271 [gr-qc]. Based on IMRPhenomD, arXiv:1508.07250 and arXiv:1508.07253.
* @remarks Implemented in lalsimulation (frequency domain). */
IMRPhenomPv2_NRTidal, /**< Frequency domain tidal version of IMRPhenomPv2, using NRTidal framework from arXiv:1706.02969 */
IMRPhenomPv2_NRTidalv2, /**< Frequency domain tidal version; based on https://arxiv.org/abs/1905.06011
* @remarks Implemented in lalsimulation (time domain and frequency domain).*/
IMRPhenomFC, /**< Frequency domain (non-precessing spins) inspiral-merger-ringdown templates of Santamaria et al [Santamaria:2010yb] with phenomenological coefficients defined in the Table II of [Santamaria:2010yb]
* @attention Not implemented in lalsimulation.*/
TaylorEt, /**< UNDOCUMENTED
......
......@@ -1042,7 +1042,7 @@ int XLALSimInspiralChooseFDWaveformSequence(
ABORT_NONZERO_TIDES(LALpars);
ret = XLALSimIMRSEOBNRv4ROMFrequencySequence(hptilde, hctilde, frequencies,
phiRef, f_ref, distance, inclination, m1, m2, S1z, S2z, -1);
phiRef, f_ref, distance, inclination, m1, m2, S1z, S2z, -1, LALpars, NoNRT_V);
break;
case SEOBNRv4_ROM_NRTidal:
......@@ -1052,8 +1052,25 @@ int XLALSimInspiralChooseFDWaveformSequence(
if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
ABORT_NONZERO_TRANSVERSE_SPINS(LALpars);
ret = XLALSimInspiralSetQuadMonParamsFromLambdas(LALpars);
XLAL_CHECK(XLAL_SUCCESS == ret, ret, "Failed to set QuadMon from Lambdas for SEOBNRv4_ROM_NRTidal");
ret = XLALSimIMRSEOBNRv4ROMNRTidalFrequencySequence(hptilde, hctilde, frequencies,
phiRef, f_ref, distance, inclination, m1, m2, S1z, S2z, lambda1, lambda2);
phiRef, f_ref, distance, inclination, m1, m2, S1z, S2z, lambda1, lambda2, LALpars, NRTidal_V);
break;
case SEOBNRv4_ROM_NRTidalv2:
/* Waveform-specific sanity checks */
if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALpars) )
ABORT_NONDEFAULT_LALDICT_FLAGS(LALpars);
if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
ABORT_NONZERO_TRANSVERSE_SPINS(LALpars);
ret = XLALSimInspiralSetQuadMonParamsFromLambdas(LALpars);
XLAL_CHECK(XLAL_SUCCESS == ret, ret, "Failed to set QuadMon from Lambdas for SEOBNRv4_ROM_NRTidalv2");
ret = XLALSimIMRSEOBNRv4ROMNRTidalFrequencySequence(hptilde, hctilde, frequencies,
phiRef, f_ref, distance, inclination, m1, m2, S1z, S2z, lambda1, lambda2, LALpars, NRTidalv2_V);
break;
case SEOBNRv4T_surrogate:
......@@ -1100,7 +1117,7 @@ int XLALSimInspiralChooseFDWaveformSequence(
/* Call the waveform driver routine */
ret = XLALSimIMRPhenomPFrequencySequence(hptilde, hctilde, frequencies,
chi1_l, chi2_l, chip, thetaJN,
m1, m2, distance, alpha0, phi_aligned, f_ref, IMRPhenomPv1_V, NULL);
m1, m2, distance, alpha0, phi_aligned, f_ref, IMRPhenomPv1_V, NoNRT_V, NULL);
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
cos_2zeta = cos(2.0*zeta_polariz);
sin_2zeta = sin(2.0*zeta_polariz);
......@@ -1132,7 +1149,7 @@ int XLALSimInspiralChooseFDWaveformSequence(
/* Call the waveform driver routine */
ret = XLALSimIMRPhenomPFrequencySequence(hptilde, hctilde, frequencies,
chi1_l, chi2_l, chip, thetaJN,
m1, m2, distance, alpha0, phi_aligned, f_ref, IMRPhenomPv2_V, NULL);
m1, m2, distance, alpha0, phi_aligned, f_ref, IMRPhenomPv2_V, NoNRT_V, NULL);
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
cos_2zeta = cos(2.0*zeta_polariz);
sin_2zeta = sin(2.0*zeta_polariz);
......@@ -1154,7 +1171,7 @@ int XLALSimInspiralChooseFDWaveformSequence(
ABORT_NONZERO_TIDES(LALpars);
ret = XLALSimIMRPhenomDFrequencySequence(hptilde, frequencies,
phiRef, f_ref, m1, m2, S1z, S2z, distance, LALpars);
phiRef, f_ref, m1, m2, S1z, S2z, distance, LALpars, NoNRT_V);
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
/* Produce both polarizations */
*hctilde = XLALCreateCOMPLEX16FrequencySeries("FD hcross",
......@@ -1174,7 +1191,30 @@ int XLALSimInspiralChooseFDWaveformSequence(
ABORT_NONZERO_TRANSVERSE_SPINS(LALpars);
ret = XLALSimIMRPhenomDNRTidalFrequencySequence(hptilde, frequencies,
phiRef, f_ref, distance, m1, m2, S1z, S2z, lambda1, lambda2, LALpars);
phiRef, f_ref, distance, m1, m2, S1z, S2z, lambda1, lambda2, LALpars, NRTidal_V);
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
/* Produce both polarizations */
*hctilde = XLALCreateCOMPLEX16FrequencySeries("FD hcross",
&((*hptilde)->epoch), (*hptilde)->f0, (*hptilde)->deltaF,
&((*hptilde)->sampleUnits), (*hptilde)->data->length);
for(j = 0; j < (*hptilde)->data->length; j++) {
(*hctilde)->data->data[j] = -I*cfac * (*hptilde)->data->data[j];
(*hptilde)->data->data[j] *= pfac;
}
break;
case IMRPhenomD_NRTidalv2:
/* Waveform-specific sanity checks */
if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALpars) )
ABORT_NONDEFAULT_LALDICT_FLAGS(LALpars);
if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
ABORT_NONZERO_TRANSVERSE_SPINS(LALpars);
ret = XLALSimInspiralSetQuadMonParamsFromLambdas(LALpars);
XLAL_CHECK(XLAL_SUCCESS == ret, ret, "Failed to set QuadMon from Lambdas for IMRPhenomD_NRTidalv2");
ret = XLALSimIMRPhenomDNRTidalFrequencySequence(hptilde, frequencies,
phiRef, f_ref, distance, m1, m2, S1z, S2z, lambda1, lambda2, LALpars, NRTidalv2_V);
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
/* Produce both polarizations */
*hctilde = XLALCreateCOMPLEX16FrequencySeries("FD hcross",
......@@ -1203,7 +1243,34 @@ int XLALSimInspiralChooseFDWaveformSequence(
/* Call the waveform driver routine */
ret = XLALSimIMRPhenomPFrequencySequence(hptilde, hctilde, frequencies,
chi1_l, chi2_l, chip, thetaJN,
m1, m2, distance, alpha0, phi_aligned, f_ref, IMRPhenomPv2NRTidal_V, LALpars);
m1, m2, distance, alpha0, phi_aligned, f_ref, IMRPhenomPv2NRTidal_V, NRTidal_V, LALpars);
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
for (UINT4 idx=0;idx<(*hptilde)->data->length;idx++) {
PhPpolp=(*hptilde)->data->data[idx];
PhPpolc=(*hctilde)->data->data[idx];
(*hptilde)->data->data[idx] =cos(2.*zeta_polariz)*PhPpolp+sin(2.*zeta_polariz)*PhPpolc;
(*hctilde)->data->data[idx]=cos(2.*zeta_polariz)*PhPpolc-sin(2.*zeta_polariz)*PhPpolp;
}
break;
case IMRPhenomPv2_NRTidalv2:
/* Waveform-specific sanity checks */
if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALpars) )
ABORT_NONDEFAULT_FRAME_AXIS(LALpars);/* Default is LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L : z-axis along direction of orbital angular momentum. */
if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALpars) )
ABORT_NONDEFAULT_MODES_CHOICE(LALpars);
/* Tranform to model parameters */
if(f_ref==0.0)
f_ref = f_min; /* Default reference frequency is minimum frequency */
XLALSimIMRPhenomPCalculateModelParametersFromSourceFrame(
&chi1_l, &chi2_l, &chip, &thetaJN, &alpha0, &phi_aligned, &zeta_polariz,
m1, m2, f_ref, phiRef, inclination,
S1x, S1y, S1z,
S2x, S2y, S2z, IMRPhenomPv2NRTidal_V);
/* Call the waveform driver routine */
ret = XLALSimIMRPhenomPFrequencySequence(hptilde, hctilde, frequencies,
chi1_l, chi2_l, chip, thetaJN,
m1, m2, distance, alpha0, phi_aligned, f_ref, IMRPhenomPv2NRTidal_V, NRTidalv2_V, LALpars);
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
for (UINT4 idx=0;idx<(*hptilde)->data->length;idx++) {
PhPpolp=(*hptilde)->data->data[idx];
......