Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-5e288d3
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimIMRSpinEOBHcapExactDerivativePrec_v3opt.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2014 Craig Robinson, Enrico Barausse, Yi Pan, Prayush Kumar,
3 * Stanislav Babak, Andrea Taracchini
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with with program; see the file COPYING. If not, write to the
17 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18 * MA 02110-1301 USA
19 */
20/**
21 * \author Craig Robinson, Yi Pan, Stas Babak, Prayush Kumar, Andrea Taracchini
22 * \brief Functions to compute the RHSs of the SEOBNRv3 ODEs
23 */
24
25#ifndef _LALSIMIMRSPINPRECEOBHCAPEXACTDERIVATIVEPREC_V3OPT_C
26#define _LALSIMIMRSPINPRECEOBHCAPEXACTDERIVATIVEPREC_V3OPT_C
27
28#ifdef __GNUC__
29#define UNUSED __attribute__ ((unused))
30#else
31#define UNUSED
32#endif
33
34#include <math.h>
35#include <gsl/gsl_deriv.h>
36#include <lal/LALSimInspiral.h>
37#include <lal/LALSimIMR.h>
38
39#include "LALSimIMRSpinEOB.h"
40
47
49
50/*------------------------------------------------------------------------------------------
51 *
52 * Prototypes of functions defined in this code.
53 *
54 *------------------------------------------------------------------------------------------
55 */
56
58 double UNUSED t, /**<< UNUSED */
59 const REAL8 values[], /**<< Dynamical variables */
60 REAL8 dvalues[], /**<< Time derivatives of variables (returned) */
61 void *funcParams /**<< EOB parameters */
62);
63
64
65/*------------------------------------------------------------------------------------------
66 *
67 * Defintions of functions.
68 *
69 *------------------------------------------------------------------------------------------
70 */
71
72/**
73 * Function to calculate numerical derivatives of the spin EOB Hamiltonian,
74 * which correspond to time derivatives of the dynamical variables in conservative dynamcis.
75 * All derivatives are returned together.
76 * The derivatives are combined with energy flux to give right hand side of the ODEs
77 * of a generic spin EOB model, as decribed in Eqs. A4, A5, 26 and 27 of
78 * Pan et al. PRD 81, 084041 (2010)
79 * Later on we use BB1, that is Barausse and Buonanno PRD 81, 084024 (2010)
80 * This function is not used by the spin-aligned model.
81 */
83 double UNUSED t, /**<< UNUSED */
84 const REAL8 values[], /**<< Dynamical variables */
85 REAL8 dvalues[], /**<< Time derivatives of variables (returned) */
86 void *funcParams /**<< EOB parameters */
87)
88{
89 int debugPK = 0;
90
91 /** lMax: l index up to which h_{lm} modes are included in the computation of the GW enegy flux: see Eq. in 13 in PRD 86, 024011 (2012) */
92 static const INT4 lMax = 8;
94
95 /* Since we take numerical derivatives wrt dynamical variables */
96 /* but we want them wrt time, we use this temporary vector in */
97 /* the conversion */
98 REAL8 tmpDValues[14] = {0};
99
100 //Hamiltonian
101 REAL8 H = 0.0;
102 REAL8 flux;
103
104 /* Set up pointers for GSL */
105 params.values = values;
106 params.params = (SpinEOBParams *) funcParams;
107
108 UINT4 SpinAlignedEOBversion;
109 UINT4 i, j, k, l;
110
111 REAL8Vector rVec, pVec;
112 REAL8 rData[3], pData[3];
113
114 /* We need r, phi, pr, pPhi to calculate the flux */
115 REAL8 r;
116 REAL8Vector polarDynamics, cartDynamics;
117 REAL8 polData[4];
118
119 REAL8 mass1, mass2, eta;
120 UNUSED REAL8 rrTerm2, pDotS1, pDotS2;
121 REAL8Vector s1, s2, s1norm, s2norm, sKerr, sStar;
122 REAL8 s1Data[3], s2Data[3], s1DataNorm[3], s2DataNorm[3];
123 REAL8 sKerrData[3], sStarData[3];
124 REAL8 chiS, chiA, a, tplspin;
125 REAL8 s1dotLN, s2dotLN;
126
127 /* Orbital angular momentum */
128 REAL8 Lx, Ly, Lz, magL;
129 REAL8 Lhatx, Lhaty, Lhatz;
130 REAL8 dLx, dLy, dLz;
131 REAL8 dLhatx, dLhaty, dMagL;
132
133 REAL8 alphadotcosi;
134
135 REAL8 rCrossV_x, rCrossV_y, rCrossV_z, omega;
136
137 REAL8 tmpP[3], rMag, rMag2, prT;
138 REAL8 u, u2, u3, u4, u5, w2, a2;
140 REAL8 eobD_r, deltaU_u, deltaU_r;
141 REAL8 dcsi, csi;
142
143 REAL8 tmpValues[12];
144 REAL8 Tmatrix[3][3], invTmatrix[3][3], dTijdXk[3][3][3];
145 REAL8 tmpPdotT1[3], tmpPdotT2[3], tmpPdotT3[3];
146 //3 terms of Eq.A5
147
148 /* NQC coefficients container */
149 EOBNonQCCoeffs * nqcCoeffs = NULL;
150
151 nqcCoeffs = params.params->nqcCoeffs;
152
153 mass1 = params.params->eobParams->m1;
154 mass2 = params.params->eobParams->m2;
155 eta = params.params->eobParams->eta;
156 SpinAlignedEOBversion = params.params->seobCoeffs->SpinAlignedEOBversion;
157 SpinEOBHCoeffs *coeffs = (SpinEOBHCoeffs *) params.params->seobCoeffs;
158
159 /*
160 * For precessing binaries, the effective spin of the Kerr background
161 * evolves with time. The coefficients used to compute the
162 * Hamiltonian depend on the Kerr spin, and hence need to be updated
163 * for the current spin values
164 */
165
166 /*
167 * Set the position/momenta vectors to point to the appropriate
168 * things.
169 * Here pvec is the reduced tortoise p^* vector of Pan et al.
170 * PRD 81, 084041 (2010)
171 */
172 rVec.length = pVec.length = 3;
173 rVec.data = rData;
174 pVec.data = pData;
175 memcpy(rData, values, sizeof(rData));
176 memcpy(pData, values + 3, sizeof(pData));
177
178 /*
179 * We need to re-calculate the parameters at each step as precessing
180 * spins will not be constant
181 */
182
183 /*
184 * We cannot point to the values vector directly as it leads to a
185 * warning
186 */
187 s1.length = s2.length = s1norm.length = s2norm.length = 3;
188 s1.data = s1Data;
189 s2.data = s2Data;
190 s1norm.data = s1DataNorm;
191 s2norm.data = s2DataNorm;
192
193 memcpy(s1Data, values + 6, 3 * sizeof(REAL8));
194 memcpy(s2Data, values + 9, 3 * sizeof(REAL8));
195 memcpy(s1DataNorm, values + 6, 3 * sizeof(REAL8));
196 memcpy(s2DataNorm, values + 9, 3 * sizeof(REAL8));
197
198 const REAL8 mass1mass2inv = 1.0/(mass1*mass2);
199 const REAL8 mass1inv = mass1mass2inv*mass2;
200 const REAL8 mass2inv = mass1mass2inv*mass1;
201 const REAL8 totalMass = mass1 + mass2;
202 const REAL8 totalMassinv = 1.0/totalMass;
203 const REAL8 totalMass2 = totalMass*totalMass;
204 const REAL8 totalMass2inv = totalMassinv*totalMassinv;
205
206 for (i = 0; i < 3; i++) {
207 s1.data[i] *= totalMass2;
208 s2.data[i] *= totalMass2;
209 }
210
211 sKerr.length = 3;
212 sKerr.data = sKerrData;
213
214 sStar.length = 3;
215 sStar.data = sStarData;
216
217 for ( i = 0; i < 3; i++ )
218 {
219 sKerr.data[i] = (s1.data[i] + s2.data[i])*totalMass2inv;
220 sStar.data[i] = (mass2*mass1inv * s1.data[i] + mass1*mass2inv * s2.data[i])*totalMass2inv;
221 }
222
223 a = sqrt(sKerr.data[0] * sKerr.data[0] + sKerr.data[1] * sKerr.data[1]
224 + sKerr.data[2] * sKerr.data[2]);
225
226 /* Convert momenta to p Eq. A3 of PRD 81, 084041 (2010) */
227 rMag = sqrt(rVec.data[0] * rVec.data[0] + rVec.data[1] * rVec.data[1] + rVec.data[2] * rVec.data[2]);
228
229 rMag2 = rMag * rMag;
230 u = 1. / rMag;
231 u2 = u * u;
232 u3 = u2 * u;
233 u4 = u2 * u2;
234 u5 = u4 * u;
235 a2 = a * a;
236 w2 = rMag2 + a2;
237
238 /* This is p^*.r/|r| */
239 prT = pVec.data[0] * (rVec.data[0] * u) + pVec.data[1] * (rVec.data[1] * u) + pVec.data[2] * (rVec.data[2] * u);
240
241 /* Eq. 5.83 of BB1, inverse */
242 D = 1. + log(1. + 6. * eta * u2 + 2. * (26. - 3. * eta) * eta * u3);
243 /* d(Eq. 5.83 of BB1)/dr */
244 eobD_r = (u2 / (D * D)) * (12. * eta * u + 6. * (26. - 3. * eta) * eta * u2) / (1.
245 + 6. * eta * u2 + 2. * (26. - 3. * eta) * eta * u3);
246 m1PlusetaKK = -1. + eta * coeffs->KK;
247 const double m1PlusetaKKinv = 1.0/m1PlusetaKK;
248 const double logu = log(u);
249
250 /* Eq. 5.75 of BB1 */
251 /* a as returned by XLALSimIMRSpinEOBCalculateSigmaKerr is S/M^2 so that a is unitless, i.e. 1/M^2 is absorbed in a2. Also, u = M/|r| is unitless */
252 bulk = m1PlusetaKKinv*m1PlusetaKKinv + (2. * u) * m1PlusetaKKinv + a2 * u2;
253 /* Eq. 5.73 of BB1 */
254 /* The 4PN term with coefficients k5 and k5l are defined in the SEOBNRv2 review document https://dcc.ligo.org/T1400476 */
255 logarg = 1. + coeffs->k1 * u + coeffs->k2 * u2 + coeffs->k3 * u3 + coeffs->k4 * u4 + coeffs->k5 * u5 + coeffs->k5l * u5 * log(u);
256 logTerms = 1. + eta * coeffs->k0 + eta * log(logarg);
257
258 /* Eq. 5.73 of BB1 */
259 deltaU = bulk * logTerms;
260 /* Eq. 5.71 of BB1 */
261 deltaT = rMag2 * deltaU;
262 /* ddeltaU/du */
263 /* The log(u) is treated as a constant when taking the derivative wrt u */
264 deltaU_u = 2. * (1.*m1PlusetaKKinv + a2 * u) * logTerms + bulk * (eta * (coeffs->k1 + u * (2. * coeffs->k2 + u * (3. * coeffs->k3
265 + u * (4. * coeffs->k4 + 5. * (coeffs->k5 + coeffs->k5l * logu) * u)))))/(1. + coeffs->k1 * u + coeffs->k2 * u2 + coeffs->k3 * u3
266 + coeffs->k4 * u4 + (coeffs->k5 + coeffs->k5l * logu) * u5);
267 deltaU_r = -u2 * deltaU_u;
268 /* Eq. 5.38 of BB1 */
269 deltaR = deltaT * D;
270 if (params.params->tortoise)
271 csi = sqrt(deltaT * deltaR) / w2; /* Eq. 28 of Pan et al.
272 * PRD 81, 084041 (2010) */
273 else
274 csi = 1.0;
275
276 /* This is A3 of PRD 81, 084041 (2010) explicitly evaluated */
277 for (i = 0; i < 3; i++) {
278 tmpP[i] = pVec.data[i] - (rVec.data[i] * u) * prT * (csi - 1.) / csi;
279 }
280
281 if (debugPK) {
282 XLAL_PRINT_INFO("csi = %.12e\n", csi);
283 for (i = 0; i < 3; i++)
284 XLAL_PRINT_INFO("p,p*: %.12e\t%.12e\n", pData[i], tmpP[i]);
285 }
286 /*
287 * Calculate the T-matrix, required to convert P from tortoise to
288 * non-tortoise coordinates, and/or vice-versa. This is given
289 * explicitly in Eq. A3 of 0912.3466
290 */
291 for (i = 0; i < 3; i++)
292 for (j = 0; j <= i; j++) {
293 Tmatrix[i][j] = Tmatrix[j][i] = (rVec.data[i] * rVec.data[j] * u2) * (csi - 1.);
294 invTmatrix[i][j] = invTmatrix[j][i] = -(csi - 1) / csi * (rVec.data[i] * rVec.data[j] * u2);
295 if (i == j) {
296 Tmatrix[i][j]++;
297 invTmatrix[i][j]++;
298 }
299 }
300 /* This is dcsi/dr: this is needed in the last term of Eq. A5 of PRD 81, 084041 (2010) */
301 dcsi = csi * (2. * u + deltaU_r / deltaU) + csi * csi * csi
302 / (2. * rMag2 * rMag2 * deltaU * deltaU) * (rMag * (-4. * w2) / D - eobD_r * (w2 * w2));
303
304 for (i = 0; i < 3; i++)
305 for (j = 0; j < 3; j++)
306 for (k = 0; k < 3; k++) {
307 dTijdXk[i][j][k] = (rVec.data[i] * KRONECKER_DELTA(j, k) + KRONECKER_DELTA(i, k) * rVec.data[j])
308 * (csi - 1.) * u2
309 + rVec.data[i] * rVec.data[j] * rVec.data[k] * u2 * u * (-2. * u * (csi - 1.) + dcsi);
310 }
311
312 //Print out the T - matrix for comparison
313 if (debugPK) {
314 XLAL_PRINT_INFO("\nT-Matrix:\n");
315 for (i = 0; i < 3; i++)
316 XLAL_PRINT_INFO("%le\t%le\t%le\n", Tmatrix[i][0], Tmatrix[i][1], Tmatrix[i][2]);
317
318 for (i = 0; i < 3; i++) {
319 XLAL_PRINT_INFO("dT[%d][j]dX[k]:\n", i);
320 for (j = 0; j < 3; j++)
321 XLAL_PRINT_INFO("%.12e\t%.12e\t%.12e\n", dTijdXk[i][j][0],
322 dTijdXk[i][j][1], dTijdXk[i][j][2]);
323 XLAL_PRINT_INFO("\n");
324 }
325 }
326 INT4 updateHCoeffsOld = params.params->seobCoeffs->updateHCoeffs;
327
328 memcpy(tmpValues, params.values, sizeof(tmpValues));
329 tmpValues[3] = tmpP[0];
330 tmpValues[4] = tmpP[1];
331 tmpValues[5] = tmpP[2];
332 params.params->seobCoeffs->updateHCoeffs = 1;
333
334 XLALSEOBNRv3_opt_ComputeHamiltonianDerivatives(tmpValues,values,params.params->seobCoeffs,tmpDValues,params.params,&H);
335
336 params.params->seobCoeffs->updateHCoeffs = updateHCoeffsOld;
337 /* Now make the conversion */
338 /* rVectorDot */
339 // Eq. A4 of PRD 81, 084041 (2010). Note that dvalues[i] = \dot{X^i} but tmpDValues[j] = dH/dvalues[j]
340 for (i = 0; i < 3; i++)
341 for (j = 0, dvalues[i] = 0.; j < 3; j++)
342 dvalues[i] += tmpDValues[j + 3] * Tmatrix[i][j];
343
344 /* Calculate the orbital angular momentum */
345 Lx = values[1] * values[5] - values[2] * values[4];
346 Ly = values[2] * values[3] - values[0] * values[5];
347 Lz = values[0] * values[4] - values[1] * values[3];
348
349 magL = sqrt(Lx * Lx + Ly * Ly + Lz * Lz);
350
351 Lhatx = Lx / magL;
352 Lhaty = Ly / magL;
353 Lhatz = Lz / magL;
354
355 /* Calculate the polar data */
356 polarDynamics.length = 4;
357 polarDynamics.data = polData;
358
359 r = polarDynamics.data[0] = sqrt(values[0] * values[0] + values[1] * values[1]
360 + values[2] * values[2]);
361 polarDynamics.data[1] = 0;
362 /* Note that poldata[2] and poldata[3] differ from v2 in which one is normalized by polData[0]. Behaviour is reverted. */
363 polarDynamics.data[2] = (values[0] * values[3] + values[1] * values[4]
364 + values[2] * values[5]) / polData[0];
365 polarDynamics.data[3] = magL;
366
367
368 /* Now calculate rCrossRdot and omega */
369 rCrossV_x = values[1] * dvalues[2] - values[2] * dvalues[1];
370 rCrossV_y = values[2] * dvalues[0] - values[0] * dvalues[2];
371 rCrossV_z = values[0] * dvalues[1] - values[1] * dvalues[0];
372
373 omega = sqrt(rCrossV_x * rCrossV_x + rCrossV_y * rCrossV_y + rCrossV_z * rCrossV_z) / (r * r);
374
375 //
376 // Compute \vec{L_N} = \vec{r} \times \.{\vec{r}}, \vec{S_i} \dot
377 // \vec{L_N} and chiS and chiA
378 //
379
380 /* Eq. 16 of PRD 89, 084006 (2014): it's S_{1,2}/m_{1,2}^2.LNhat */
381 s1dotLN = (s1.data[0] * rCrossV_x + s1.data[1] * rCrossV_y + s1.data[2] * rCrossV_z) /
382 (r * r * omega * mass1 * mass1);
383 s2dotLN = (s2.data[0] * rCrossV_x + s2.data[1] * rCrossV_y + s2.data[2] * rCrossV_z) /
384 (r * r * omega * mass2 * mass2);
385
386 chiS = 0.5 * (s1dotLN + s2dotLN);
387 chiA = 0.5 * (s1dotLN - s2dotLN);
388
389 if (debugPK) {
390 XLAL_PRINT_INFO("chiS = %.12e, chiA = %.12e\n", chiS, chiA);
391 fflush(NULL);
392 }
393 if (debugPK) {
394 XLAL_PRINT_INFO("Computing derivatives for values\n%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\n\n",
395 (double)values[0], (double)values[1], (double)values[2],
396 (double)values[3], (double)values[4], (double)values[5],
397 (double)values[6], (double)values[7], (double)values[8],
398 (double)values[9], (double)values[10], (double)values[11]);
399 XLAL_PRINT_INFO("tmpDvalues\n%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t\n",
400 (double)tmpDValues[0], (double)tmpDValues[1], (double)tmpDValues[2],
401 (double)tmpDValues[3], (double)tmpDValues[4], (double)tmpDValues[5],
402 (double)tmpDValues[6], (double)tmpDValues[7], (double)tmpDValues[8],
403 (double)tmpDValues[9], (double)tmpDValues[10], (double)tmpDValues[11]);
404 }
405 /*
406 * Compute the test-particle limit spin of the deformed-Kerr
407 * background
408 */
409 /* /\* See below Eq. 4 of PRD 89, 061502(R) (2014)*\/ */
410 tplspin = (1. - 2. * eta) * chiS + (mass1 - mass2) * totalMassinv * chiA;
411
412 memcpy(params.params->s1Vec->data, s1norm.data, 3*sizeof(*params.params->s1Vec->data));
413 memcpy(params.params->s2Vec->data, s2norm.data, 3*sizeof(*params.params->s2Vec->data));
414 memcpy(params.params->sigmaStar->data, sStar.data, 3*sizeof(*params.params->sigmaStar->data));
415 memcpy(params.params->sigmaKerr->data, sKerr.data, 3*sizeof(*params.params->sigmaKerr->data));
416
417 params.params->a = a;
418
419 XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients(params.params->eobParams->hCoeffs, mass1, mass2, eta, tplspin,
420 chiS, chiA, 3);
421
422 H = H * (mass1 + mass2);
423
424 /* Now we have the ingredients to compute the flux */
425 memcpy(tmpValues, values, 12 * sizeof(REAL8));
426 cartDynamics.data = tmpValues;
427 if (debugPK) {
428 XLAL_PRINT_INFO("params.params->a = %.12e, %.12e\n", a, params.params->a);
429 fflush(NULL);
430 }
431 /* Eq. 13 of PRD 86, 024011 (2012) */
432 if ( params.params->ignoreflux == 0 ) {
433 flux = XLALInspiralPrecSpinFactorizedFlux_exact(&polarDynamics, &cartDynamics,
434 nqcCoeffs, omega, params.params, H * totalMassinv, lMax, SpinAlignedEOBversion);
435 }
436 else if ( params.params->ignoreflux == 1) {
437 flux = 0.;
438 }
439 else {
440 XLAL_PRINT_INFO("Wrong ignorflux option in XLALSpinPrecHcapExactDerivative!\n");
442 }
443 /*
444 * Looking at consistency with the non-spinning model, we have to divide the
445 * flux by eta
446 */
447 flux = flux / eta;
448 pDotS1 = pData[0] * s1.data[0] + pVec.data[1] * s1.data[1] + pVec.data[2] * s1.data[2];
449 pDotS2 = pVec.data[0] * s2.data[0] + pVec.data[1] * s2.data[1] + pVec.data[2] * s2.data[2];
450 const double omega2=omega*omega;
451 const double omega4=omega2*omega2;
452 const double omega8=omega4*omega4;
453 const double omega8o3=cbrt(omega8);
454 rrTerm2 = 8. / 15. * eta * eta * omega8o3 / (magL * magL * r) * ((61. + 48. * mass2 * mass1inv) * pDotS1 + (61. + 48. * mass1 * mass2inv) * pDotS2);
455
456 if (debugPK) {
457 XLAL_PRINT_INFO("omega = %.12e \n flux = %.12e \n Lmag = %.12e\n", omega, flux, magL);
458 XLAL_PRINT_INFO("rrForce = %.12e %.12e %.12e\n", -flux * values[3] / (omega * magL), -flux * values[4] / (omega * magL), -flux * values[5] / (omega * magL));
459 }
460 /* Now pDot */
461 /* Compute the first and second terms in eq. A5 of PRD 81, 084041 (2010) */
462 for (i = 0; i < 3; i++) {
463 for (j = 0, tmpPdotT1[i] = 0.; j < 3; j++)
464 tmpPdotT1[i] += -tmpDValues[j] * Tmatrix[i][j];
465
466 tmpPdotT2[i] = -flux * values[i + 3] / (omega * magL);
467 }
468
469 /* Compute the third term in eq. A5 */
470 REAL8 tmpPdotT3T11[3][3][3], tmpPdotT3T12[3][3], tmpPdotT3T2[3];
471
472 for (i = 0; i < 3; i++)
473 for (j = 0; j < 3; j++)
474 for (l = 0; l < 3; l++)
475 for (k = 0, tmpPdotT3T11[i][j][l] = 0.; k < 3; k++)
476 tmpPdotT3T11[i][j][l] += dTijdXk[i][k][j] * invTmatrix[k][l];
477
478 if (debugPK) {
479 for (i = 0; i < 3; i++)
480 for (j = 0; j < 1; j++)
481 for (l = 0; l < 3; l++) {
482 double sum = 0;
483 for (k = 0; k < 3; k++)
484 sum += dTijdXk[i][k][j] * invTmatrix[k][l];
485 XLAL_PRINT_INFO("\n sum[%d][%d][%d] = %.12e", i, j, l, sum);
486 }
487 XLAL_PRINT_INFO("\n\n Printing dTdX * Tinverse:\n");
488 for (l = 0; l < 3; l++) {
489 for (i = 0; i < 3; i++)
490 for (j = 0; j < 3; j++) {
491 double sum = 0;
492 for (k = 0; k < 3; k++) {
493 sum += dTijdXk[i][k][l] * invTmatrix[k][j];
494 XLAL_PRINT_INFO("\n sum[%d][%d][%d] = %.12e", l, i, j, sum);
495 }
496 }
497 }
498 }
499 if (debugPK)
500 XLAL_PRINT_INFO("\npData: {%.12e, %.12e, %.12e}\n", pData[0], pData[1], pData[2]);
501 for (i = 0; i < 3; i++)
502 for (j = 0; j < 3; j++)
503 for (k = 0, tmpPdotT3T12[i][j] = 0.; k < 3; k++)
504 tmpPdotT3T12[i][j] += tmpPdotT3T11[i][j][k] * pData[k];
505
506 for (i = 0; i < 3; i++)
507 for (j = 0, tmpPdotT3T2[i] = 0.; j < 3; j++)
508 tmpPdotT3T2[i] += tmpDValues[j + 3] * Tmatrix[i][j];
509
510 for (i = 0; i < 3; i++)
511 for (j = 0, tmpPdotT3[i] = 0.; j < 3; j++)
512 tmpPdotT3[i] += tmpPdotT3T12[i][j] * tmpPdotT3T2[j];
513
514 /* Add them to obtain pDot */
515 for (i = 0; i < 3; i++)
516 dvalues[i + 3] = tmpPdotT1[i] + tmpPdotT2[i] + tmpPdotT3[i];
517
518 if (debugPK) {
519 XLAL_PRINT_INFO("\ntmpPdotT3 = ");
520 for (i = 0; i < 3; i++)
521 XLAL_PRINT_INFO("%.12e ", tmpPdotT3[i]);
522 XLAL_PRINT_INFO("\n");
523 }
524
525 /* Eqs. 11c-11d of PRD 89, 084006 (2014) */
526 /* spin1 */
527 if (debugPK) {
528 XLAL_PRINT_INFO("Raw spin1 derivatives = %.12e %.12e %.12e\n", tmpDValues[6], tmpDValues[7], tmpDValues[8]);
529 XLAL_PRINT_INFO("Raw spin2 derivatives = %.12e %.12e %.12e\n", tmpDValues[9], tmpDValues[10], tmpDValues[11]);
530 }
531 /* The factor eta is there becasue Hreal is normalized by eta */
532 dvalues[6] = eta * (tmpDValues[7] * values[8] - tmpDValues[8] * values[7]);
533 dvalues[7] = eta * (tmpDValues[8] * values[6] - tmpDValues[6] * values[8]);
534 dvalues[8] = eta * (tmpDValues[6] * values[7] - tmpDValues[7] * values[6]);
535
536 /* spin2 */
537 /* The factor eta is there becasue Hreal is normalized by eta */
538 dvalues[9] = eta * (tmpDValues[10] * values[11] - tmpDValues[11] * values[10]);
539 dvalues[10] = eta * (tmpDValues[11] * values[9] - tmpDValues[9] * values[11]);
540 dvalues[11] = eta * (tmpDValues[9] * values[10] - tmpDValues[10] * values[9]);
541
542 /* phase and precessing bit */
543 dLx = dvalues[1] * values[5] - dvalues[2] * values[4]
544 + values[1] * dvalues[5] - values[2] * dvalues[4];
545
546 dLy = dvalues[2] * values[3] - dvalues[0] * values[5]
547 + values[2] * dvalues[3] - values[0] * dvalues[5];
548
549 dLz = dvalues[0] * values[4] - dvalues[1] * values[3]
550 + values[0] * dvalues[4] - values[1] * dvalues[3];
551
552 dMagL = (Lx * dLx + Ly * dLy + Lz * dLz) / magL;
553
554 dLhatx = (dLx * magL - Lx * dMagL) / (magL * magL);
555 dLhaty = (dLy * magL - Ly * dMagL) / (magL * magL);
556
557 /*
558 * Finn Chernoff convention is used here.
559 */
560 /* Eqs. 19-20 of PRD 89, 084006 (2014) */
561 if (Lhatx == 0.0 && Lhaty == 0.0) {
562 alphadotcosi = 0.0;
563 } else {
564 alphadotcosi = Lhatz * (Lhatx * dLhaty - Lhaty * dLhatx) / (Lhatx * Lhatx + Lhaty * Lhaty);
565 }
566 /* These are ODEs for the phase that enters the h_{lm}: see Eq. 3.11 of PRD 79, 104023 (2009) */
567 dvalues[12] = omega - alphadotcosi;
568 dvalues[13] = alphadotcosi;
569 if (debugPK) {
570 XLAL_PRINT_INFO("\nIn XLALSpinPrecHcapExactDerivative:\n");
571 /* Print out all mass parameters */
572 XLAL_PRINT_INFO("m1 = %12.12lf, m2 = %12.12lf, eta = %12.12lf\n",
573 (double)mass1, (double)mass2, (double)eta);
574 /* Print out all spin parameters */
575 XLAL_PRINT_INFO("spin1 = {%12.12lf,%12.12lf,%12.12lf}, spin2 = {%12.12lf,%12.12lf,%12.12lf}\n",
576 (double)s1.data[0], (double)s1.data[1], (double)s1.data[2],
577 (double)s2.data[0], (double)s2.data[1], (double)s2.data[2]);
578 XLAL_PRINT_INFO("sigmaStar = {%12.12lf,%12.12lf,%12.12lf}, sigmaKerr = {%12.12lf,%12.12lf,%12.12lf}\n",
579 (double)sStar.data[0], (double)sStar.data[1],
580 (double)sStar.data[2], (double)sKerr.data[0],
581 (double)sKerr.data[1], (double)sKerr.data[2]);
582 XLAL_PRINT_INFO("L = {%12.12lf,%12.12lf,%12.12lf}, |L| = %12.12lf\n",
583 (double)Lx, (double)Ly, (double)Lz, (double)magL);
584 XLAL_PRINT_INFO("dLdt = {%12.12lf,%12.12lf,%12.12lf}, d|L|dt = %12.12lf\n",
585 (double)dLx, (double)dLy, (double)dLz, (double)dMagL);
586 XLAL_PRINT_INFO("Polar coordinates = {%12.12lf, %12.12lf, %12.12lf, %12.12lf}\n",
587 (double)polData[0], (double)polData[1], (double)polData[2],
588 (double)polData[3]);
589
590 XLAL_PRINT_INFO("Cartesian coordinates: {%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf}\n",
591 (double)values[0], (double)values[1], (double)values[2],
592 (double)values[3], (double)values[4], (double)values[5],
593 (double)values[6], (double)values[7], (double)values[8],
594 (double)values[9], (double)values[10], (double)values[11],
595 (double)values[12], (double)values[13]);
596 XLAL_PRINT_INFO("Cartesian derivatives: {%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf,%12.12lf}\n",
597 (double)dvalues[0], (double)dvalues[1], (double)dvalues[2],
598 (double)dvalues[3], (double)dvalues[4], (double)dvalues[5],
599 (double)dvalues[6], (double)dvalues[7], (double)dvalues[8],
600 (double)dvalues[9], (double)dvalues[10], (double)dvalues[11],
601 (double)dvalues[12], (double)dvalues[13]);
602
603 XLAL_PRINT_INFO("Hamiltonian = %12.12lf, Flux = %12.12lf, Omega = %12.12lf\n",
604 H/ (mass1 + mass2), eta*flux, omega);
605 fflush(NULL);
606 }
607
608 if(debugPK){
609 for(i=0; i<14; i++)
610 if(dvalues[i] > 1e3){
611 XLAL_PRINT_INFO("\nIn XLALSpinPrecHcapExactDerivative:\n");
612 XLAL_PRINT_INFO("Derivatives have blown up!\n");
613 for(j=0; j<14; XLAL_PRINT_INFO("dvalues[%d] = %3.12f\n", j, dvalues[j]), j++);
614 XLAL_PRINT_INFO("Flux = %3.12f\n\n", flux);
615 break;
616 }
617 }
618 return XLAL_SUCCESS;
619}
620
621
622
623
624#endif // _LALSIMIMRSPINPRECEOBHCAPEXACTDERIVATIVEPREC_V3OPT_C
#define KRONECKER_DELTA(i, j)
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 int XLALSpinPrecHcapExactDerivative(double UNUSED t, const REAL8 values[], REAL8 dvalues[], void *funcParams)
Function to calculate numerical derivatives of the spin EOB Hamiltonian, which correspond to time der...
static INT4 XLALSEOBNRv3_opt_ComputeHamiltonianDerivatives(const REAL8 *valuestort1, const REAL8 *valuestort2, SpinEOBHCoeffs *coeffs, REAL8 *derivs, SpinEOBParams *params, REAL8 *HReal)
static REAL8 sum(REAL8 *arr1, REAL8 *arr2, int k)
int l
Definition: bh_qnmode.c:135
double i
Definition: bh_ringdown.c:118
const double prT
const double H
const double deltaU
const double deltaR
const double u
const double w2
const double logu
const double u3
const double bulk
const double u5
const double logarg
const double a2
const double m1PlusetaKK
const double u2
const double u4
const double csi
const double logTerms
double REAL8
uint32_t UINT4
int32_t INT4
static const INT4 r
static const INT4 a
#define XLAL_PRINT_INFO(...)
#define XLAL_ERROR(...)
XLAL_SUCCESS
XLAL_EFUNC
The coefficients which are used in calculating the non-quasicircular correction to the EOBNRv2 model.
REAL8 * data
Parameters for the spinning EOB model, used in calculating the Hamiltonian.
Parameters for the spinning EOB model.
Definition: burst.c:245
LALDict * params
Definition: inspiral.c:248
double deltaT
Definition: unicorn.c:24