20 #ifndef _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORM_V3OPT_C
21 #define _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORM_V3OPT_C
22 #include <gsl/gsl_deriv.h>
24 #include <lal/LALSimInspiral.h>
25 #include <lal/LALSimIMR.h>
118 for (
INT4 l = 2;
l <= lMax;
l++)
130 REAL8 UNUSED
r ,
pp, Omega, v2, k, hathatk, eulerlogxabs;
132 REAL8 UNUSED rcrossp_x, rcrossp_y, rcrossp_z;
133 REAL8 Slm , rholm, rholmPwrl;
135 REAL8 hathatksq4, hathatk4pi, Tlmprefac, Tlmprodfac;
152 eta =
params->eobParams->eta;
156 XLALPrintError(
"XLAL Error - %s: Eta seems to be > 0.25 - this isn't allowed!\n", __func__);
168 pp = values->data[3];
170 rcrossp_x = cartvalues->data[1] * cartvalues->data[5] - cartvalues->data[2] * cartvalues->data[4];
171 rcrossp_y = cartvalues->data[2] * cartvalues->data[3] - cartvalues->data[0] * cartvalues->data[5];
172 rcrossp_z = cartvalues->data[0] * cartvalues->data[4] - cartvalues->data[1] * cartvalues->data[3];
185 vPhi =
r * cbrt(vPhi);
188 XLAL_PRINT_INFO(
"In XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole, getting rW = %.12e\n",
203 XLAL_PRINT_INFO(
"Calculating hNewton, with v = %.12e, vPhi = %.12e, r = %.12e, Phi = %.12e, l = %d, m = %d\n",
212 if (((
l +
m) % 2) == 0) {
219 XLAL_PRINT_INFO(
"In XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform: Hreal = %e, Slm = %e, eta = %e\n",
Hreal, Slm, eta);
228 hathatksq4 = 4. * hathatk * hathatk;
229 hathatk4pi = 4. *
LAL_PI * hathatk;
237 if (
status != GSL_SUCCESS) {
238 XLALPrintError(
"XLAL Error - %s: Error in GSL function\n", __func__);
247 Tlmprefac = sqrt(hathatk4pi / (1. - exp(-hathatk4pi))) /
z2.val;
250 for (Tlmprodfac = 1.,
i = 1;
i <=
l;
i++) {
251 Tlmprodfac *= (hathatksq4 + (
REAL8)
i *
i);
254 Tlm = Tlmprefac * sqrt(Tlmprodfac);
282 XLAL_PRINT_INFO(
"PK:: rho22v2 = %.12e, rho22v3 = %.12e, rho22v4 = %.12e,\n rho22v5 = %.16e, rho22v6 = %.16e, rho22v6LOG = %.16e, \n rho22v7 = %.12e, rho22v8 = %.16e, rho22v8LOG = %.16e, \n rho22v10 = %.16e, rho22v10LOG = %.16e\n, rho22v6 = %.12e, rho22v8 = %.12e, rho22v10 = %.12e\n",
293 rholm = 1. + v * (hCoeffs->
rho21v1
299 auxflm = v * hCoeffs->
f21v1 + v2 * v * hCoeffs->
f21v3;
313 auxflm = v * v2 * hCoeffs->
f33v3;
316 rholm = 1. + v * (hCoeffs->
rho32v
325 auxflm = v * v2 * hCoeffs->
f31v3;
336 rholm = 1. + v2 * (hCoeffs->
rho44v2
339 + hCoeffs->
rho44v6l * eulerlogxabs) * v))));
342 rholm = 1. + v * (hCoeffs->
rho43v
346 auxflm = v * hCoeffs->
f43v;
349 rholm = 1. + v2 * (hCoeffs->
rho42v2
354 rholm = 1. + v * (hCoeffs->
rho41v
358 auxflm = v * hCoeffs->
f41v;
368 rholm = 1. + v2 * (hCoeffs->
rho55v2
377 rholm = 1. + v2 * (hCoeffs->
rho53v2
385 rholm = 1. + v2 * (hCoeffs->
rho51v2
427 rholm = 1. + hCoeffs->
rho76v2 * v2;
433 rholm = 1. + hCoeffs->
rho74v2 * v2;
439 rholm = 1. + hCoeffs->
rho72v2 * v2;
452 rholm = 1. + hCoeffs->
rho88v2 * v2;
455 rholm = 1. + hCoeffs->
rho87v2 * v2;
458 rholm = 1. + hCoeffs->
rho86v2 * v2;
461 rholm = 1. + hCoeffs->
rho85v2 * v2;
464 rholm = 1. + hCoeffs->
rho84v2 * v2;
467 rholm = 1. + hCoeffs->
rho83v2 * v2;
470 rholm = 1. + hCoeffs->
rho82v2 * v2;
473 rholm = 1. + hCoeffs->
rho81v2 * v2;
502 if (eta == 0.25 &&
m % 2) {
508 if (
r > 0.0 && debugPK) {
509 XLAL_PRINT_INFO(
"YP::dynamics variables in waveform: %i, %i, r = %.12e, v = %.12e\n",
l,
m,
r, v);
510 XLAL_PRINT_INFO(
"rholm^l = %.16e, Tlm = %.16e + i %.16e, \nSlm = %.16e, hNewton = %.16e + i %.16e, delta = %.16e\n", rholmPwrl, Tlm, 0.0, Slm, creal(hNewton), cimag(hNewton), 0.0);
513 *hlm = Tlm * Slm * rholmPwrl;
549 REAL8 UNUSED
r ,
pp, Omega, v2, Omegav2, vh, vh3, k, hathatk, eulerlogxabs;
551 REAL8 UNUSED rcrossp_x, rcrossp_y, rcrossp_z;
552 REAL8 Slm , deltalm, rholm, rholmPwrl;
554 REAL8 hathatksq4, hathatk4pi, Tlmprefac, Tlmprodfac;
572 eta =
params->eobParams->eta;
576 XLALPrintError(
"XLAL Error - %s: Eta seems to be > 0.25 - this isn't allowed!\n", __func__);
586 pp = values->data[3];
588 rcrossp_x = cartvalues->data[1] * cartvalues->data[5] - cartvalues->data[2] * cartvalues->data[4];
589 rcrossp_y = cartvalues->data[2] * cartvalues->data[3] - cartvalues->data[0] * cartvalues->data[5];
590 rcrossp_z = cartvalues->data[0] * cartvalues->data[4] - cartvalues->data[1] * cartvalues->data[3];
603 if (
params->alignedSpins) {
610 vPhi =
r * cbrt(vPhi);
613 XLAL_PRINT_INFO(
"In XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole, getting rW = %.12e\n",
624 vPhi =
r * cbrt(vPhi);
627 XLAL_PRINT_INFO(
"In XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole, getting rW = %.12e\n",
638 XLAL_PRINT_INFO(
"\nValues inside XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform:\n");
639 for (
i = 0;
i < 11;
i++)
642 XLAL_PRINT_INFO(
"Calculating hNewton, with v = %.12e, vPhi = %.12e, r = %.12e, Phi = %.12e, l = %d, m = %d\n",
646 vPhi2,
r, cartvalues->data[12] + cartvalues->data[13], (
UINT4)
l,
m,
params->eobParams);
651 if (((
l +
m) % 2) == 0) {
658 XLAL_PRINT_INFO(
"In XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform: Hreal = %e, Slm = %e, eta = %e\n",
Hreal, Slm, eta);
668 gsl_sf_result lnr1, arg1;
670 if (
status != GSL_SUCCESS) {
671 XLALPrintError(
"XLAL Error - %s: Error in GSL function\n", __func__);
675 if (
status != GSL_SUCCESS) {
676 XLALPrintError(
"XLAL Error - %s: Error in GSL function\n", __func__);
679 Tlm = cexp((lnr1.val +
LAL_PI * hathatk) + I * (
680 arg1.val + 2.0 * hathatk * log(4.0 * k / sqrt(
LAL_E))));
686 hathatksq4 = 4. * hathatk * hathatk;
687 hathatk4pi = 4. *
LAL_PI * hathatk;
689 Tlmprefac = sqrt(hathatk4pi / (1. - exp(-hathatk4pi))) /
z2.val;
692 for (Tlmprodfac = 1.,
i = 1;
i <=
l;
i++)
693 Tlmprodfac *= (hathatksq4 + (
REAL8)
i *
i);
696 Tlmold = Tlmprefac * sqrt(Tlmprodfac);
697 XLAL_PRINT_INFO(
"Tlm = %e + i%e, |Tlm| = %.16e (should be %.16e)\n", creal(Tlm), cimag(Tlm), cabs(Tlm), Tlmold);
729 XLAL_PRINT_INFO(
"PK:: rho22v2 = %.12e, rho22v3 = %.12e, rho22v4 = %.12e,\n rho22v5 = %.16e, rho22v6 = %.16e, rho22v6LOG = %.16e, \n rho22v7 = %.12e, rho22v8 = %.16e, rho22v8LOG = %.16e, \n rho22v10 = %.16e, rho22v10LOG = %.16e\n, rho22v6 = %.12e, rho22v8 = %.12e, rho22v10 = %.12e\n",
743 rholm = 1. + v * (hCoeffs->
rho21v1
749 auxflm = v * hCoeffs->
f21v1 + v2 * v * hCoeffs->
f21v3;
765 auxflm = v * v2 * hCoeffs->
f33v3;
770 rholm = 1. + v * (hCoeffs->
rho32v
782 auxflm = v * v2 * hCoeffs->
f31v3;
795 rholm = 1. + v2 * (hCoeffs->
rho44v2
798 + hCoeffs->
rho44v6l * eulerlogxabs) * v))));
803 rholm = 1. + v * (hCoeffs->
rho43v
807 auxflm = v * hCoeffs->
f43v;
811 rholm = 1. + v2 * (hCoeffs->
rho42v2
818 rholm = 1. + v * (hCoeffs->
rho41v
822 auxflm = v * hCoeffs->
f41v;
833 rholm = 1. + v2 * (hCoeffs->
rho55v2
844 rholm = 1. + v2 * (hCoeffs->
rho53v2
854 rholm = 1. + v2 * (hCoeffs->
rho51v2
904 rholm = 1. + hCoeffs->
rho76v2 * v2;
912 rholm = 1. + hCoeffs->
rho74v2 * v2;
920 rholm = 1. + hCoeffs->
rho72v2 * v2;
935 rholm = 1. + hCoeffs->
rho88v2 * v2;
939 rholm = 1. + hCoeffs->
rho87v2 * v2;
943 rholm = 1. + hCoeffs->
rho86v2 * v2;
947 rholm = 1. + hCoeffs->
rho85v2 * v2;
951 rholm = 1. + hCoeffs->
rho84v2 * v2;
955 rholm = 1. + hCoeffs->
rho83v2 * v2;
959 rholm = 1. + hCoeffs->
rho82v2 * v2;
963 rholm = 1. + hCoeffs->
rho81v2 * v2;
979 XLAL_PRINT_INFO(
"exp(delta_%d_%d) = %.16e + i%.16e (abs=%e)\n",
l,
m, creal(cexp(I * deltalm)),
980 cimag(cexp(I * deltalm)), cabs(cexp(I * deltalm)));
995 if (eta == 0.25 &&
m % 2) {
1001 if (
r > 0.0 && debugPK) {
1002 XLAL_PRINT_INFO(
"YP::dynamics variables in waveform: %i, %i, r = %.12e, v = %.12e\n",
l,
m,
r, v);
1003 XLAL_PRINT_INFO(
"rholm^l = %.16e, Tlm = %.16e + i %.16e, \nSlm = %.16e, hNewton = %.16e + i %.16e, delta = %.16e\n", rholmPwrl, creal(Tlm), cimag(Tlm), Slm, creal(hNewton), cimag(hNewton), 0.0);
1006 *hlm = Tlm * cexp(I * deltalm) * Slm * rholmPwrl;
1008 if (
r > 8.5 && debugPK) {
1009 XLAL_PRINT_INFO(
"YP::FullWave: %.16e,%.16e, %.16e\n", creal(*hlm), cimag(*hlm), cabs(*hlm));
1036 const REAL8 values[],
1042 for(
int i =0;
i < 12;
i++)
1043 if( isnan(values[
i]) ) {
1044 XLAL_PRINT_INFO(
"XLALSimIMRSpinPrecEOBCalcOmega_exact::values %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f\n", values[0], values[1], values[2], values[3], values[4], values[5], values[6], values[7], values[8], values[9], values[10], values[11]);
1045 XLALPrintError(
"XLAL Error - %s: nan in input values \n", __func__);
1053 static const REAL8 STEP_SIZE = 1.0e-4;
1059 REAL8 cartValues[14] = {0.}, dvalues[14] = {0.};
1060 REAL8 cartvalues[14] = {0.}, polarvalues[6] = {0.};
1061 REAL8 polarRPcartSvalues[14] = {0.};
1062 memcpy( cartValues, values, 14 *
sizeof(
REAL8) );
1066 REAL8 rvec[3] = {0.,0,0}, pvec[3] = {0.,0,0};
1067 REAL8 s1vec[3] = {0.,0,0}, s2vec[3] = {0.,0,0};
1069 REAL8 rdotvec[3] = {0.,0,0};
1070 REAL8 rvecprime[3] = {0.,0,0}, pvecprime[3] = {0.,0,0},
1071 s1vecprime[3]= {0.,0,0}, s2vecprime[3]= {0.,0,0};
1072 REAL8 rvectmp[3] = {0.,0,0}, pvectmp[3] = {0.,0,0},
1073 s1vectmp[3] = {0.,0,0}, s2vectmp[3]= {0.,0,0};
1074 REAL8 LNhatprime[3]= {0.,0,0}, LNhatTmp[3]= {0.,0,0};
1075 REAL8 rcrossrdot[3] = {0.,0,0};
1077 REAL8 Rot1[3][3] ={{0.}};
1078 REAL8 Rot2[3][3] ={{0.}} ;
1079 REAL8 LNhat[3] = {0.,0,0};
1081 REAL8 Xhat[3] = {1, 0, 0};
1082 UNUSED
REAL8 Yhat[3] = {0, 1, 0};
1083 UNUSED
REAL8 Zhat[3] = {0, 0, 1};
1085 REAL8 Xprime[3] = {0.,0,0}, Yprime[3] = {0.,0,0}, Zprime[3] = {0.,0,0};
1100 memcpy( rvec, values, 3*
sizeof(
REAL8) );
1101 memcpy( pvec, values+3, 3*
sizeof(
REAL8) );
1102 memcpy( s1vec, values+6, 3*
sizeof(
REAL8) );
1103 memcpy( s2vec, values+9, 3*
sizeof(
REAL8) );
1106 memset( dvalues, 0, 14 *
sizeof(
REAL8) );
1112 memcpy( rdotvec, dvalues, 3*
sizeof(
REAL8) );
1115 for(
int ii =0; ii < 12; ii++)
1116 if( isnan(dvalues[ii]) ) {
1117 XLAL_PRINT_INFO(
"XLALSimIMRSpinPrecEOBCalcOmega::dvalues %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f\n", dvalues[0], dvalues[1], dvalues[2], dvalues[3], dvalues[4], dvalues[5], dvalues[6], dvalues[7], dvalues[8], dvalues[9], dvalues[10], dvalues[11]);
1126 for(
i = 0;
i < 3;
i++ )
1127 rcrossrdot[
i] /= rcrossrdotNorm;
1128 memcpy( LNhat, rcrossrdot, 3 *
sizeof(
REAL8) );
1144 Rot1[0][0] = 1.; Rot1[0][1] = 0; Rot1[0][2] = 0;
1145 Rot1[1][0] = 0.; Rot1[1][1] = 1; Rot1[1][2] = 0;
1146 Rot1[2][0] = 0.; Rot1[2][1] = 0; Rot1[2][2] = 1;
1148 memcpy(Xprime, LNhat, 3 *
sizeof(
REAL8));
1152 Rot1[0][0] = 1./sqrt(2); Rot1[0][1] = -1/sqrt(2); Rot1[0][2] = 0;
1153 Rot1[1][0] = 1./sqrt(2); Rot1[1][1] = 1./sqrt(2); Rot1[1][2] = 0;
1154 Rot1[2][0] = 0.; Rot1[2][1] = 0; Rot1[2][2] = 1;
1155 LNhatTmp[0] = LNhatTmp[1] = LNhatTmp[2] = 0.;
1159 LNhatTmp[
i] += Rot1[
i][j]*LNhat[j];
1161 memcpy(Xprime, LNhatTmp, 3*
sizeof(
REAL8));
1167 for(
i=0;
i<3;
i++) Yprime[
i] /= tmpvar;
1171 for(
i=0;
i<3;
i++) Zprime[
i] /= tmpvar;
1173 Rot2[0][0] = Xprime[0]; Rot2[0][1] = Xprime[1]; Rot2[0][2] = Xprime[2];
1174 Rot2[1][0] = Yprime[0]; Rot2[1][1] = Yprime[1]; Rot2[1][2] = Yprime[2];
1175 Rot2[2][0] = Zprime[0]; Rot2[2][1] = Zprime[1]; Rot2[2][2] = Zprime[2];
1177 memset( rvectmp, 0, 3 *
sizeof(
REAL8) );
1178 memset( pvectmp, 0, 3 *
sizeof(
REAL8) );
1179 memset( s1vectmp, 0, 3 *
sizeof(
REAL8) );
1180 memset( s2vectmp, 0, 3 *
sizeof(
REAL8) );
1181 memset( rvecprime, 0, 3 *
sizeof(
REAL8) );
1182 memset( pvecprime, 0, 3 *
sizeof(
REAL8) );
1183 memset( s1vecprime, 0, 3 *
sizeof(
REAL8) );
1184 memset( s2vecprime, 0, 3 *
sizeof(
REAL8) );
1185 memset( LNhatprime, 0, 3 *
sizeof(
REAL8) );
1186 memset( LNhatTmp, 0, 3 *
sizeof(
REAL8) );
1192 rvectmp[
i] += Rot1[
i][j]*rvec[j];
1193 pvectmp[
i] += Rot1[
i][j]*pvec[j];
1194 s1vectmp[
i] += Rot1[
i][j]*s1vec[j];
1195 s2vectmp[
i] += Rot1[
i][j]*s2vec[j];
1196 LNhatTmp[
i] += Rot1[
i][j]*LNhat[j];
1201 rvecprime[
i] += Rot2[
i][j]*rvectmp[j];
1202 pvecprime[
i] += Rot2[
i][j]*pvectmp[j];
1203 s1vecprime[
i] += Rot2[
i][j]*s1vectmp[j];
1204 s2vecprime[
i] += Rot2[
i][j]*s2vectmp[j];
1205 LNhatprime[
i] += Rot2[
i][j]*LNhatTmp[j];
1208 memcpy(cartvalues, rvecprime, 3*
sizeof(
REAL8));
1209 memcpy(cartvalues+3, pvecprime, 3*
sizeof(
REAL8));
1210 memcpy(cartvalues+6, s1vecprime, 3*
sizeof(
REAL8));
1211 memcpy(cartvalues+9, s2vecprime, 3*
sizeof(
REAL8));
1226 polarvalues[1] = acos(rvecprime[0] / polarvalues[0]);
1227 polarvalues[2] = atan2(-rvecprime[1], rvecprime[2]);
1232 REAL8 rvecprime_x_xhat[3] = {0.}, rvecprime_x_xhat_x_rvecprime[3] = {0.};
1234 cross_product(rvecprime_x_xhat, rvecprime, rvecprime_x_xhat_x_rvecprime);
1236 polarvalues[4] = -
inner_product(rvecprime_x_xhat_x_rvecprime, pvecprime)
1237 / polarvalues[0] / sin(polarvalues[1]);
1238 polarvalues[5] = -
inner_product(rvecprime_x_xhat, pvecprime);
1246 memcpy( polarRPcartSvalues, cartvalues, 12*
sizeof(
REAL8));
1247 memcpy( polarRPcartSvalues, polarvalues, 6*
sizeof(
REAL8));
1250 params.values = polarRPcartSvalues;
1259 XLAL_CALLGSL( gslStatus = gsl_deriv_central( &F, polarvalues[5],
1260 STEP_SIZE, &omega, &absErr ) );
1262 if ( gslStatus != GSL_SUCCESS )
1264 XLALPrintError(
"XLAL Error - %s: Failure in GSL function\n", __func__ );
1282 const REAL8 values[],
1288 for(
int i =0;
i < 12;
i++)
1289 if( isnan(values[
i]) ) {
1290 XLAL_PRINT_INFO(
"XLALSimIMRSpinPrecEOBNonKeplerCoeff_exact::values %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f\n",
1291 values[0], values[1], values[2], values[3], values[4], values[5],
1292 values[6], values[7], values[8], values[9], values[10], values[11]);
1298 REAL8 omegaCirc = 0;
1299 REAL8 tmpValues[14]= {0.};
1303 memcpy( tmpValues, values,
sizeof(tmpValues) );
1314 return 1.0/(omegaCirc*omegaCirc*r3);
static UNUSED int XLALSimIMRSpinEOBCalculateNewtonianMultipole(COMPLEX16 *multipole, REAL8 x, UNUSED REAL8 r, REAL8 phi, UINT4 l, INT4 m, EOBParams *params)
This function calculates the Newtonian multipole part of the factorized waveform for the SEOBNRv1 mod...
static UNUSED int XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole(COMPLEX16 *multipole, REAL8 x, UNUSED REAL8 r, UNUSED REAL8 phi, UINT4 l, INT4 m, EOBParams *params)
This function calculates the Newtonian multipole part of the factorized waveform for the SEOBNRv1 mod...
static REAL8 XLALSimIMRSpinAlignedEOBNonKeplerCoeff(const REAL8 values[], SpinEOBParams *funcParams)
Function to calculate the non-Keplerian coefficient for the spin-aligned EOB model.
static REAL8 * cross_product(const REAL8 values1[], const REAL8 values2[], REAL8 result[])
static double GSLSpinPrecHamiltonianWrapperFordHdpphi(double x, void *params)
Wrapper for GSL to call the Hamiltonian function.
static REAL8 inner_product(const REAL8 values1[], const REAL8 values2[])
Functions to compute the inner product and cross products between vectors.
static REAL8 UNUSED z2(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
static int XLALSpinPrecHcapRvecDerivative_exact(double UNUSED t, const REAL8 values[], REAL8 dvalues[], void *funcParams)
Function to calculate exact derivatives of the spin EOB Hamiltonian, to get , as decribed in Eqs.
#define XLAL_CALLGSL(statement)
#define XLAL_ERROR_REAL8(...)
#define XLAL_PRINT_INFO(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_IS_REAL8_FAIL_NAN(val)
Parameters for the spinning EOB model.