22 #include <gsl/gsl_math.h>
24 #include <lal/Sequence.h>
116 const REAL8 distance,
117 LALDict *extraParams,
135 const REAL8 q = (m1 > m2) ? (m1 / m2) : (m2 / m1);
138 XLAL_PRINT_WARNING(
"Warning: The model is not supported for high mass ratio, see MAX_ALLOWED_MASS_RATIO\n");
140 if (chi1 > 1.0 || chi1 < -1.0 || chi2 > 1.0 || chi2 < -1.0)
144 REAL8 fRef = (fRef_in == 0.0) ?
f_min : fRef_in;
156 f_max_prime = (f_max_prime > fCut) ? fCut : f_max_prime;
157 if (f_max_prime <=
f_min)
165 freqs->
data[1] = f_max_prime;
168 distance, extraParams, NRTidal_version);
172 if (f_max_prime <
f_max) {
175 size_t n = (*htilde)->data->length;
178 XLAL_CHECK ( *htilde,
XLAL_ENOMEM,
"Failed to resize waveform COMPLEX16FrequencySeries of length %zu (for internal fCut=%f) to new length %zu (for user-requested f_max=%f).", n, fCut, n_full,
f_max );
209 const REAL8 distance,
210 LALDict *extraParams,
226 const REAL8 q = (m1 > m2) ? (m1 / m2) : (m2 / m1);
229 XLAL_PRINT_WARNING(
"Warning: The model is not supported for high mass ratio, see MAX_ALLOWED_MASS_RATIO\n");
231 if (chi1 > 1.0 || chi1 < -1.0 || chi2 > 1.0 || chi2 < -1.0)
235 REAL8 fRef = (fRef_in == 0.0) ? freqs->
data[0] : fRef_in;
239 distance, extraParams, NRTidal_version);
267 const REAL8 distance,
268 LALDict *extraParams,
276 LALDict *extraParams_in = extraParams;
278 REAL8 dquadmon1_in = 0., dquadmon2_in = 0., lambda1_in = 0, lambda2_in = 0.;
286 REAL8 chi1, chi2, m1, m2, dquadmon1, dquadmon2, lambda1, lambda2;
292 dquadmon1 = dquadmon1_in;
293 dquadmon2 = dquadmon2_in;
294 lambda1 = lambda1_in;
295 lambda2 = lambda2_in;
301 dquadmon1 = dquadmon2_in;
302 dquadmon2 = dquadmon1_in;
303 lambda1 = lambda2_in;
304 lambda2 = lambda1_in;
326 if (eta > 0.25 || eta < 0.0)
342 XLAL_CHECK (
XLALGPSAdd(&ligotimegps_zero, -1. / deltaF),
XLAL_EFUNC,
"Failed to shift coalescence time to t=0, tried to apply shift of -1.0/deltaF with deltaF=%g.", deltaF);
344 XLAL_CHECK ( *htilde,
XLAL_ENOMEM,
"Failed to allocated waveform COMPLEX16FrequencySeries of length %zu for f_max=%f, deltaF=%g.", npts,
f_max, deltaF);
346 size_t iStart = (size_t) (
f_min / deltaF);
347 size_t iStop = (size_t) (
f_max / deltaF);
348 XLAL_CHECK ( (iStop<=npts) && (iStart<=iStop),
XLAL_EDOM,
"minimum freq index %zu and maximum freq index %zu do not fulfill 0<=ind_min<=ind_max<=htilde->data>length=%zu.", iStart, iStop, npts);
353 freqs->
data[
i-iStart] =
i*deltaF;
358 XLAL_CHECK ( *htilde,
XLAL_ENOMEM,
"Failed to allocated waveform COMPLEX16FrequencySeries of length %zu from sequence.", npts);
367 memset((*htilde)->data->data, 0, npts *
sizeof(
COMPLEX16));
375 the model might misbehave here.", finspin);
381 if (extraParams==NULL)
416 const REAL8 MfRef = M_sec * fRef;
420 const REAL8 phifRef =
IMRPhenDPhase(MfRef, pPhi,
pn, &powers_of_fRef, &phi_prefactors, 1.0, 1.0);
423 const REAL8 phi_precalc = 2.*phi0 + phifRef;
432 XLAL_CHECK(
XLAL_SUCCESS == ret, ret,
"Failed to generate tidal amplitude series to construct IMRPhenomD_NRTidalv2 waveform.");
434 #pragma omp parallel for
436 double Mf = M_sec * freqs->
data[
i];
437 double ampT = amp_tidal->
data[
i];
444 XLALPrintError(
"init_useful_powers failed for Mf, status_in_for=%d", status_in_for);
451 phi -= t0*(Mf-MfRef) + phi_precalc;
452 ((*htilde)->data->data)[j] = amp0 * (amp+2*sqrt(
LAL_PI/5.)*ampT) * cexp(-I * phi);
456 #pragma omp parallel for
458 double Mf = M_sec * freqs->
data[
i];
465 XLALPrintError(
"init_useful_powers failed for Mf, status_in_for=%d", status_in_for);
472 phi -= t0*(Mf-MfRef) + phi_precalc;
473 ((*htilde)->data->data)[j] = amp0 * amp * cexp(-I * phi);
487 if (extraParams && !extraParams_in) {
514 REAL8 chi1, chi2, m1, m2;
534 if (eta > 0.25 || eta < 0.0)
542 the model might misbehave here.", finspin);
580 double DPhiMRDval =
DPhiMRD(Mf,
p, 1.0, 1.0) +
p->C2MRD;
607 if (chi1_in > 1.0 || chi1_in < -1.0 || chi2_in > 1.0 || chi2_in < -1.0)
614 REAL8 chi1, chi2, m1, m2;
629 if (fHzSt > fHzPeak){
630 XLAL_PRINT_WARNING(
"Starting frequency = %f Hz is higher IMRPhenomD peak frequency %f Hz. Results may be unreliable.", fHzSt, fHzPeak);
641 if (eta > 0.25 || eta < 0.0)
646 const REAL8 MfSt = M_sec * fHzSt;
653 the model might misbehave here.", finspin);
654 LALDict *extraParams = NULL;
655 if (extraParams == NULL)
691 const REAL8 dphidiff = dphifRD - dphifSt;
694 const REAL8 ChirpTimeSec = dphidiff / 2. /
LAL_PI * M_sec;
716 REAL8 chi1, chi2, m1, m2;
734 if (eta > 0.25 || eta < 0.0)
741 the model might misbehave here.", finspin);
779 &pD, m1, m2, chi1x, chi1y, chi1z,
780 chi2x, chi2y, chi2z, Rholm, Taulm, extraParams);
783 XLALPrintError(
"XLAL Error - IMRPhenomDSetupAmpAndPhaseCoefficients failed\n");
789 #pragma omp parallel for
790 for (
size_t i = ind_min;
i < ind_max;
i++)
798 XLALPrintError(
"init_useful_powers failed for Mf, status_in_for=%d\n", status_in_for);
799 retcode = status_in_for;
846 const REAL8 Mtot = m1 + m2;
847 const REAL8 eta = m1 * m2 / (Mtot * Mtot);
857 the model might misbehave here.",
873 #pragma omp parallel for
874 for (
size_t i = ind_min;
i < ind_max;
i++)
882 XLALPrintError(
"init_useful_powers failed for Mf, status_in_for=%d", status_in_for);
883 retcode = status_in_for;
908 if (extraParams == NULL)
961 LALDict *extraParams_in = extraParams;
971 const REAL8 Mtot = m1 + m2;
972 const REAL8 eta = m1 * m2 / (Mtot * Mtot);
982 the model might misbehave here.",
986 if (extraParams == NULL)
1006 REAL8 testGRcor = 1.0;
1008 pn->v[6] -= (
Subtract3PNSS(m1, m2, Mtot, eta, chi1z, chi2z) *
pn->v[0]) * testGRcor;
1033 pDPreComp->
pn = *
pn;
1034 pDPreComp->
pPhi = *pPhi;
1037 pDPreComp->
pAmp = *pAmp;
1046 if (extraParams && !extraParams_in) {
void XLALDestroyDict(LALDict *dict)
LALDict * XLALCreateDict(void)
int XLALSimNRTunedTidesFDTidalAmplitudeFrequencySeries(const REAL8Sequence *amp_tidal, const REAL8Sequence *fHz, REAL8 m1_SI, REAL8 m2_SI, REAL8 lambda1, REAL8 lambda2)
Function to call amplitude tidal series only; done for convenience to use for PhenomD_NRTidalv2 and S...
static size_t NextPow2(const size_t n)
double XLALSimIMRPhenomDFinalSpin(const REAL8 m1_in, const REAL8 m2_in, const REAL8 chi1_in, const REAL8 chi2_in)
Function to return the final spin (spin of the remnant black hole) as predicted by the IMRPhenomD mod...
double XLALSimIMRPhenomDChirpTime(const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi1_in, const REAL8 chi2_in, const REAL8 fHzSt)
Estimates the length of the time domain IMRPhenomD signal This does NOT taking into account any taper...
UNUSED REAL8 IMRPhenomDPhase_OneFrequency(REAL8 Mf, PhenDAmpAndPhasePreComp pD, REAL8 Rholm, REAL8 Taulm)
Function to return the phenomD phase using the IMRPhenomDSetupAmpAndPhaseCoefficients struct.
int IMRPhenomDPhaseFrequencySequence(REAL8Sequence *phases, REAL8Sequence *freqs, size_t ind_min, size_t ind_max, REAL8 m1, REAL8 m2, REAL8 chi1x, REAL8 chi1y, REAL8 chi1z, REAL8 chi2x, REAL8 chi2y, REAL8 chi2z, REAL8 Rholm, REAL8 Taulm, LALDict *extraParams)
Helper function used in PhenomHM and PhenomPv3HM Returns the phenomD phase, with modified QNM.
static int IMRPhenomDGenerateFD(COMPLEX16FrequencySeries **htilde, const REAL8Sequence *freqs_in, double deltaF, const REAL8 phi0, const REAL8 fRef, const REAL8 m1_in, const REAL8 m2_in, const REAL8 chi1_in, const REAL8 chi2_in, const REAL8 distance, LALDict *extraParams, NRTidal_version_type NRTidal_version)
static double PhenDPhaseDerivFrequencyPoint(double Mf, IMRPhenomDPhaseCoefficients *p, PNPhasingSeries *pn)
Helper function to return the value of the frequency derivative of the Fourier domain phase.
double XLALIMRPhenomDGetPeakFreq(const REAL8 m1_in, const REAL8 m2_in, const REAL8 chi1_in, const REAL8 chi2_in)
Function to return the frequency (in Hz) of the peak of the frequency domain amplitude for the IMRPhe...
int IMRPhenomDAmpFrequencySequence(REAL8Sequence *amps, REAL8Sequence *freqs, size_t ind_min, size_t ind_max, REAL8 m1, REAL8 m2, REAL8 chi1x, REAL8 chi1y, REAL8 chi1z, REAL8 chi2x, REAL8 chi2y, REAL8 chi2z)
Helper function used in PhenomHM and PhenomPv3HM Returns the phenomD amplitude.
REAL8 IMRPhenomDComputet0(REAL8 eta, REAL8 chi1z, REAL8 chi2z, REAL8 finspin, LALDict *extraParams)
computes the time shift as the approximate time of the peak of the 22 mode.
UsefulPowers powers_of_pi
useful powers of LAL_PI, calculated once and kept constant - to be initied with a call to init_useful...
int IMRPhenomDSetupAmpAndPhaseCoefficients(PhenDAmpAndPhasePreComp *pDPreComp, REAL8 m1, REAL8 m2, REAL8 chi1x, REAL8 chi1y, REAL8 chi1z, REAL8 chi2x, REAL8 chi2y, REAL8 chi2z, const REAL8 Rholm, const REAL8 Taulm, LALDict *extraParams)
Function to compute the amplitude and phase coefficients for PhenomD Used to optimise the calls to IM...
#define MAX_ALLOWED_MASS_RATIO
A large mass ratio causes memory over-runs.
#define f_CUT
Dimensionless frequency (Mf) at which define the end of the waveform.
#define MIN_FINAL_SPIN
Minimal final spin value below which the waveform might behave pathological because the ISCO frequenc...
Internal function for IMRPhenomD phenomenological waveform model. See LALSimIMRPhenom....
static void ComputeIMRPhenDPhaseConnectionCoefficients(IMRPhenomDPhaseCoefficients *p, PNPhasingSeries *pn, PhiInsPrefactors *prefactors, double Rholm, double Taulm)
This function aligns the three phase parts (inspiral, intermediate and merger-rindown) such that they...
static double IMRPhenDPhase(double f, IMRPhenomDPhaseCoefficients *p, PNPhasingSeries *pn, UsefulPowers *powers_of_f, PhiInsPrefactors *prefactors, double Rholm, double Taulm)
This function computes the IMR phase given phenom coefficients.
static int init_phi_ins_prefactors(PhiInsPrefactors *prefactors, IMRPhenomDPhaseCoefficients *p, PNPhasingSeries *pn)
static double FinalSpin0815(double eta, double chi1, double chi2)
Wrapper function for FinalSpin0815_s.
static void ComputeIMRPhenomDPhaseCoefficients(IMRPhenomDPhaseCoefficients *p, double eta, double chi1, double chi2, double finspin, LALDict *extraParams)
A struct containing all the parameters that need to be calculated to compute the phenomenological pha...
static double DPhiIntAnsatz(double Mf, IMRPhenomDPhaseCoefficients *p)
First frequency derivative of PhiIntAnsatz (this time with 1.
static void ComputeIMRPhenomDAmplitudeCoefficients(IMRPhenomDAmplitudeCoefficients *p, double eta, double chi1, double chi2, double finspin)
A struct containing all the parameters that need to be calculated to compute the phenomenological amp...
static double Subtract3PNSS(double m1, double m2, double M, double eta, double chi1, double chi2)
Subtract 3PN spin-spin term below as this is in LAL's TaylorF2 implementation (LALSimInspiralPNCoeffi...
static double IMRPhenDAmplitude(double f, IMRPhenomDAmplitudeCoefficients *p, UsefulPowers *powers_of_f, AmpInsPrefactors *prefactors)
This function computes the IMR amplitude given phenom coefficients.
static double DPhiInsAnsatzInt(double Mf, IMRPhenomDPhaseCoefficients *p, PNPhasingSeries *pn)
First frequency derivative of PhiInsAnsatzInt.
static bool StepFunc_boolean(const double t, const double t1)
**
static int init_amp_ins_prefactors(AmpInsPrefactors *prefactors, IMRPhenomDAmplitudeCoefficients *p)
static double DPhiMRD(double f, IMRPhenomDPhaseCoefficients *p, double Rholm, double Taulm)
First frequency derivative of PhiMRDAnsatzInt Rholm was added when IMRPhenomHM (high mode) was added.
static int init_useful_powers(UsefulPowers *p, REAL8 number)
int PhenomInternal_PrecessingSpinEnforcePrimaryIsm1(REAL8 *m1, REAL8 *m2, REAL8 *chi1x, REAL8 *chi1y, REAL8 *chi1z, REAL8 *chi2x, REAL8 *chi2y, REAL8 *chi2z)
Given m1 with spins (chi1x, chi1y, chi1z) and m2 with spins (chi2x,chi2y,chi2z).
void PhenomInternal_nudge(REAL8 *x, REAL8 X, REAL8 epsilon)
If x and X are approximately equal to relative accuracy epsilon then set x = X.
REAL8 XLALSimPhenomUtilsPhenomPv2FinalSpin(const REAL8 m1, const REAL8 m2, const REAL8 chi1_l, const REAL8 chi2_l, const REAL8 chip)
Wrapper for final-spin formula based on:
REAL8 XLALSimPhenomUtilsChiP(const REAL8 m1, const REAL8 m2, const REAL8 s1x, const REAL8 s1y, const REAL8 s2x, const REAL8 s2y)
Function to compute the effective precession parameter chi_p (1408.1810)
COMPLEX16FrequencySeries * XLALResizeCOMPLEX16FrequencySeries(COMPLEX16FrequencySeries *series, int first, size_t length)
COMPLEX16FrequencySeries * XLALCreateCOMPLEX16FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
void * XLALMalloc(size_t n)
@ NRTidalv2_V
version NRTidalv2: https://arxiv.org/abs/1905.06011
int XLALSimIMRPhenomDGenerateFD(COMPLEX16FrequencySeries **htilde, const REAL8 phi0, const REAL8 fRef_in, const REAL8 deltaF, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi1, const REAL8 chi2, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, LALDict *extraParams, NRTidal_version_type NRTidal_version)
Driver routine to compute the spin-aligned, inspiral-merger-ringdown phenomenological waveform IMRPhe...
int XLALSimIMRPhenomDFrequencySequence(COMPLEX16FrequencySeries **htilde, const REAL8Sequence *freqs, const REAL8 phi0, const REAL8 fRef_in, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi1, const REAL8 chi2, const REAL8 distance, LALDict *extraParams, NRTidal_version_type NRTidal_version)
Compute waveform in LAL format at specified frequencies for the IMRPhenomD model.
@ LAL_SIM_INSPIRAL_SPIN_ORDER_35PN
@ LAL_SIM_INSPIRAL_SPIN_ORDER_ALL
int XLALSimInspiralTaylorF2AlignedPhasing(PNPhasingSeries **pfa, const REAL8 m1, const REAL8 m2, const REAL8 chi1, const REAL8 chi2, LALDict *extraPars)
Returns structure containing TaylorF2 phasing coefficients for given physical parameters.
void XLALDestroyREAL8Sequence(REAL8Sequence *sequence)
REAL8Sequence * XLALCreateREAL8Sequence(size_t length)
const LALUnit lalStrainUnit
const LALUnit lalSecondUnit
LALUnit * XLALUnitMultiply(LALUnit *output, const LALUnit *unit1, const LALUnit *unit2)
#define XLAL_CHECK(assertion,...)
#define XLAL_PRINT_WARNING(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
used to cache the recurring (frequency-independent) prefactors of AmpInsAnsatz.
Structure holding all coefficients for the amplitude.
Structure holding all coefficients for the phase.
A struct the store the amplitude and phase structs for phenomD.
AmpInsPrefactors amp_prefactors
IMRPhenomDPhaseCoefficients pPhi
PhiInsPrefactors phi_prefactors
IMRPhenomDAmplitudeCoefficients pAmp
used to cache the recurring (frequency-independent) prefactors of PhiInsAnsatzInt.
useful powers in GW waveforms: 1/6, 1/3, 2/3, 4/3, 5/3, 2, 7/3, 8/3, -1, -1/6, -7/6,...