Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimIMRSpinEOBFactorizedFluxPrec_v3opt.c
Go to the documentation of this file.
1/**
2 * \author Craig Robinson, Yi Pan, Prayush Kumar, Stas Babak, Andrea Taracchini
3 */
4
5#ifndef _LALSIMIMRSPINPRECEOBFACTORIZEDFLUX_V3OPT_C
6#define _LALSIMIMRSPINPRECEOBFACTORIZEDFLUX_V3OPT_C
7
8
9#include <complex.h>
10#include <lal/LALSimInspiral.h>
11#include <lal/LALSimIMR.h>
12
13#include "LALSimIMREOBNRv2.h"
14#include "LALSimIMRSpinEOB.h"
15
19
20
21/*------------------------------------------------------------------------------------------
22 *
23 * Prototypes of functions defined in this code.
24 *
25 *------------------------------------------------------------------------------------------
26 */
27
28static REAL8
30 REAL8Vector * polvalues, /**< \f$(r,\phi,p_r,p_\phi)\f$ */
31 REAL8Vector * values, /**< dynamical variables */
32 EOBNonQCCoeffs * nqcCoeffs, /**< pre-computed NQC coefficients */
33 const REAL8 omega, /**< orbital frequency */
34 SpinEOBParams * ak, /**< physical parameters */
35 const REAL8 H, /**< real Hamiltonian */
36 const INT4 lMax, /**< upper limit of the summation over l */
37 const UINT4 SpinAlignedEOBversion /**< 1 for SEOBNRv1, 2 for SEOBNRv2 */
38);
39/*------------------------------------------------------------------------------------------
40 *
41 * Defintions of functions.
42 *
43 *------------------------------------------------------------------------------------------
44 */
45/**
46 * This function calculates the spin factorized-resummed GW energy flux
47 * for given dynamical variables.
48 * Eq. 12 of PRD 86, 024011 (2012)
49 */
50static REAL8
52 REAL8Vector * polvalues, /**< \f$(r,\phi,p_r,p_\phi)\f$ */
53 REAL8Vector * values, /**< dynamical variables */
54 EOBNonQCCoeffs * nqcCoeffs, /**< pre-computed NQC coefficients */
55 const REAL8 omega, /**< orbital frequency */
56 SpinEOBParams * ak, /**< physical parameters */
57 const REAL8 H, /**< real Hamiltonian */
58 const INT4 lMax, /**< upper limit of the summation over l */
59 const UINT4 SpinAlignedEOBversion /**< 1 for SEOBNRv1, 2 for SEOBNRv2 */
60)
61{
62 int debugPK = 0;
63 int i = 0;
64 double radius = sqrt(values->data[0]*values->data[0] + values->data[1] *values->data[1] + values->data[2] *values->data[2] );
65 if (radius < 1.) {
66 return 0.;
67 }
68 if (1){
69 for( i =0; i < 4; i++)
70 if( isnan(polvalues->data[i]) ) {
71 XLAL_PRINT_INFO("XLALInspiralPrecSpinFactorizedFlux (from input)::polvalues %3.10f %3.10f %3.10f %3.10f\n", polvalues->data[0], polvalues->data[1], polvalues->data[2], polvalues->data[3]);
72 XLALPrintError( "XLAL Error - %s: nan polvalues: %3.10f %3.10f %3.10f %3.10f \n", __func__, polvalues->data[0], polvalues->data[1], polvalues->data[2], polvalues->data[3] );
74 }
75
76 for( i =0; i < 12; i++)
77 if( isnan(values->data[i]) ) {
78 XLAL_PRINT_INFO("XLALInspiralPrecSpinFactorizedFlux (from input)::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->data[0], values->data[1], values->data[2], values->data[3], values->data[4], values->data[5], values->data[6], values->data[7], values->data[8], values->data[9], values->data[10], values->data[11]);
79 XLALPrintError( "XLAL Error - %s: nan in input 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", __func__, values->data[0], values->data[1], values->data[2], values->data[3], values->data[4], values->data[5], values->data[6], values->data[7], values->data[8], values->data[9], values->data[10], values->data[11] );
81 }
82 }
83
84 REAL8 flux = 0.0;
85 REAL8 v;
86 REAL8 omegaSq;
87 COMPLEX16 hLM;
88 INT4 l , m;
89
90 //EOBNonQCCoeffs nqcCoeffs;
91
92 if (!values || !ak) {
94 }
95
96 if (lMax < 2) {
98 }
99 /* Omega is the derivative of phi */
100 omegaSq = omega * omega;
101
102 v = cbrt(omega);
103
104 /* Update the factorized multipole coefficients, w.r.t. new spins */
105 if (0) { /* {{{ */
106 XLAL_PRINT_INFO("\nValues inside Flux:\n");
107 for (i = 0; i < 11; i++)
108 XLAL_PRINT_INFO("values[%d] = %.12e\n", i, values->data[i]);
109 /*
110 * Assume that initial conditions are available at this
111 * point, to compute the chiS and chiA parameters. Calculate
112 * the values of chiS and chiA, as given in Eq.16 of
113 * Precessing EOB paper Pan et.al. arXiv:1307.6232 (or PRD 89, 084006 (2014)). Assuming \vec{L} to be pointing in
114 * the direction of \vec{r}\times\vec{p}
115 */
116 REAL8 rcrossp [3], rcrosspMag, s1dotL, s2dotL;
117 REAL8 chiS , chiA, tplspin;
118
119 rcrossp[0] = values->data[1] * values->data[5] - values->data[2] * values->data[4];
120 rcrossp[1] = values->data[2] * values->data[3] - values->data[0] * values->data[5];
121 rcrossp[2] = values->data[0] * values->data[4] - values->data[1] * values->data[3];
122 rcrosspMag = sqrt(rcrossp[0] * rcrossp[0] + rcrossp[1] * rcrossp[1] +
123 rcrossp[2] * rcrossp[2]);
124
125 rcrossp[0] /= rcrosspMag;
126 rcrossp[1] /= rcrosspMag;
127 rcrossp[2] /= rcrosspMag;
128
129 s1dotL = values->data[6] * rcrossp[0] + values->data[7] * rcrossp[1]
130 + values->data[8] * rcrossp[2];
131 s2dotL = values->data[9] * rcrossp[0] + values->data[10] * rcrossp[1]
132 + values->data[11] * rcrossp[2];
133
134 chiS = 0.5 * (s1dotL + s2dotL);
135 chiA = 0.5 * (s1dotL - s2dotL);
136
137 /*
138 * Compute the test-particle limit spin of the deformed-Kerr
139 * background
140 */
141 switch (SpinAlignedEOBversion) {
142 case 1:
143 tplspin = 0.0;
144 break;
145 case 2:
146 tplspin = (1. - 2. * ak->eobParams->eta) * chiS + (ak->eobParams->m1
147 - ak->eobParams->m2) / (ak->eobParams->m1 + ak->eobParams->m2) * chiA;
148 break;
149 default:
150 XLALPrintError("XLAL Error - %s: Unknown SEOBNR version!\nAt present only v1 and v2 are available.\n", __func__);
152 break;
153 }
154
155 /* ************************************************* */
156 /* Re-Populate the Waveform structures */
157 /* ************************************************* */
158
159 /* Re-compute the spinning coefficients for hLM */
160 //debugPK
161 XLAL_PRINT_INFO("Re-calculating waveform coefficients in the Flux function with chiS, chiA = %e, %e!\n", chiS, chiA);
162 chiS = 0.3039435650957116;
163 chiA = -0.2959424290852973;
164 XLAL_PRINT_INFO("Changed them to the correct values = %e, %e!\n", chiS, chiA);
165
166 if (ak->alignedSpins==1) {
168 ak->eobParams->m1, ak->eobParams->m2, ak->eobParams->eta,
169 tplspin, chiS, chiA, SpinAlignedEOBversion) == XLAL_FAILURE) {
172 }
173 }
174 else {
176 ak->eobParams->m1, ak->eobParams->m2, ak->eobParams->eta,
177 tplspin, chiS, chiA, 3) == XLAL_FAILURE) {
180 }
181 }
182 } /* }}} */
183 //XLAL_PRINT_INFO("v = %.16e\n", v);
184 COMPLEX16 hLMTab[lMax+1][lMax+1];
185 if (XLALSimIMRSpinEOBFluxGetPrecSpinFactorizedWaveform_exact(&hLMTab[0][0], polvalues, values, v, H,
186 lMax, ak) == XLAL_FAILURE) {
188 }
189 for (l = 2; l <= lMax; l++) {
190
191 for (m = 1; m <= l; m++) {
192
193 if (debugPK)
194 XLAL_PRINT_INFO("\nGetting (%d, %d) mode for flux!\n", l, m);
195 //XLAL_PRINT_INFO("Stas, computing the waveform l = %d, m =%d\n", l, m);
196 hLM = hLMTab[l][m];
197 //XLAL_PRINT_INFO("Stas: done\n");
198 /*
199 * For the 2,2 mode, we apply NQC correction to the
200 * flux
201 */
202 if (l == 2 && m == 2) {
203 COMPLEX16 hNQC;
204 /*
205 * switch ( SpinAlignedEOBversion ) { case 1:
206 * XLALSimIMRGetEOBCalibratedSpinNQC(
207 * &nqcCoeffs, l, m, ak->eobParams->eta,
208 * ak->a ); break; case 2: //
209 * XLALSimIMRGetEOBCalibratedSpinNQCv2(
210 * &nqcCoeffs, l, m, ak->eobParams->eta,
211 * ak->a );
212 * XLALSimIMRGetEOBCalibratedSpinNQC3D(
213 * &nqcCoeffs, l, m, ak->eobParams->eta,
214 * ak->a, (ak->chi1 - ak->chi2)/2. ); break;
215 * default: XLALPrintError( "XLAL Error - %s:
216 * Unknown SEOBNR version!\nAt present only
217 * v1 and v2 are available.\n", __func__);
218 * XLAL_ERROR( XLAL_EINVAL ); break; }
219 */
220 if (debugPK)
221 XLAL_PRINT_INFO("\tl = %d, m = %d, NQC: a1 = %.16e, a2 = %.16e, a3 = %.16e, a3S = %.16e, a4 = %.16e, a5 = %.16e\n\tb1 = %.16e, b2 = %.16e, b3 = %.16e, b4 = %.16e\n",
222 2, 2, nqcCoeffs->a1, nqcCoeffs->a2, nqcCoeffs->a3, nqcCoeffs->a3S, nqcCoeffs->a4, nqcCoeffs->a5,
223 nqcCoeffs->b1, nqcCoeffs->b2, nqcCoeffs->b3, nqcCoeffs->b4);
224 XLALSimIMREOBNonQCCorrection(&hNQC, polvalues, omega, nqcCoeffs);
225 if (debugPK)
226 XLAL_PRINT_INFO("\tl = %d, m = %d, hNQC = %.16e + i%.16e, |hNQC| = %.16e\n", l, m,
227 creal(hNQC), cimag(hNQC), sqrt(creal(hNQC) * creal(hNQC) + cimag(hLM) * cimag(hLM)));
228
229 if((m * m) * omegaSq * (creal(hLM) * creal(hLM) + cimag(hLM) * cimag(hLM)) > 5.) {
230
231 XLAL_PRINT_INFO("\tl = %d, m = %d, mag(hLM) = %.17e, mag(hNQC) = %.17e, omega = %.16e\n",
232 l, m, sqrt(creal(hLM) * creal(hLM) + cimag(hLM) * cimag(hLM)),
233 sqrt(creal(hNQC) * creal(hNQC) + cimag(hNQC) * cimag(hNQC)), omega);
234
235 XLAL_PRINT_INFO("XLALInspiralPrecSpinFactorizedFlux (from input)::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->data[0], values->data[1], values->data[2], values->data[3], values->data[4], values->data[5], values->data[6], values->data[7], values->data[8], values->data[9], values->data[10], values->data[11]);
236 }
237
238 /* Eq. 16 */
239 //FIXME
240 hLM *= hNQC;
241 }
242 if (debugPK)
243 XLAL_PRINT_INFO("\tl = %d, m = %d, mag(hLM) = %.17e, omega = %.16e\n", l, m, sqrt(creal(hLM) * creal(hLM) + cimag(hLM) * cimag(hLM)), omega);
244
245 /* Eq. 13 */
246 flux += (REAL8) (m * m) * omegaSq * (creal(hLM) * creal(hLM) + cimag(hLM) * cimag(hLM));
247 }
248 }
249 if( (omegaSq > 1 || flux > 5) ) {
250 if(debugPK) {
251 XLAL_PRINT_INFO("In XLALInspiralPrecSpinFactorizedFlux: omegaSq = %3.12f, FLUX = %3.12f, r = %3.12f\n",
252 omegaSq, flux,radius);
253 }
254 flux = 0.;
255 }
256
257 if (debugPK)
258 XLAL_PRINT_INFO("\tStas, FLUX = %.16e\n", flux * LAL_1_PI / 8.0);
259 return flux * LAL_1_PI / 8.0;
260}
261#endif /* _LALSIMIMRSPINPRECEOBFACTORIZEDFLUX_C */
static UNUSED int XLALSimIMREOBNonQCCorrection(COMPLEX16 *restrict nqc, REAL8Vector *restrict values, const REAL8 omega, EOBNonQCCoeffs *restrict coeffs)
This function calculates the non-quasicircular correction to apply to the waveform.
static REAL8 XLALInspiralPrecSpinFactorizedFlux_exact(REAL8Vector *polvalues, REAL8Vector *values, EOBNonQCCoeffs *nqcCoeffs, const REAL8 omega, SpinEOBParams *ak, const REAL8 H, const INT4 lMax, const UINT4 SpinAlignedEOBversion)
This function calculates the spin factorized-resummed GW energy flux for given dynamical variables.
static int XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients(FacWaveformCoeffs *const coeffs, const REAL8 m1, const REAL8 m2, const REAL8 eta, const REAL8 a, const REAL8 chiS, const REAL8 chiA, const UINT4 SpinAlignedEOBversion)
Waveform expressions are given by Taracchini et al.
static INT4 XLALSimIMRSpinEOBFluxGetPrecSpinFactorizedWaveform_exact(COMPLEX16 *restrict hlmTab, REAL8Vector *restrict values, REAL8Vector *restrict cartvalues, const REAL8 v, const REAL8 Hreal, const INT4 lMax, SpinEOBParams *restrict params)
FOR PRECESSING EOB This function calculates hlm mode factorized-resummed waveform for given dynamical...
int l
Definition: bh_qnmode.c:135
double i
Definition: bh_ringdown.c:118
const double H
#define LAL_1_PI
double complex COMPLEX16
double REAL8
uint32_t UINT4
int32_t INT4
static const INT4 m
void XLALDestroyREAL8Vector(REAL8Vector *vector)
#define XLAL_ERROR_REAL8(...)
#define XLAL_PRINT_INFO(...)
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_EFAULT
XLAL_EFUNC
XLAL_EINVAL
XLAL_FAILURE
The coefficients which are used in calculating the non-quasicircular correction to the EOBNRv2 model.
FacWaveformCoeffs * hCoeffs
REAL8 * data
Parameters for the spinning EOB model.
EOBParams * eobParams