30 #include <lal/LALStdlib.h>
31 #include <lal/LALConstants.h>
32 #include <lal/LALDatatypes.h>
33 #include <lal/XLALError.h>
41 #ifndef PHENOMXHMDEBUG
63 if (Mf <= Mf_beta_lower)
73 if (Mf >= Mf_beta_upper)
116 double msa_window = 1-pnr_window;
122 if (Mf <= Mf_beta_lower)
132 if (Mf >= Mf_beta_upper)
195 const double omega =
LAL_PI * Mf;
196 const double omega_cbrt = cbrt(omega);
206 REAL8 L =
XLALSimIMRPhenomXLPNAnsatz(omega_cbrt, pWF->
eta / omega_cbrt, pPrec->
L0, pPrec->
L1, pPrec->
L2, pPrec->
L3, pPrec->
L4, pPrec->
L5, pPrec->
L6, pPrec->
L7, pPrec->
L8, pPrec->
L8L);
215 beta = acos(copysign(1.0, L + pPrec->
SL) / sqrt(1.0 + s2));
225 vector vangles = {0., 0., 0.};
230 REAL8 beta_full = acos(vangles.
z);
240 REAL8 beta_SingleSpin = acos(vangles_SingleSpin.
z);
243 if (Mf <= Mf_beta_lower)
246 REAL8 oscillations = beta_full - beta_SingleSpin;
247 REAL8 envelope = cos(2.0 *
LAL_PI * Mf / (4.0 * Mf_beta_lower)) * cos(2.0 *
LAL_PI * Mf / (4.0 * Mf_beta_lower));
248 beta = (beta_SingleSpin + oscillations * envelope);
252 beta = beta_SingleSpin;
264 XLAL_ERROR(
XLAL_EINVAL,
"Error: IMRPhenomXPrecessionVersion not recognized in IMRPhenomX_PNR_GetPNBetaAtFreq.\n");
287 const double omega =
LAL_PI * Mf;
288 const double omega_cbrt = cbrt(omega);
299 const REAL8 L =
XLALSimIMRPhenomXLPNAnsatz(omega_cbrt, pWF->
eta / omega_cbrt, pPrec->
L0, pPrec->
L1, pPrec->
L2, pPrec->
L3, pPrec->
L4, pPrec->
L5, pPrec->
L6, pPrec->
L7, pPrec->
L8, pPrec->
L8L);
307 beta = acos(copysign(1.0, L + pPrec->
SL) / sqrt(1.0 + s2));
317 vector vangles = {0., 0., 0.};
322 beta = acos(vangles.
z);
328 XLAL_ERROR(
XLAL_EINVAL,
"Error: IMRPhenomXPrecessionVersion not recognized in IMRPhenomX_PNR_GetPNBetaAtFreq_fulltwospin.\n");
359 REAL8 chi_parr = (m1 + m2) * chi_eff / m1;
364 REAL8 L_3PN =
M*
M*
XLALSimIMRPhenomXLPNAnsatz(v, L_norm, pPrec->
L0, pPrec->
L1, pPrec->
L2, pPrec->
L3, pPrec->
L4, pPrec->
L5, pPrec->
L6, pPrec->
L7, pPrec->
L8, pPrec->
L8L);
368 REAL8 costheta_temp = chi_parr / chi_temp;
392 REAL8 S_perp = (J0 - (L0 - L)) * sin(
beta);
393 REAL8 chi_p = S_perp / (m1 * m1);
396 REAL8 chi = sqrt(chi_parr * chi_parr + chi_p * chi_p);
421 REAL8 nu = (m1 * m2) / (
M *
M);
429 REAL8 v = pow(w_orb, 1.0 / 3.0);
434 REAL8 cos_iota = cos(iota);
435 REAL8 sin_iota = sin(iota);
437 REAL8 cos_half_iota = cos(iota / 2.0);
438 REAL8 sin_half_iota = sin(iota / 2.0);
444 REAL8 N0 = 84.0 * sin_iota;
445 REAL8 N2 = 2.0 * (55.0 * nu - 107.0) * sin_iota;
446 REAL8 N3 = -7.0 * (5.0 * nu + 6.0 *
delta + 6.0) * chi * (2.0 * cos_iota - 1.0) * sin_theta + 56.0 * (3.0 *
LAL_PI - (1.0 +
delta - nu) * chi * cos_theta) * sin_iota;
448 REAL8 N = (N0 + N2 * v2 + N3 * v3) / cos_half_iota;
451 REAL8 D0 = 84.0 * cos_half_iota;
452 REAL8 D2 = 2.0 * (55.0 * nu - 107.0) * cos_half_iota;
453 REAL8 D3 = 14.0 * 4.0 * (3.0 *
LAL_PI + (nu - 1.0 -
delta) * chi * cos_theta) * cos_half_iota + 14.0 * (6.0 + 6.0 *
delta + 5.0 * nu) * chi * sin_theta * sin_half_iota;
524 REAL8 D = beta1 * beta1 * Mf;
525 REAL8 N = -2.0 * beta1 * (beta2 - beta1) + (beta1 * dbeta2 - beta2 * dbeta1) * Mf;
539 REAL8 D = (beta1 * Mf) * (beta1 * Mf);
540 REAL8 N = beta1 * (beta2 - beta1) - (beta1 * dbeta2 - beta2 * dbeta1) * Mf;
554 return 1.0 +
b1 * Mf +
b2 * Mf * Mf;
580 return B0 + (B1 + B2 * Mf + B3 * Mf * Mf) / (1.0 + B4 * (Mf + B5) * (Mf + B5));
598 return ((2.0 * B3 * B4 * B5 - B2 * B4) * Mf * Mf + (2.0 * B3 - 2.0 * B1 * B4 + 2.0 * B3 * B4 * B5 * B5) * Mf + (B2 - 2.0 * B1 * B4 * B5 + B2 * B4 * B5 * B5)) / pow((1.0 + B4 * (Mf + B5) * (Mf + B5)), 2);
616 REAL8 a = B2 * B4 * B4 - 2.0 * B3 * B4 * B4 * B5;
617 REAL8 b = -3.0 * B3 * B4 + 3.0 * B1 * B4 * B4 - 3.0 * B3 * B4 * B4 * B5 * B5;
618 REAL8 c = -3.0 * B2 * B4 + 6.0 * B1 * B4 * B4 * B5 - 3.0 * B2 * B4 * B4 * B5 * B5;
619 REAL8 d = B3 - B1 * B4 - 2.0 * B2 * B4 * B5 + 2.0 * B3 * B4 * B5 * B5 + 3.0 * B1 * B4 * B4 * B5 * B5 - 2.0 * B2 * B4 * B4 * B5 * B5 * B5 + B3 * B4 * B4 * B5 * B5 * B5 * B5;
621 return 2.0 * (
a * Mf * Mf * Mf + b * Mf * Mf +
c * Mf +
d) / pow((1.0 + B4 * (B5 + Mf) * (B5 + Mf)), 3);
639 REAL8 a = -B2 * B4 * B4 + 2.0 * B3 * B4 * B4 * B5;
640 REAL8 b = 4.0 * B3 * B4 - 4.0 * B1 * B4 * B4 + 4.0 * B3 * B4 * B4 * B5 * B5;
641 REAL8 c = 6.0 * B2 * B4 - 12.0 * B1 * B4 * B4 * B5 + 6.0 * B2 * B4 * B4 * B5 * B5;
642 REAL8 d = -4.0 * B3 + 4.0 * B1 * B4 + 8.0 * B2 * B4 * B5 - 8.0 * B3 * B4 * B5 * B5 - 12.0 * B1 * B4 * B4 * B5 * B5 + 8.0 * B2 * B4 * B4 * B5 * B5 * B5 - 4.0 * B3 * B4 * B4 * B5 * B5 * B5 * B5;
643 REAL8 e = -B2 - 2.0 * B3 * B5 + 4.0 * B1 * B4 * B5 + 2.0 * B2 * B4 * B5 * B5 - 4.0 * B3 * B4 * B5 * B5 * B5 - 4.0 * B1 * B4 * B4 * B5 * B5 * B5 + 3.0 * B2 * B4 * B4 * B5 * B5 * B5 * B5 - 2.0 * B3 * B4 * B4 * B5 * B5 * B5 * B5 * B5;
645 return 6.0 * B4 * (
a * Mf * Mf * Mf * Mf + b * Mf * Mf * Mf +
c * Mf * Mf +
d * Mf +
e) / pow((1.0 + B4 * (B5 + Mf) * (B5 + Mf)), 4);
677 REAL8 ddbeta_Mf_plus;
685 REAL8 dbeta_inflection;
699 root_term = B4 * (B2 - 2.0 * B3 * B5) * (B2 - 2.0 * B1 * B4 * B5 + B2 * B4 * B5 * B5) + (B3 - B1 * B4 + B3 * B4 * B5 * B5) * (B3 - B1 * B4 + B3 * B4 * B5 * B5);
701 Mf_plus = ((B3 - B1 * B4 + B3 * B4 * B5 * B5) + sqrt(root_term)) / (B4 * (B2 - 2.0 * B3 * B5));
702 Mf_minus = ((B3 - B1 * B4 + B3 * B4 * B5 * B5) - sqrt(root_term)) / (B4 * (B2 - 2.0 * B3 * B5));
707 if (ddbeta_Mf_plus > 0.0)
710 Mf_at_minimum = Mf_plus;
711 Mf_at_maximum = Mf_minus;
716 Mf_at_minimum = Mf_minus;
717 Mf_at_maximum = Mf_plus;
727 if (dbeta_inflection > 0.0)
736 chosen_dbeta = sign * ((dbeta_inflection / 100.0) * (dbeta_inflection / 100.0)) * 25.0;
739 d1 = 1.0 / dddbeta * (-ddbeta + sqrt(ddbeta * ddbeta + 2.0 * dddbeta * chosen_dbeta));
740 d2 = 1.0 / dddbeta * (-ddbeta - sqrt(ddbeta * ddbeta + 2.0 * dddbeta * chosen_dbeta));
751 delta = (d1 < d2) ? d1 : d2;
766 if (Mf_inflection >= 0.06)
768 Mf_low = Mf_inflection - 0.03;
772 Mf_low = 3.0 * Mf_inflection / 5.0;
781 if ((Mf_at_minimum > Mf_at_maximum) || (Mf_inflection > Mf_at_maximum))
785 if (Mf_at_maximum >= Mf_low)
787 Mf_IM = Mf_at_maximum +
delta;
797 if (Mf_at_minimum > 0.06)
799 Mf_IM = Mf_at_minimum - 0.03;
803 Mf_IM = 3.0 * Mf_at_minimum / 5.0;
811 if (Mf_at_minimum > Mf_inflection)
813 Mf_MR = Mf_at_minimum;
821 if ((Mf_IM < 0.0) || isnan(Mf_IM))
861 a = 2 * (B2 * B4 * B4 - 2 * B3 * B4 * B4 * B5);
862 b = 6 * (-B3 * B4 + B1 * B4 * B4 - B3 * B4 * B4 * B5 * B5);
863 c = 6 * (-B2 * B4 + 2 * B1 * B4 * B4 * B5 - B2 * B4 * B4 * B5 * B5);
864 d = 2 * (B3 - B1 * B4 - 2 * B2 * B4 * B5 + 2 * B3 * B4 * B5 * B5 + 3 * B1 * B4 * B4 * B5 * B5 - 2 * B2 * B4 * B4 * B5 * B5 * B5 + B3 * B4 * B4 * B5 * B5 * B5 * B5);
867 r = (3 *
a *
c - b * b) / (3 *
a *
a);
868 s = (2 * b * b * b - 9 *
a * b *
c + 27 *
a *
a *
d) / (27 *
a *
a *
a);
870 phi = acos(((3 *
s) / (2 *
r)) * csqrt(-3 /
r));
872 for (
int n = 0; n < 3; n++)
874 f[n] = 2 * csqrt(-
r / 3) * cos((phi - 2 * n *
LAL_PI) / 3) - b / (3 *
a);
903 b = 6 * (-B3 * B4 + B1 * B4 * B4 - B3 * B4 * B4 * B5 * B5);
904 c = 6 * (-B2 * B4 + 2 * B1 * B4 * B4 * B5 - B2 * B4 * B4 * B5 * B5);
905 d = 2 * (B3 - B1 * B4 - 2 * B2 * B4 * B5 + 2 * B3 * B4 * B5 * B5 + 3 * B1 * B4 * B4 * B5 * B5 - 2 * B2 * B4 * B4 * B5 * B5 * B5 + B3 * B4 * B4 * B5 * B5 * B5 * B5);
907 f[0] = (1 / (2 * b)) * (-
c + csqrt(
c *
c - 4 * b *
d));
908 f[1] = (1 / (2 * b)) * (-
c - csqrt(
c *
c - 4 * b *
d));
940 REAL8 f_inflection = 0.0;
943 a = 2.0 * (B2 * B4 * B4 - 2.0 * B3 * B4 * B4 * B5);
944 b = 6.0 * (-B3 * B4 + B1 * B4 * B4 - B3 * B4 * B4 * B5 * B5);
945 c = 6.0 * (-B2 * B4 + 2.0 * B1 * B4 * B4 * B5 - B2 * B4 * B4 * B5 * B5);
946 d = 2.0 * (B3 - B1 * B4 - 2.0 * B2 * B4 * B5 + 2.0 * B3 * B4 * B5 * B5 + 3.0 * B1 * B4 * B4 * B5 * B5 - 2.0 * B2 * B4 * B4 * B5 * B5 * B5 + B3 * B4 * B4 * B5 * B5 * B5 * B5);
955 f_inflection = -
d /
c;
965 for (
i = 0;
i < 2;
i++)
970 f_inflection = f_inf[
i];
987 for (
i = 0;
i < 3;
i++)
989 f_IM = cimag(f_inf[
i]);
993 f_temp = creal(f_inf[
i]);
999 f_inflection = f_temp;
1007 f_inflection = creal(f_inf[1]);
1012 if (b / (3.0 *
a) > B5 / 2.0 - (2141.0 / 90988.0))
1014 f_inflection = creal(f_inf[0]);
1018 f_inflection = creal(f_inf[2]);
1024 return f_inflection;
1071 REAL8 derivative_beta_lower = (b3 -
b1) / (2.0 * dMf);
1072 REAL8 derivative_beta_upper = (b6 - b4) / (2.0 * dMf);
1078 beta_rescale_1 = isnan(beta_rescale_1) ? 0.0 : beta_rescale_1;
1079 beta_rescale_2 = isnan(beta_rescale_2) ? 0.0 : beta_rescale_2;
1102 REAL8 window_border = 0.01;
1109 REAL8 PI_by_2 = 1.570796326794897;
1110 REAL8 PI_by_2_1mp = 1.569378278348018;
1111 REAL8 PI_by_2_1oq = 7.308338225719002e97;
1112 REAL8 sign = (
beta < PI_by_2) ? -1.0 : 1.0;
1114 return sign * PI_by_2_1mp * pow(atan2(pow(
beta - PI_by_2, 1.0 /
p), PI_by_2_1oq),
p) + PI_by_2;
#define MAX_TOL_ATAN
Tolerance used below which numbers are treated as zero for the calculation of atan2.
static double double delta
REAL8 IMRPhenomX_PNR_MR_dbeta_expression(REAL8 Mf, const IMRPhenomX_PNR_beta_parameters *betaParams)
expression for first derivative of beta in merger-ringdown regime
UINT4 IMRPhenomX_PNR_AttachMRBeta(const IMRPhenomX_PNR_beta_parameters *betaParams)
Determine whether to attach the MR contributions to beta.
COMPLEX16 * IMRPhenomX_PNR_three_inflection_points(const IMRPhenomX_PNR_beta_parameters *betaParams)
Compute the three roots of a depressed cubic given by Eq.
REAL8 IMRPhenomX_PNR_MR_ddbeta_expression(REAL8 Mf, const IMRPhenomX_PNR_beta_parameters *betaParams)
expression for second derivative of beta in merger-ringdown regime
REAL8 IMRPhenomX_PNR_GeneratePNRBetaNoMR(REAL8 Mf, IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec)
This function evaluates only the rescaled inspiral beta given in Eq.
REAL8 IMRPhenomX_PNR_MR_dddbeta_expression(REAL8 Mf, const IMRPhenomX_PNR_beta_parameters *betaParams)
expression for third derivative of beta in merger-ringdown regime
REAL8 IMRPhenomX_PNR_GenerateRingdownPNRBeta(IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec)
We evaluate beta at the final Mf_beta_upper connection frequency to approximate the final value of be...
REAL8 IMRPhenomX_PNR_GetPNBetaAtFreq(REAL8 Mf, const IMRPhenomX_PNR_beta_parameters *betaParams, IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec, IMRPhenomXWaveformStruct *pWF_SingleSpin, IMRPhenomXPrecessionStruct *pPrec_SingleSpin)
A wrapper to produce either the NNLO or MSA beta depending on the IMRPhenomXPrecVersion.
REAL8 IMRPhenomX_PNR_arctan_window(REAL8 beta)
Utility function to compute the arctan windowing described in Eq.
int IMRPhenomX_PNR_precompute_beta_coefficients(IMRPhenomX_PNR_beta_parameters *betaParams, IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec)
This function evaluates the Ansatz coefficients of beta outlined in Eq.
REAL8 IMRPhenomX_PNR_PNWaveformBetaWrapper(REAL8 Mf, REAL8 pn_beta, IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec)
A wrapper to generate the "waveform" PN beta from Eq.
REAL8 IMRPhenomX_PNR_MR_beta_expression(REAL8 Mf, const IMRPhenomX_PNR_beta_parameters *betaParams)
These four functions produce the MR Ansatz of beta described in Sec.
REAL8 IMRPhenomX_PNR_GetPNBetaAtFreq_fulltwospin(REAL8 Mf, IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec)
A wrapper to produce either the NNLO or MSA beta depending on the IMRPhenomXPrecVersion.
REAL8 IMRPhenomX_PNR_beta_rescaling_2(REAL8 Mf, REAL8 beta1, REAL8 beta2, REAL8 dbeta1, REAL8 dbeta2)
REAL8 IMRPhenomX_PNR_GenerateMergedPNRBetaAtMf(REAL8 Mf, const IMRPhenomX_PNR_beta_parameters *betaParams, IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec, IMRPhenomXWaveformStruct *pWF_SingleSpin, IMRPhenomXPrecessionStruct *pPrec_SingleSpin)
This function generates beta with the tuned angles and PN expressions blended during merger-ringdown.
REAL8 IMRPhenomX_PNR_chi_calc(REAL8 m1, REAL8 L, REAL8 J0, REAL8 L0, REAL8 chi_parr, REAL8 beta)
The magnitude of the effective total spin is computed from the total and orbital angular momenta,...
int IMRPhenomX_PNR_BetaConnectionFrequencies(IMRPhenomX_PNR_beta_parameters *betaParams)
Here we work through the construction of the connection frequency for beta, outlined in Sec.
int IMRPhenomX_PNR_beta_connection_parameters(IMRPhenomX_PNR_beta_parameters *betaParams, IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec, IMRPhenomXWaveformStruct *pWF_SingleSpin, IMRPhenomXPrecessionStruct *pPrec_SingleSpin)
This function combines several functions together to compute the connection frequencies and the inspi...
REAL8 IMRPhenomX_PNR_rescale_beta_expression(REAL8 Mf, const IMRPhenomX_PNR_beta_parameters *betaParams)
REAL8 IMRPhenomX_PNR_GeneratePNRBetaAtMf(REAL8 Mf, const IMRPhenomX_PNR_beta_parameters *betaParams, IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec, IMRPhenomXWaveformStruct *pWF_SingleSpin, IMRPhenomXPrecessionStruct *pPrec_SingleSpin)
This function evaluates Eqs.
REAL8 IMRPhenomX_PNR_beta_rescaling_1(REAL8 Mf, REAL8 beta1, REAL8 beta2, REAL8 dbeta1, REAL8 dbeta2)
These three functions produce the inspiral rescaling of beta described in Sec.
COMPLEX16 * IMRPhenomX_PNR_two_inflection_points(const IMRPhenomX_PNR_beta_parameters *betaParams)
Compute the two roots of a depressed cubic given by Eq.
REAL8 IMRPhenomX_PNR_PNWaveformBeta(REAL8 Mf, REAL8 iota, REAL8 m1, REAL8 m2, REAL8 chi, REAL8 costheta)
The "waveform" PN beta from Eq.
REAL8 IMRPhenomX_PNR_single_inflection_point(const IMRPhenomX_PNR_beta_parameters *betaParams)
Compute the roots of a depressed cubic given by Eq.
REAL8 IMRPhenomX_PNR_beta_B4_coefficient(REAL8 eta, REAL8 chi, REAL8 costheta)
REAL8 IMRPhenomX_PNR_beta_B1_coefficient(REAL8 eta, REAL8 chi, REAL8 costheta)
REAL8 IMRPhenomX_PNR_beta_B5_coefficient(REAL8 eta, REAL8 chi, REAL8 costheta)
REAL8 IMRPhenomX_PNR_beta_B3_coefficient(REAL8 eta, REAL8 chi, REAL8 costheta)
REAL8 IMRPhenomX_PNR_beta_B2_coefficient(REAL8 eta, REAL8 chi, REAL8 costheta)
REAL8 IMRPhenomX_PNR_beta_Bf_coefficient(REAL8 eta, REAL8 chi, REAL8 costheta)
REAL8 IMRPhenomX_PNR_beta_B0_coefficient(REAL8 eta, REAL8 chi, REAL8 costheta)
expressions for each of the co-efficients that appear in the merger-ringdown expression for beta
REAL8 IMRPhenomX_PNR_AnglesWindow(IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec)
INT4 IMRPhenomX_PNR_CheckTwoSpin(IMRPhenomXPrecessionStruct *pPrec)
This function quickly checks to see if we expect two-spin effects to be present in the inspiral prece...
REAL8 XLALSimIMRPhenomXLPNAnsatz(REAL8 v, REAL8 LNorm, REAL8 L0, REAL8 L1, REAL8 L2, REAL8 L3, REAL8 L4, REAL8 L5, REAL8 L6, REAL8 L7, REAL8 L8, REAL8 L8L)
This is a convenient wrapper function for PN orbital angular momentum.
vector IMRPhenomX_Return_phi_zeta_costhetaL_MSA(const double v, IMRPhenomXWaveformStruct *pWF, IMRPhenomXPrecessionStruct *pPrec)
Wrapper to generate , and at a given frequency.
REAL8 XLALSimIMRPhenomXatan2tol(REAL8 a, REAL8 b, REAL8 tol)
INT4 D0(REAL8 *f, REAL8 dx, INT4 n, REAL8 *df)
XLALWignerdMatrix, XLALWignerDMatrix in <lal/SphericalHarmonics.h>
INT4 D2(REAL8 *f, REAL8 dx, INT4 n, REAL8 *d2f)
static vector d(const double L_norm, const double J_norm, const vector roots)
Internal function that returns the coefficients "d_0", "d_2" and "d_4" from 1703.03967 corresponding ...
static double beta(const double a, const double b, const sysq *system)
Internal function that computes the spin-orbit couplings.
#define XLAL_CHECK(assertion,...)
REAL8 derivative_beta_lower
REAL8 derivative_beta_upper
REAL8 costheta_singleSpin
Polar angle of effective single spin, Eq.
INT4 IMRPhenomXPrecVersion
Flag to set version of Euler angles used.
REAL8 chi_singleSpin
Magnitude of effective single spin used for tapering two-spin angles, Eq.
REAL8 costheta_final_singleSpin
Polar angle of approximate final spin, see technical document FIXME: add reference.