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
LALSimIMRSpinEOBInitialConditionsPrec.c
Go to the documentation of this file.
1#include <lal/LALSimInspiral.h>
2#include <lal/LALSimIMR.h>
3
4#include <gsl/gsl_vector.h>
5#include <gsl/gsl_matrix.h>
6#include <gsl/gsl_blas.h>
7
8#include <gsl/gsl_multiroots.h>
9#include <gsl/gsl_deriv.h>
10
11#include "LALSimIMRSpinEOB.h"
17
19
20#ifndef _LALSIMIMRSPINPRECEOBINITIALCONDITIONS_C
21#define _LALSIMIMRSPINPRECEOBINITIALCONDITIONS_C
22
23/* The following values are used to scale the inputs to the equations being
24 * solved to get the initial spherical orbit. By scaling we bring different
25 * inputs to the same scale, and drastically increase the rate of convergence.
26 * */
27static REAL8 scale1 = 1, scale2 = 2, scale3 = 200;
28
29/**
30 * Calculate the dot product of two vectors
31 */
32static inline
33 REAL8
35 const REAL8 a[], /**<< first vector */
36 const REAL8 b[] /**<< second vector */)
37{
38 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
39}
40
41
42/**
43 * Calculates the ith component of the cross product of two vectors,
44 * where i is in the range 0-2.
45 */
46static inline
47 REAL8
49 const INT4 i, /**<< i-th component of cross product */
50 const REAL8 a[], /**<< first vector */
51 const REAL8 b[] /**<< second vector */)
52{
53 return a[(i + 1) % 3] * b[(i + 2) % 3] - a[(i + 2) % 3] * b[(i + 1) % 3];
54}
55
56
57/**
58 * Normalizes the given vector
59 */
60static inline
61int
63{
64 REAL8 norm = sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
65
66 a[0] /= norm;
67 a[1] /= norm;
68 a[2] /= norm;
69
70 return XLAL_SUCCESS;
71}
72
73/**
74 * Calculate the rotation matrix and its inverse.
75 * Rotate the ex-ey-ez frame to the r-v-L frame.
76 * This static function is called only by XLALSimIMRSpinEOBInitialConditionsPrec
77 */
78static int
80 gsl_matrix * rotMatrix, /**< OUTPUT, rotation matrix */
81 gsl_matrix * rotInverse, /**< OUTPUT, rotation matrix inversed */
82 REAL8 r[], /**< position vector */
83 REAL8 v[], /**< velocity vector */
84 REAL8 L[] /**< orbital angular momentum */
85)
86{
87
88 /** a, b, g are the angles alpha, beta and gamma */
89 /* Use a, b and g to avoid shadowing gamma and getting a warning */
90 REAL8 a , b, g;
91 REAL8 cosa , sina, cosb, sinb, cosg, sing;
92
93 /* Calculate the Euclidean Euler angles */
94
95 /* We need to avoid a singularity if L is along z axis */
96 if (L[2] > 0.9999) {
97 a = b = g = 0.0;
98 } else {
99 a = atan2(L[0], -L[1]);
100 b = atan2(sqrt(L[0] * L[0] + L[1] * L[1]), L[2]);
101 g = atan2(r[2], v[2]);
102 }
103
104 if ((cosa = cos(a)) < 1.0e-16)
105 cosa = 0.0;
106 if ((sina = sin(a)) < 1.0e-16)
107 sina = 0.0;
108 if ((cosb = cos(b)) < 1.0e-16)
109 cosb = 0.0;
110 if ((sinb = sin(b)) < 1.0e-16)
111 sinb = 0.0;
112 if ((cosg = cos(g)) < 1.0e-16)
113 cosg = 0.0;
114 if ((sing = sin(g)) < 1.0e-16)
115 sing = 0.0;
116
117 /**
118 * Implement the Rotation Matrix following the "x-convention"
119 * 1. rotate about the z-axis by an angle a, rotation matrix Rz(a);
120 * 2. rotate about the former x-axis (now x') by an angle b, rotation matrix Rx(b);
121 * 3. rotate about the former z-axis (now z') by an algle g, rotation matrix Rz(g);
122 */
123 /* populate the matrix */
124 /*
125 * gsl_matrix_set( rotMatrix, 0, 0, cosg*cosa - cosb*sina*sing );
126 * gsl_matrix_set( rotMatrix, 0, 1, cosg*sina + cosb*cosa*sing );
127 * gsl_matrix_set( rotMatrix, 0, 2, sing*sinb ); gsl_matrix_set(
128 * rotMatrix, 1, 0, -sing*cosa - cosb*sina*cosg ); gsl_matrix_set(
129 * rotMatrix, 1, 1, -sing*sina + cosb*cosa*cosg ); gsl_matrix_set(
130 * rotMatrix, 1, 2, cosg*sinb ); gsl_matrix_set( rotMatrix, 2, 0,
131 * sinb*sina ); gsl_matrix_set( rotMatrix, 2, 1, -sinb*cosa );
132 * gsl_matrix_set( rotMatrix, 2, 2, cosb );
133 */
134 //Better to avoid angle singularity
135 gsl_matrix_set(rotMatrix, 0, 0, r[0]);
136 gsl_matrix_set(rotMatrix, 0, 1, r[1]);
137 gsl_matrix_set(rotMatrix, 0, 2, r[2]);
138 gsl_matrix_set(rotMatrix, 1, 0, v[0]);
139 gsl_matrix_set(rotMatrix, 1, 1, v[1]);
140 gsl_matrix_set(rotMatrix, 1, 2, v[2]);
141 gsl_matrix_set(rotMatrix, 2, 0, L[0]);
142 gsl_matrix_set(rotMatrix, 2, 1, L[1]);
143 gsl_matrix_set(rotMatrix, 2, 2, L[2]);
144
145
146 /* Now populate the transpose (which should be the inverse) */
147 gsl_matrix_transpose_memcpy(rotInverse, rotMatrix);
148
149 /* Test that the code does what it should do */
150
151 /*
152 * gsl_matrix *ab = gsl_matrix_alloc( 3, 3 ); gsl_blas_dgemm(
153 * CblasNoTrans, CblasNoTrans, 1., rotMatrix, rotInverse, 0., ab );
154 */
155
156 /*
157 * XLAL_PRINT_INFO( "The generated rotation matrix: \n" ); for ( int i = 0; i
158 * < 3; i++ ) { for ( int j = 0; j < 3; j++ ) { XLAL_PRINT_INFO( "%.16e\t",
159 * gsl_matrix_get( rotMatrix, i, j ) ); } XLAL_PRINT_INFO( "\n" ); }
160 */
161
162 /*
163 * XLAL_PRINT_INFO( "Is the transpose of the rotation matrix the inverse?\n"
164 * ); for ( int i = 0; i < 3; i++ ) { for ( int j = 0; j < 3; j++ ) {
165 * XLAL_PRINT_INFO( "%.16e\t", gsl_matrix_get( ab, i, j ) ); } XLAL_PRINT_INFO( "\n" );
166 * }
167 */
168 /* gsl_matrix_free( ab ); */
169
170 return XLAL_SUCCESS;
171}
172
173
174/**
175 * Applies a rotation matrix to a given vector
176 */
177static int
179 gsl_matrix * rotMatrix, /**< rotation matrix */
180 REAL8 a[] /**< OUTPUT, vector rotated */
181)
182{
183
184 double b[3];
185 gsl_vector_view tmpView1 = gsl_vector_view_array(a, 3);
186 gsl_vector_view tmpView2 = gsl_vector_view_array(b, 3);
187 gsl_vector *tmpVec1 = &tmpView1.vector;
188 gsl_vector *tmpVec2 = &tmpView2.vector;
189
190 gsl_blas_dgemv(CblasNoTrans, 1.0, rotMatrix, tmpVec1, 0.0, tmpVec2);
191
192 gsl_vector_memcpy(tmpVec1, tmpVec2);
193
194 return XLAL_SUCCESS;
195}
196
197/**
198 * Performs a co-ordinate transformation from spherical to Cartesian co-ordinates.
199 * In the code from Tyson Littenberg, this was only applicable to the special theta=pi/2, phi=0 case.
200 * This special transformation is a static function called only by
201 * GSLSpinHamiltonianDerivWrapperPrec and XLALSimIMRSpinEOBInitialConditionsPrec
202 */
203static int
205 REAL8 qCart[], /**<< OUTPUT, position vector in Cartesean coordinates */
206 REAL8 pCart[], /**<< OUTPUT, momentum vector in Cartesean coordinates */
207 const REAL8 qSph[], /**<< position vector in spherical coordinates */
208 const REAL8 pSph[] /**<< momentum vector in spherical coordinates */
209)
210{
211
212 REAL8 r;
213 REAL8 pr , pTheta, pPhi;
214
215 REAL8 x , y, z;
216 REAL8 pX , pY, pZ;
217
218 r = qSph[0];
219 pr = pSph[0];
220 pTheta = pSph[1];
221 pPhi = pSph[2];
222
223 x = r;
224 y = 0.0;
225 z = 0.0;
226
227 pX = pr;
228 pY = pPhi / r;
229 pZ = -pTheta / r;
230
231 /* Copy these values to the output vectors */
232 qCart[0] = x;
233 qCart[1] = y;
234 qCart[2] = z;
235 pCart[0] = pX;
236 pCart[1] = pY;
237 pCart[2] = pZ;
238
239 return XLAL_SUCCESS;
240}
241
242
243/**
244 * Perform a co-ordinate transformation from Cartesian to spherical co-ordinates.
245 * In the code from Tyson Littenberg, this was only applicable to the special theta=pi/2, phi=0 case.
246 * This special transformation is a static function called only by XLALSimIMRSpinEOBInitialConditionsPrec
247 */
248static int
250 REAL8 qSph[], /**<< OUTPUT, position vector in spherical coordinates */
251 REAL8 pSph[], /**<< OUTPUT, momentum vector in Cartesean coordinates */
252 const REAL8 qCart[], /**<< position vector in spherical coordinates */
253 const REAL8 pCart[] /**<< momentum vector in Cartesean coordinates */
254)
255{
256
257 REAL8 r;
258 REAL8 pr , pTheta, pPhi;
259
260 REAL8 x;
261 //, y, z;
262 REAL8 pX , pY, pZ;
263
264 x = qCart[0];
265 //y = qCart[1];
266 //z = qCart[2];
267 pX = pCart[0];
268 pY = pCart[1];
269 pZ = pCart[2];
270
271 r = x;
272 pr = pX;
273 pTheta = -r * pZ;
274 pPhi = r * pY;
275
276 /* Fill the output vectors */
277 qSph[0] = r;
278 qSph[1] = LAL_PI_2;
279 qSph[2] = 0.0;
280 pSph[0] = pr;
281 pSph[1] = pTheta;
282 pSph[2] = pPhi;
283
284 return XLAL_SUCCESS;
285}
286
287
288/**
289 * Root function for gsl root finders.
290 * The gsl root finder will look for gsl_vector *x position in parameter space
291 * where the returned values in gsl_vector *f are zero.
292 * The functions on which we search for roots are:
293 * dH/dr, dH/dptheta and dH/dpphi-omega,
294 * namely, Eqs. 4.8 and 4.9 of Buonanno et al. PRD 74, 104005 (2006).
295 */
296static int
298 const gsl_vector * x, /**<< Parameters requested by gsl root finder */
299 void *params, /**<< Spin EOB parameters */
300 gsl_vector * f /**<< Function values for the given parameters*/
301)
302{
303 int debugPK = 0;
304 SEOBRootParams *rootParams = (SEOBRootParams *) params;
305 REAL8 mTotal = rootParams->params->eobParams->m1 + rootParams->params->eobParams->m2;
306 if (debugPK)
307 XLAL_PRINT_INFO("\n\nmtotal in XLALFindSphericalOrbitPrec = %f\n", (float)mTotal);
308
309 REAL8 py , pz, r, ptheta, pphi;
310
311 /* Numerical derivative of Hamiltonian wrt given value */
312 REAL8 dHdx , dHdpy, dHdpz;
313 REAL8 dHdr , dHdptheta, dHdpphi;
314
315 /* Populate the appropriate values */
316 /* In the special theta=pi/2 phi=0 case, r is x */
317
318 REAL8 temp = gsl_vector_get(x, 0)/scale1;
319 REAL8 prefactor = 1.0;
320 if (temp < 0.0){
321 prefactor=-1.0;
322 }
323
324 rootParams->values[0] = r = sqrt(temp*temp+36.0);
325 rootParams->values[4] = py = gsl_vector_get(x, 1)/scale2;
326 rootParams->values[5] = pz = gsl_vector_get(x, 2)/scale3;
327 //printf("r is %.17f\n",r);
328 if(isnan(rootParams->values[0])) {
329 rootParams->values[0] = 100.;
330 }
331 if(isnan(rootParams->values[4])) {
332 rootParams->values[4] = 0.1;
333 }
334 if(isnan(rootParams->values[5])) {
335 rootParams->values[5] = 0.01;
336 }
337 if (debugPK) {
338 XLAL_PRINT_INFO("%3.10f %3.10f %3.10f %3.10f %3.10f %3.10f\n",
339 rootParams->values[0], rootParams->values[1], rootParams->values[2],
340 rootParams->values[3], rootParams->values[4], rootParams->values[5]);
341 XLAL_PRINT_INFO("Values r = %.16e, py = %.16e, pz = %.16e\n", r, py, pz);
342 fflush(NULL);
343 }
344
345 ptheta = -r * pz;
346 pphi = r * py;
347
348 if (debugPK)
349 XLAL_PRINT_INFO("Input Values r = %.16e, py = %.16e, pz = %.16e\n pthetha = %.16e pphi = %.16e\n", r, py, pz, ptheta, pphi);
350
351 /* dH by dR and dP */
352 REAL8 tmpDValues[14];
353 int UNUSED status;
354 for (int i = 0; i < 3; i++) {
355 rootParams->values[i + 6] /= mTotal * mTotal;
356 rootParams->values[i + 9] /= mTotal * mTotal;
357 }
358 UINT4 oldignoreflux = rootParams->params->ignoreflux;
359 rootParams->params->ignoreflux = 1;
360 status = XLALSpinPrecHcapNumericalDerivative(0, rootParams->values, tmpDValues, rootParams->params);
361 rootParams->params->ignoreflux = oldignoreflux;
362 for (int i = 0; i < 3; i++) {
363 rootParams->values[i + 6] *= mTotal * mTotal;
364 rootParams->values[i + 9] *= mTotal * mTotal;
365 }
366 REAL8 rvec[3] =
367 {rootParams->values[0], rootParams->values[1], rootParams->values[2]};
368 REAL8 pvec[3] =
369 {rootParams->values[3], rootParams->values[4], rootParams->values[5]};
370 REAL8 chi1vec[3] =
371 {rootParams->values[6], rootParams->values[7], rootParams->values[8]};
372 REAL8 chi2vec[3] =
373 {rootParams->values[9], rootParams->values[10], rootParams->values[11]};
374 REAL8 Lvec[3] = {CalculateCrossProductPrec(0, rvec, pvec),
375 CalculateCrossProductPrec(1, rvec, pvec), CalculateCrossProductPrec(2, rvec, pvec)};
376 REAL8 theta1 = acos(CalculateDotProductPrec(chi1vec, Lvec)
377 / sqrt(CalculateDotProductPrec(chi1vec, chi1vec))
378 / sqrt(CalculateDotProductPrec(Lvec, Lvec)));
379 REAL8 theta2 = acos(CalculateDotProductPrec(chi2vec, Lvec)
380 / sqrt(CalculateDotProductPrec(chi2vec, chi2vec))
381 / sqrt(CalculateDotProductPrec(Lvec, Lvec)));
382
383 if (debugPK) {
384 XLAL_PRINT_INFO("rvec = %.16e %.16e %.16e\n", rvec[0], rvec[1], rvec[2]);
385 XLAL_PRINT_INFO("pvec = %.16e %.16e %.16e\n", pvec[0], pvec[1], pvec[2]);
386 XLAL_PRINT_INFO("theta1 = %.16e\n", theta1);
387 XLAL_PRINT_INFO("theta2 = %.16e\n", theta2);
388 }
389
390 if (theta1 > 1.0e-6 && theta2 >= 1.0e-6) {
391 dHdx = -tmpDValues[3];
392 dHdpy = tmpDValues[1];
393 dHdpz = tmpDValues[2];
394 } else {
395 //rootParams->values[5] = 0.;
396 //rootParams->values[6] = 0.;
397 //rootParams->values[7] = 0.;
398 //rootParams->values[8] = sqrt(CalculateDotProductPrec(chi1vec, chi1vec));
399 //rootParams->values[9] = 0.;
400 //rootParams->values[10]= 0.;
401 //rootParams->values[11]= sqrt(CalculateDotProductPrec(chi2vec, chi2vec));
403 rootParams->values, rootParams->params);
405 rootParams->values, rootParams->params);
407 rootParams->values, rootParams->params);
408 }
412 if (debugPK)
413 XLAL_PRINT_INFO("dHdx = %.16e, dHdpy = %.16e, dHdpz = %.16e\n", dHdx, dHdpy, dHdpz);
414
415 /* Now convert to spherical polars */
416 dHdr = dHdx - dHdpy * pphi / (r * r) + dHdpz * ptheta / (r * r);
417 dHdptheta = -dHdpz / r;
418 dHdpphi = dHdpy / r;
419
420 if (debugPK)
421 XLAL_PRINT_INFO("dHdr = %.16e dHdptheta = %.16e dHdpphi = %.16e\n",
422 dHdr, dHdptheta, dHdpphi);
423 /* populate the function vector */
424 gsl_vector_set(f, 0, dHdr);
425 gsl_vector_set(f, 1, dHdptheta);
426 gsl_vector_set(f, 2, dHdpphi - rootParams->omega);
427
428 if (debugPK)
429 XLAL_PRINT_INFO("Current funcvals = %.16e %.16e %.16e\n",
430 gsl_vector_get(f, 0), gsl_vector_get(f, 1), gsl_vector_get(f, 2) );
431
432 /* Rescale back */
433
434
435 rootParams->values[0] = prefactor*scale1*(sqrt(rootParams->values[0]*rootParams->values[0]-36.0));
436 rootParams->values[4] *= scale2;
437 rootParams->values[5] *= scale3;
438
439 return XLAL_SUCCESS;
440}
441
442
443/**
444 * Wrapper for calculating specific derivative of the Hamiltonian in spherical co-ordinates,
445 * dH/dr, dH/dptheta and dH/dpphi.
446 * It only works for the specific co-ord system we use here
447 */
448static double
449GSLSpinHamiltonianDerivWrapperPrec(double x, /**<< Derivative at x */
450 void *params /**<< Function parameters */
451)
452{
453 int debugPK = 0;
455 REAL8 mTotal = dParams->params->eobParams->m1 + dParams->params->eobParams->m2;
456 REAL8 sphValues[12];
457 REAL8 cartValues[12];
458
459 REAL8 dHdr , dHdx, dHdpy, dHdpz;
460 REAL8 r , ptheta, pphi;
461
462 memcpy(sphValues, dParams->sphValues, sizeof(sphValues));
463 sphValues[dParams->varyParam1] = x;
464
465 SphericalToCartesianPrec(cartValues, cartValues + 3, sphValues, sphValues + 3);
466 memcpy(cartValues + 6, sphValues + 6, 6 * sizeof(REAL8));
467
468 r = sphValues[0];
469 ptheta = sphValues[4];
470 pphi = sphValues[5];
471
472 /* New code to compute derivatives w.r.t. cartesian variables */
473 REAL8 tmpDValues[14];
474 int UNUSED status;
475 for (int i = 0; i < 3; i++) {
476 cartValues[i + 6] /= mTotal * mTotal;
477 cartValues[i + 9] /= mTotal * mTotal;
478 }
479 UINT4 oldignoreflux = dParams->params->ignoreflux;
480 dParams->params->ignoreflux = 1;
481 status = XLALSpinPrecHcapNumericalDerivative(0, cartValues, tmpDValues, dParams->params);
482 dParams->params->ignoreflux = oldignoreflux;
483 for (int i = 0; i < 3; i++) {
484 cartValues[i + 6] *= mTotal * mTotal;
485 cartValues[i + 9] *= mTotal * mTotal;
486 }
487
488 double rvec [3] = {cartValues[0], cartValues[1], cartValues[2]};
489 double pvec [3] = {cartValues[3], cartValues[4], cartValues[5]};
490 double chi1vec [3] = {cartValues[6], cartValues[7], cartValues[8]};
491 double chi2vec [3] = {cartValues[9], cartValues[10], cartValues[11]};
492 double Lvec [3] = {CalculateCrossProductPrec(0, rvec, pvec), CalculateCrossProductPrec(1, rvec, pvec), CalculateCrossProductPrec(2, rvec, pvec)};
493 double theta1 = acos(CalculateDotProductPrec(chi1vec, Lvec) / sqrt(CalculateDotProductPrec(chi1vec, chi1vec)) / sqrt(CalculateDotProductPrec(Lvec, Lvec)));
494 double theta2 = acos(CalculateDotProductPrec(chi2vec, Lvec) / sqrt(CalculateDotProductPrec(chi2vec, chi2vec)) / sqrt(CalculateDotProductPrec(Lvec, Lvec)));
495 if (debugPK) {
496 XLAL_PRINT_INFO("In GSLSpinHamiltonianDerivWrapperPrec");
497 XLAL_PRINT_INFO("rvec = %.16e %.16e %.16e\n", rvec[0], rvec[1], rvec[2]);
498 XLAL_PRINT_INFO("pvec = %.16e %.16e %.16e\n", pvec[0], pvec[1], pvec[2]);
499 XLAL_PRINT_INFO("theta1 = %.16e\n", theta1);
500 XLAL_PRINT_INFO("theta2 = %.16e\n", theta2);
501 }
502 if (theta1 > 1.0e-6 && theta2 >= 1.0e-6) {
503 switch (dParams->varyParam2) {
504 case 0:
505 /* dHdr */
506 dHdx = -tmpDValues[3];
507 //XLALSpinPrecHcapNumDerivWRTParam(0, cartValues, dParams->params);
508 dHdpy = tmpDValues[1];
509 //XLALSpinPrecHcapNumDerivWRTParam(4, cartValues, dParams->params);
510 dHdpz = tmpDValues[2];
511 //XLALSpinPrecHcapNumDerivWRTParam(5, cartValues, dParams->params);
512
513 dHdr = dHdx - dHdpy * pphi / (r * r) + dHdpz * ptheta / (r * r);
514 //XLAL_PRINT_INFO("dHdr = %.16e\n", dHdr);
515 return dHdr;
516
517 break;
518 case 4:
519 /* dHdptheta */
520 dHdpz = tmpDValues[2];
521 //XLALSpinPrecHcapNumDerivWRTParam(5, cartValues, dParams->params);
522 return -dHdpz / r;
523 break;
524 case 5:
525 /* dHdpphi */
526 dHdpy = tmpDValues[1];
527 //XLALSpinPrecHcapNumDerivWRTParam(4, cartValues, dParams->params);
528 return dHdpy / r;
529 break;
530 default:
531 XLALPrintError("This option is not supported in the second derivative function!\n");
533 break;
534 }
535 } else {
536 switch (dParams->varyParam2) {
537 case 0:
538 /* dHdr */
539 dHdx = XLALSpinPrecHcapNumDerivWRTParam(0, cartValues, dParams->params);
540 dHdpy = XLALSpinPrecHcapNumDerivWRTParam(4, cartValues, dParams->params);
541 dHdpz = XLALSpinPrecHcapNumDerivWRTParam(5, cartValues, dParams->params);
542 dHdr = dHdx - dHdpy * pphi / (r * r) + dHdpz * ptheta / (r * r);
543 //XLAL_PRINT_INFO("dHdr = %.16e\n", dHdr);
544 return dHdr;
545
546 break;
547 case 4:
548 /* dHdptheta */
549 dHdpz = XLALSpinPrecHcapNumDerivWRTParam(5, cartValues, dParams->params);
550 return -dHdpz / r;
551 break;
552 case 5:
553 /* dHdpphi */
554 dHdpy = XLALSpinPrecHcapNumDerivWRTParam(4, cartValues, dParams->params);
555 return dHdpy / r;
556 break;
557 default:
558 XLALPrintError("This option is not supported in the second derivative function!\n");
560 break;
561 }
562 }
563}
564
565static int XLALRobustDerivative(const gsl_function * F, double x, double h,
566 double *result, double *absErr){
567 // Take the derivative
568 gsl_deriv_central(F, x,h, result, absErr);
569 // We check the estimate of the error
570 REAL8 frac= 0.01; // Adjust this to change how accurtate we demand the derivative to be
571 if (fabs(*absErr)>frac*fabs(*result)){
572 REAL8 temp1 = 0.0;
573 REAL8 temp2 =0.0;
574 REAL8 absErr1 = 0.0;
575 REAL8 absErr2 = 0.0;
576 UINT4 deriv_ok = 0;
577 UINT4 n = 1;
578 while (!deriv_ok && n<=10){
579 XLAL_PRINT_WARNING("Warning: second derivative computation went wrong. Trying again");
580 gsl_deriv_central(F, x,2*n*h, &temp1, &absErr1);
581 gsl_deriv_central(F, x,h/(2*n), &temp2, &absErr2);
582 if (fabs(absErr1)/fabs(temp1)<fabs(absErr2)/fabs(temp2) && fabs(absErr1)<frac*fabs(temp1)){
583 *result = temp1;
584 *absErr = absErr1;
585 deriv_ok = 1;
586 break;
587 }
588 else if (fabs(absErr1)/fabs(temp1)>fabs(absErr2)/fabs(temp2) && fabs(absErr2)<frac*fabs(temp2)){
589 *result = temp2;
590 *absErr = absErr2;
591 deriv_ok = 1;
592 break;
593 }
594 else{
595 n+=1;
596 }
597 }
598 if (!deriv_ok){
599 XLAL_PRINT_ERROR("The computation of the second derivative of H in initial data has failed!");
601 }
602 }
603 return XLAL_SUCCESS;
604}
605
606/**
607 * Function to calculate the second derivative of the Hamiltonian.
608* The derivatives are taken with respect to indices idx1, idx2 */
609static REAL8
611 const int idx1, /**<< Derivative w.r.t. index 1 */
612 const int idx2, /**<< Derivative w.r.t. index 2 */
613 const REAL8 values[], /**<< Dynamical variables in spherical coordinates */
614 SpinEOBParams * params, /**<< Spin EOB Parameters */
615 INT4 use_optimized /**<< use_optimized=1 -> use EXACT instead of finite difference derivatives... UNSUPPORTED AT MOMENT. */
616)
617{
618
619 static const REAL8 STEP_SIZE = 3.0e-3;
620
621 REAL8 result;
622 REAL8 UNUSED absErr;
623
624 HcapSphDeriv2Params dParams;
625
626 gsl_function F;
627 INT4 UNUSED gslStatus;
628
629 dParams.sphValues = values;
630 dParams.varyParam1 = idx1;
631 dParams.varyParam2 = idx2;
632 dParams.params = params;
633 dParams.use_optimized = use_optimized;
634
635 /*
636 * XLAL_PRINT_INFO( " In second deriv function: values\n" ); for ( int i = 0;
637 * i < 12; i++ ) { XLAL_PRINT_INFO( "%.16e ", values[i] ); } XLAL_PRINT_INFO( "\n" );
638 */
640 F.params = &dParams;
641
642 /* GSL seemed to give weird answers - try my own code */
643 /*
644 * result = GSLSpinHamiltonianDerivWrapperPrec( values[idx1] + STEP_SIZE,
645 * &dParams ) - GSLSpinHamiltonianDerivWrapperPrec( values[idx1] -
646 * STEP_SIZE, &dParams ); XLAL_PRINT_INFO( "%.16e - %.16e / 2h\n",
647 * GSLSpinHamiltonianDerivWrapperPrec( values[idx1] + STEP_SIZE, &dParams
648 * ), GSLSpinHamiltonianDerivWrapperPrec( values[idx1] - STEP_SIZE,
649 * &dParams ) );
650 *
651 * result = result / ( 2.*STEP_SIZE );
652 */
653
654 //XLAL_CALLGSL(gslStatus = gsl_deriv_central(&F, values[idx1],
655 // STEP_SIZE, &result, &absErr));
656 XLALRobustDerivative(&F, values[idx1],
657 STEP_SIZE, &result, &absErr);
658 /*
659 if (gslStatus != GSL_SUCCESS) {
660 XLALPrintError("XLAL Error %s - Failure in GSL function\n", __func__);
661 XLAL_ERROR_REAL8(XLAL_EDOM);
662 }
663 */
664 //XLAL_PRINT_INFO("Second deriv abs err = %.16e\n", absErr);
665
666 //XLAL_PRINT_INFO("RESULT = %.16e\n", result);
667 return result;
668}
669
670/**
671 * Main function for calculating the spinning EOB initial conditions, following the
672 * quasi-spherical initial conditions described in Sec. IV A of
673 * Buonanno, Chen & Damour PRD 74, 104005 (2006).
674 * All equation numbers in the comments of this function refer to this paper.
675 *
676 * Inputs are binary parameters (masses, spin vectors and inclination),
677 * EOB model parameters and initial frequency.
678 *
679 * Outputs are initial dynamical variables in a vector
680 * (x, y, z, px, py, pz, S1x, S1y, S1z, S2x, S2y, S2z).
681 *
682 * In the paper, the initial conditions are solved for a given initial radius,
683 * while in this function, they are solved for a given inital orbital frequency.
684 *
685 * This difference is reflected in solving Eq. (4.8).
686 * The calculation is done in 5 steps:
687 * STEP 1) Rotate to LNhat0 along z-axis and N0 along x-axis frame, where
688 * LNhat0 and N0 are initial normal to orbital plane and initial orbital separation;
689 * STEP 2) After rotation in STEP 1, in spherical coordinates,
690 * phi0 and theta0 are given directly in Eq. (4.7),
691 * r0, pr0, ptheta0 and pphi0 are obtained by solving Eqs. (4.8) and (4.9)
692 * (using gsl_multiroot_fsolver).
693 * At this step, we find initial conditions for a spherical orbit without
694 * radiation reaction.
695 * STEP 3) Rotate to L0 along z-axis and N0 along x-axis frame, where
696 * L0 is the initial orbital angular momentum and
697 * L0 is calculated using initial position and linear momentum obtained in STEP 2.
698 * STEP 4) In the L0-N0 frame, we calculate (dE/dr)|sph using Eq. (4.14),
699 * then initial dr/dt using Eq. (4.10), and finally pr0 using Eq. (4.15).
700 * STEP 5) Rotate back to the original inertial frame by inverting the rotation of STEP 3 and
701 * then inverting the rotation of STEP 1.
702 */
703
704static int
706 REAL8Vector * initConds, /**<< OUTPUT, Initial dynamical variables */
707 const REAL8 mass1, /**<< mass 1 */
708 const REAL8 mass2, /**<< mass 2 */
709 const REAL8 fMin, /**<< Initial frequency (given) */
710 const REAL8 inc, /**<< Inclination */
711 const REAL8 spin1[], /**<< Initial spin vector 1 */
712 const REAL8 spin2[], /**<< Initial spin vector 2 */
713 SpinEOBParams * params, /**<< Spin EOB parameters */
714 INT4 use_optimized /**<< use_optimized=1 -> use EXACT instead of finite difference derivatives, where needed */
715)
716{
717
718 if (!initConds) {
720 }
721
722 int debugPK = 0; int printPK = 0;
723 FILE* UNUSED out = NULL;
724
725 if (printPK) {
726 XLAL_PRINT_INFO("Inside the XLALSimIMRSpinEOBInitialConditionsPrec function!\n");
728 "Inputs: m1 = %.16e, m2 = %.16e, fMin = %.16e, inclination = %.16e\n",
729 mass1, mass2, (double)fMin, (double)inc);
730 XLAL_PRINT_INFO("Inputs: mSpin1 = {%.16e, %.16e, %.16e}\n",
731 spin1[0], spin1[1], spin1[2]);
732 XLAL_PRINT_INFO("Inputs: mSpin2 = {%.16e, %.16e, %.16e}\n",
733 spin2[0], spin2[1], spin2[2]);
734 fflush(NULL);
735 }
736 static const int UNUSED lMax = 8;
737
738 int i;
739
740 /* Variable to keep track of whether the user requested the tortoise */
741 int tmpTortoise;
742
743 UINT4 SpinAlignedEOBversion;
744
745 REAL8 mTotal;
746 REAL8 eta;
747 REAL8 omega , v0; /* Initial velocity and angular
748 * frequency */
749
750 REAL8 ham; /* Hamiltonian */
751
752 REAL8 LnHat [3]; /* Initial orientation of angular
753 * momentum */
754 REAL8 rHat [3]; /* Initial orientation of radial
755 * vector */
756 REAL8 vHat [3]; /* Initial orientation of velocity
757 * vector */
758 REAL8 Lhat [3]; /* Direction of relativistic ang mom */
759 REAL8 qHat [3];
760 REAL8 pHat [3];
761
762 /* q and p vectors in Cartesian and spherical coords */
763 REAL8 qCart [3], pCart[3];
764 REAL8 qSph [3], pSph[3];
765
766 /* We will need to manipulate the spin vectors */
767 /* We will use temporary vectors to do this */
768 REAL8 tmpS1 [3];
769 REAL8 tmpS2 [3];
770 REAL8 tmpS1Norm[3];
771 REAL8 tmpS2Norm[3];
772
773 REAL8Vector qCartVec, pCartVec;
774 REAL8Vector s1Vec, s2Vec, s1VecNorm, s2VecNorm;
775 REAL8Vector sKerr, sStar;
776 REAL8 sKerrData[3], sStarData[3];
777 REAL8 a = 0.;
778 //, chiS, chiA;
779 //REAL8 chi1, chi2;
780
781 /*
782 * We will need a full values vector for calculating derivs of
783 * Hamiltonian
784 */
785 REAL8 sphValues[12];
786 REAL8 cartValues[12];
787
788 /* Matrices for rotating to the new basis set. */
789 /* It is more convenient to calculate the ICs in a simpler basis */
790 gsl_matrix *rotMatrix = NULL;
791 gsl_matrix *invMatrix = NULL;
792 gsl_matrix *rotMatrix2 = NULL;
793 gsl_matrix *invMatrix2 = NULL;
794
795 /* Root finding stuff for finding the spherical orbit */
796 SEOBRootParams rootParams;
797 const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids;
798 gsl_multiroot_fsolver *rootSolver = NULL;
799
800 gsl_multiroot_function rootFunction;
801 gsl_vector *initValues = NULL;
802 gsl_vector *finalValues = NULL;
803 INT4 gslStatus;
804 INT4 cntGslNoProgress = 0, MAXcntGslNoProgress = 5;
805 //INT4 cntGslNoProgress = 0, MAXcntGslNoProgress = 50;
806 REAL8 multFacGslNoProgress = 3./5.;
807 //const int maxIter = 2000;
808 const int maxIter = 10000;
809
810 memset(&rootParams, 0, sizeof(rootParams));
811
812 mTotal = mass1 + mass2;
813 eta = mass1 * mass2 / (mTotal * mTotal);
814 memcpy(tmpS1, spin1, sizeof(tmpS1));
815 memcpy(tmpS2, spin2, sizeof(tmpS2));
816 memcpy(tmpS1Norm, spin1, sizeof(tmpS1Norm));
817 memcpy(tmpS2Norm, spin2, sizeof(tmpS2Norm));
818 for (i = 0; i < 3; i++) {
819 tmpS1Norm[i] /= mTotal * mTotal;
820 tmpS2Norm[i] /= mTotal * mTotal;
821 }
822 SpinAlignedEOBversion = params->seobCoeffs->SpinAlignedEOBversion;
823 /* We compute the ICs for the non-tortoise p, and convert at the end */
824 tmpTortoise = params->tortoise;
825 params->tortoise = 0;
826
827 EOBNonQCCoeffs *nqcCoeffs = NULL;
828 nqcCoeffs = params->nqcCoeffs;
829
830 /*
831 * STEP 1) Rotate to LNhat0 along z-axis and N0 along x-axis frame,
832 * where LNhat0 and N0 are initial normal to orbital plane and
833 * initial orbital separation;
834 */
835
836 /* Set the initial orbital ang mom direction. Taken from STPN code */
837 LnHat[0] = sin(inc);
838 LnHat[1] = 0.;
839 LnHat[2] = cos(inc);
840
841 /*
842 * Set the radial direction - need to take care to avoid singularity
843 * if L is along z axis
844 */
845 if (LnHat[2] > 0.9999) {
846 rHat[0] = 1.;
847 rHat[1] = rHat[2] = 0.;
848 } else {
849 REAL8 theta0 = atan(-LnHat[2] / LnHat[0]); /* theta0 is between 0
850 * and Pi */
851 rHat[0] = sin(theta0);
852 rHat[1] = 0;
853 rHat[2] = cos(theta0);
854 }
855
856 /* Now we can complete the triad */
857 vHat[0] = CalculateCrossProductPrec(0, LnHat, rHat);
858 vHat[1] = CalculateCrossProductPrec(1, LnHat, rHat);
859 vHat[2] = CalculateCrossProductPrec(2, LnHat, rHat);
860
862
863 /* Vectors BEFORE rotation */
864 if (printPK) {
865 for (i = 0; i < 3; i++)
866 XLAL_PRINT_INFO(" LnHat[%d] = %.16e, rHat[%d] = %.16e, vHat[%d] = %.16e\n",
867 i, LnHat[i], i, rHat[i], i, vHat[i]);
868
869 XLAL_PRINT_INFO("\n\n");
870 for (i = 0; i < 3; i++)
871 XLAL_PRINT_INFO(" s1[%d] = %.16e, s2[%d] = %.16e\n", i, tmpS1[i], i, tmpS2[i]);
872 fflush(NULL);
873 }
874
875 /* Allocate and compute the rotation matrices */
876 XLAL_CALLGSL(rotMatrix = gsl_matrix_alloc(3, 3));
877 XLAL_CALLGSL(invMatrix = gsl_matrix_alloc(3, 3));
878 if (!rotMatrix || !invMatrix) {
879 if (rotMatrix)
880 gsl_matrix_free(rotMatrix);
881 if (invMatrix)
882 gsl_matrix_free(invMatrix);
884 }
885 if (CalculateRotationMatrixPrec(rotMatrix, invMatrix, rHat, vHat, LnHat) == XLAL_FAILURE) {
886 gsl_matrix_free(rotMatrix);
887 gsl_matrix_free(invMatrix);
889 }
890 /* Rotate the orbital vectors and spins */
891 ApplyRotationMatrixPrec(rotMatrix, rHat);
892 ApplyRotationMatrixPrec(rotMatrix, vHat);
893 ApplyRotationMatrixPrec(rotMatrix, LnHat);
894 ApplyRotationMatrixPrec(rotMatrix, tmpS1);
895 ApplyRotationMatrixPrec(rotMatrix, tmpS2);
896 ApplyRotationMatrixPrec(rotMatrix, tmpS1Norm);
897 ApplyRotationMatrixPrec(rotMatrix, tmpS2Norm);
898
899 /* See if Vectors have been rotated fine */
900 if (printPK) {
901 XLAL_PRINT_INFO("\nAfter applying rotation matrix:\n\n");
902 for (i = 0; i < 3; i++)
903 XLAL_PRINT_INFO(" LnHat[%d] = %.16e, rHat[%d] = %.16e, vHat[%d] = %.16e\n",
904 i, LnHat[i], i, rHat[i], i, vHat[i]);
905
906 XLAL_PRINT_INFO("\n");
907 for (i = 0; i < 3; i++)
908 XLAL_PRINT_INFO(" s1[%d] = %.16e, s2[%d] = %.16e\n", i, tmpS1[i], i, tmpS2[i]);
909
910 fflush(NULL);
911 }
912 /*
913 * STEP 2) After rotation in STEP 1, in spherical coordinates, phi0
914 * and theta0 are given directly in Eq. (4.7), r0, pr0, ptheta0 and
915 * pphi0 are obtained by solving Eqs. (4.8) and (4.9) (using
916 * gsl_multiroot_fsolver). At this step, we find initial conditions
917 * for a spherical orbit without radiation reaction.
918 */
919
920 /* Initialise the gsl stuff */
921 XLAL_CALLGSL(rootSolver = gsl_multiroot_fsolver_alloc(T, 3));
922 if (!rootSolver) {
923 gsl_matrix_free(rotMatrix);
924 gsl_matrix_free(invMatrix);
926 }
927 XLAL_CALLGSL(initValues = gsl_vector_calloc(3));
928 if (!initValues) {
929 gsl_multiroot_fsolver_free(rootSolver);
930 gsl_matrix_free(rotMatrix);
931 gsl_matrix_free(invMatrix);
933 }
934
935 rootFunction.f = XLALFindSphericalOrbitPrec;
936 rootFunction.n = 3;
937 rootFunction.params = &rootParams;
938
939 /* Set to use optimized or unoptimized code */
940 rootParams.use_optimized = use_optimized;
941
942 /* Calculate the initial velocity from the given initial frequency */
943 omega = LAL_PI * mTotal * LAL_MTSUN_SI * fMin;
944 v0 = cbrt(omega);
945
946 /* Given this, we can start to calculate the initial conditions */
947 /* for spherical coords in the new basis */
948 rootParams.omega = omega;
949 rootParams.params = params;
950
951 /* To start with, we will just assign Newtonian-ish ICs to the system */
952
953
954 rootParams.values[0] = scale1 * sqrt(1. / (v0 * v0)*1/(v0*v0) - 36.0); /* Initial r */
955 rootParams.values[4] = scale2 * v0; /* Initial p */
956 rootParams.values[5] = scale3 * 1e-3;
957 //PK
958 memcpy(rootParams.values + 6, tmpS1, sizeof(tmpS1));
959 memcpy(rootParams.values + 9, tmpS2, sizeof(tmpS2));
960
961 if (printPK) {
962 XLAL_PRINT_INFO("ICs guess: x = %.16e, py = %.16e, pz = %.16e\n",
963 rootParams.values[0]/scale1, rootParams.values[4]/scale2,
964 rootParams.values[5]/scale3);
965 fflush(NULL);
966 }
967
968 gsl_vector_set(initValues, 0, rootParams.values[0]);
969 gsl_vector_set(initValues, 1, rootParams.values[4]);
970 gsl_vector_set(initValues, 2, rootParams.values[5]);
971
972 gsl_multiroot_fsolver_set(rootSolver, &rootFunction, initValues);
973
974 /* We are now ready to iterate to find the solution */
975 i = 0;
976 if(debugPK){
977 out = fopen("ICIterations.dat", "w");
978 }
979 INT4 jittered=0;
980 REAL8 r_now = 0.0;
981 REAL8 temp = 0.0;
982 do {
983 XLAL_CALLGSL(gslStatus = gsl_multiroot_fsolver_iterate(rootSolver));
984 if (debugPK) {
985 fprintf( out, "%d\t", i );
986
987 /* Write to file */
988 r_now = sqrt(rootParams.values[0]*rootParams.values[0]+36.0);
989 fprintf( out, "%.16e\t%.16e\t%.16e\t",
990 r_now, rootParams.values[4]/scale2,
991 rootParams.values[5]/scale3 );
992
993 /* Residual Function values whose roots we are trying to find */
994 finalValues = gsl_multiroot_fsolver_f(rootSolver);
995
996 /* Write to file */
997 fprintf( out, "%.16e\t%.16e\t%.16e\t",
998 gsl_vector_get(finalValues, 0),
999 gsl_vector_get(finalValues, 1),
1000 gsl_vector_get(finalValues, 2) );
1001
1002 /* Step sizes in each of function variables */
1003 finalValues = gsl_multiroot_fsolver_dx(rootSolver);
1004 temp = gsl_vector_get(finalValues, 0)/scale1/r_now;
1005 /* Write to file */
1006 fprintf( out, "%.16e\t%.16e\t%.16e\t%d\n",
1007 temp,
1008 gsl_vector_get(finalValues, 1)/scale2,
1009 gsl_vector_get(finalValues, 2)/scale3,
1010 gslStatus );
1011 }
1012
1013 if (gslStatus == GSL_ENOPROG || gslStatus == GSL_ENOPROGJ) {
1015 "\n NO PROGRESS being made by Spherical orbit root solver\n");
1016
1017 /* Print Residual Function values whose roots we are trying to find */
1018 finalValues = gsl_multiroot_fsolver_f(rootSolver);
1019 XLAL_PRINT_INFO("Function value here given by the following:\n");
1020 XLAL_PRINT_INFO(" F1 = %.16e, F2 = %.16e, F3 = %.16e\n",
1021 gsl_vector_get(finalValues, 0),
1022 gsl_vector_get(finalValues, 1), gsl_vector_get(finalValues, 2));
1023
1024 /* Print Step sizes in each of function variables */
1025 finalValues = gsl_multiroot_fsolver_dx(rootSolver);
1026// XLAL_PRINT_INFO("Stepsizes in each dimension:\n");
1027// XLAL_PRINT_INFO(" x = %.16e, py = %.16e, pz = %.16e\n",
1028// gsl_vector_get(finalValues, 0)/scale1,
1029// gsl_vector_get(finalValues, 1)/scale2,
1030// gsl_vector_get(finalValues, 2)/scale3);
1031
1032 /* Only allow this flag to be caught MAXcntGslNoProgress no. of times */
1033 cntGslNoProgress += 1;
1034 if (cntGslNoProgress >= MAXcntGslNoProgress) {
1035 cntGslNoProgress = 0;
1036
1037 if(multFacGslNoProgress < 1.){ multFacGslNoProgress *= 1.02; }
1038 else{ multFacGslNoProgress /= 1.01; }
1039
1040 }
1041 /* Now that no progress is being made, we need to reset the initial guess
1042 * for the (r,pPhi, pTheta) and reset the integrator */
1043 rootParams.values[0] = scale1 * sqrt(1. / (v0 * v0)*1./(v0*v0) -36.0); /* Initial r */
1044 rootParams.values[4] = scale2 * v0; /* Initial p */
1045 if( cntGslNoProgress % 2 )
1046 rootParams.values[5] = scale3 * 1e-3 / multFacGslNoProgress;
1047 else
1048 rootParams.values[5] = scale3 * 1e-3 * multFacGslNoProgress;
1049 //PK
1050 memcpy(rootParams.values + 6, tmpS1, sizeof(tmpS1));
1051 memcpy(rootParams.values + 9, tmpS2, sizeof(tmpS2));
1052
1053 if (printPK) {
1054 XLAL_PRINT_INFO("New ICs guess: x = %.16e, py = %.16e, pz = %.16e\n",
1055 rootParams.values[0]/scale1, rootParams.values[4]/scale2,
1056 rootParams.values[5]/scale3);
1057 fflush(NULL);
1058 }
1059
1060 gsl_vector_set(initValues, 0, rootParams.values[0]);
1061 gsl_vector_set(initValues, 1, rootParams.values[4]);
1062 gsl_vector_set(initValues, 2, rootParams.values[5]);
1063 gsl_multiroot_fsolver_set(rootSolver, &rootFunction, initValues);
1064 }
1065 else if (gslStatus == GSL_EBADFUNC) {
1067 "Inf or Nan encountered in evaluluation of spherical orbit Eqn\n");
1068 gsl_multiroot_fsolver_free(rootSolver);
1069 gsl_vector_free(initValues);
1070 gsl_matrix_free(rotMatrix);
1071 gsl_matrix_free(invMatrix);
1073 }
1074 else if (gslStatus != GSL_SUCCESS) {
1075 XLALPrintError("Error in GSL iteration function!\n");
1076 gsl_multiroot_fsolver_free(rootSolver);
1077 gsl_vector_free(initValues);
1078 gsl_matrix_free(rotMatrix);
1079 gsl_matrix_free(invMatrix);
1081 }
1082
1083 /* different ways to test convergence of the method */
1084 XLAL_CALLGSL(gslStatus = gsl_multiroot_test_residual(rootSolver->f, 1.0e-9));
1085 /*XLAL_CALLGSL(gslStatus= gsl_multiroot_test_delta(
1086 gsl_multiroot_fsolver_dx(rootSolver),
1087 gsl_multiroot_fsolver_root(rootSolver),
1088 1.e-8, 1.e-5));*/
1089
1090 if (jittered==0) {
1091 finalValues = gsl_multiroot_fsolver_dx(rootSolver);
1092 if (isnan(gsl_vector_get(finalValues, 1))) {
1093 rootParams.values[0] = scale1 * sqrt(1. / (v0 * v0)*1/(v0*v0)*(1.+1.e-8) - 36.0); /* Jitter on initial r */
1094 rootParams.values[4] = scale2 * v0*(1.-1.e-8); /* Jitter on initial p */
1095 rootParams.values[5] = scale3 * 1e-3;
1096 memcpy(rootParams.values + 6, tmpS1, sizeof(tmpS1));
1097 memcpy(rootParams.values + 9, tmpS2, sizeof(tmpS2));
1098 gsl_vector_set(initValues, 0, rootParams.values[0]);
1099 gsl_vector_set(initValues, 1, rootParams.values[4]);
1100 gsl_vector_set(initValues, 2, rootParams.values[5]);
1101 gsl_multiroot_fsolver_set(rootSolver, &rootFunction, initValues);
1102 jittered=1;
1103 }
1104 }
1105 i++;
1106 }
1107 while (gslStatus == GSL_CONTINUE && i <= maxIter);
1108
1109 if(debugPK) { fflush(NULL); fclose(out); }
1110
1111 if (i > maxIter && gslStatus != GSL_SUCCESS) {
1112 gsl_multiroot_fsolver_free(rootSolver);
1113 gsl_vector_free(initValues);
1114 gsl_matrix_free(rotMatrix);
1115 gsl_matrix_free(invMatrix);
1116 //XLAL_ERROR(XLAL_EMAXITER);
1118 }
1119 finalValues = gsl_multiroot_fsolver_root(rootSolver);
1120
1121 if (printPK) {
1122 XLAL_PRINT_INFO("Spherical orbit conditions here given by the following:\n");
1123 XLAL_PRINT_INFO(" x = %.16e, py = %.16e, pz = %.16e\n",
1124 gsl_vector_get(finalValues, 0)/scale1,
1125 gsl_vector_get(finalValues, 1)/scale2,
1126 gsl_vector_get(finalValues, 2)/scale3);
1127 }
1128 memset(qCart, 0, sizeof(qCart));
1129 memset(pCart, 0, sizeof(pCart));
1130
1131 qCart[0] = sqrt(gsl_vector_get(finalValues, 0)*gsl_vector_get(finalValues, 0)+36.0);
1132 pCart[1] = gsl_vector_get(finalValues, 1)/scale2;
1133 pCart[2] = gsl_vector_get(finalValues, 2)/scale3;
1134
1135
1136 /* Free the GSL root finder, since we're done with it */
1137 gsl_multiroot_fsolver_free(rootSolver);
1138 gsl_vector_free(initValues);
1139
1140
1141 /*
1142 * STEP 3) Rotate to L0 along z-axis and N0 along x-axis frame, where
1143 * L0 is the initial orbital angular momentum and L0 is calculated
1144 * using initial position and linear momentum obtained in STEP 2.
1145 */
1146
1147 /* Now we can calculate the relativistic L and rotate to a new basis */
1148 memcpy(qHat, qCart, sizeof(qCart));
1149 memcpy(pHat, pCart, sizeof(pCart));
1150
1151 NormalizeVectorPrec(qHat);
1152 NormalizeVectorPrec(pHat);
1153
1154 Lhat[0] = CalculateCrossProductPrec(0, qHat, pHat);
1155 Lhat[1] = CalculateCrossProductPrec(1, qHat, pHat);
1156 Lhat[2] = CalculateCrossProductPrec(2, qHat, pHat);
1157
1158 NormalizeVectorPrec(Lhat);
1159
1160 XLAL_CALLGSL(rotMatrix2 = gsl_matrix_alloc(3, 3));
1161 XLAL_CALLGSL(invMatrix2 = gsl_matrix_alloc(3, 3));
1162
1163 if (CalculateRotationMatrixPrec(rotMatrix2, invMatrix2, qHat, pHat, Lhat) == XLAL_FAILURE) {
1164 gsl_matrix_free(rotMatrix);
1165 gsl_matrix_free(invMatrix);
1167 }
1168 ApplyRotationMatrixPrec(rotMatrix2, rHat);
1169 ApplyRotationMatrixPrec(rotMatrix2, vHat);
1170 ApplyRotationMatrixPrec(rotMatrix2, LnHat);
1171 ApplyRotationMatrixPrec(rotMatrix2, tmpS1);
1172 ApplyRotationMatrixPrec(rotMatrix2, tmpS2);
1173 ApplyRotationMatrixPrec(rotMatrix2, tmpS1Norm);
1174 ApplyRotationMatrixPrec(rotMatrix2, tmpS2Norm);
1175 ApplyRotationMatrixPrec(rotMatrix2, qCart);
1176 ApplyRotationMatrixPrec(rotMatrix2, pCart);
1177
1178 gsl_matrix_free(rotMatrix);
1179 gsl_matrix_free(rotMatrix2);
1180
1181 if (printPK) {
1182 XLAL_PRINT_INFO("qCart after rotation2 %3.10f %3.10f %3.10f\n", qCart[0], qCart[1], qCart[2]);
1183 XLAL_PRINT_INFO("pCart after rotation2 %3.10f %3.10f %3.10f\n", pCart[0], pCart[1], pCart[2]);
1184 XLAL_PRINT_INFO("S1 after rotation2 %3.10f %3.10f %3.10f\n", tmpS1Norm[0], tmpS1Norm[1], tmpS1Norm[2]);
1185 XLAL_PRINT_INFO("S2 after rotation2 %3.10f %3.10f %3.10f\n", tmpS2Norm[0], tmpS2Norm[1], tmpS2Norm[2]);
1186 }
1187 /*
1188 * STEP 4) In the L0-N0 frame, we calculate (dE/dr)|sph using Eq.
1189 * (4.14), then initial dr/dt using Eq. (4.10), and finally pr0 using
1190 * Eq. (4.15).
1191 */
1192
1193 /* Now we can calculate the flux. Change to spherical co-ords */
1194 CartesianToSphericalPrec(qSph, pSph, qCart, pCart);
1195 memcpy(sphValues, qSph, sizeof(qSph));
1196 memcpy(sphValues + 3, pSph, sizeof(pSph));
1197 memcpy(sphValues + 6, tmpS1, sizeof(tmpS1));
1198 memcpy(sphValues + 9, tmpS2, sizeof(tmpS2));
1199
1200 memcpy(cartValues, qCart, sizeof(qCart));
1201 memcpy(cartValues + 3, pCart, sizeof(pCart));
1202 memcpy(cartValues + 6, tmpS1, sizeof(tmpS1));
1203 memcpy(cartValues + 9, tmpS2, sizeof(tmpS2));
1204
1205 REAL8 dHdpphi , d2Hdr2, d2Hdrdpphi;
1206 REAL8 rDot , dHdpr, flux, dEdr;
1207
1208 d2Hdr2 = XLALCalculateSphHamiltonianDeriv2Prec(0, 0, sphValues, params, use_optimized);
1209 d2Hdrdpphi = XLALCalculateSphHamiltonianDeriv2Prec(0, 5, sphValues, params, use_optimized);
1210
1211 if (printPK)
1212 XLAL_PRINT_INFO("d2Hdr2 = %.16e, d2Hdrdpphi = %.16e\n", d2Hdr2, d2Hdrdpphi);
1213
1214 /* New code to compute derivatives w.r.t. cartesian variables */
1215
1216 REAL8 tmpDValues[14];
1217 int UNUSED status;
1218 for (i = 0; i < 3; i++) {
1219 cartValues[i + 6] /= mTotal * mTotal;
1220 cartValues[i + 9] /= mTotal * mTotal;
1221 }
1222 UINT4 oldignoreflux = params->ignoreflux;
1223 params->ignoreflux = 1;
1224 status = XLALSpinPrecHcapNumericalDerivative(0, cartValues, tmpDValues, params);
1225 params->ignoreflux = oldignoreflux;
1226 for (i = 0; i < 3; i++) {
1227 cartValues[i + 6] *= mTotal * mTotal;
1228 cartValues[i + 9] *= mTotal * mTotal;
1229 }
1230
1231 dHdpphi = tmpDValues[1] / sqrt(cartValues[0] * cartValues[0] + cartValues[1] * cartValues[1] + cartValues[2] * cartValues[2]);
1232 //XLALSpinPrecHcapNumDerivWRTParam(4, cartValues, params) / sphValues[0];
1233
1234 dEdr = -dHdpphi * d2Hdr2 / d2Hdrdpphi;
1235
1236 if (printPK)
1237 XLAL_PRINT_INFO("d2Hdr2 = %.16e d2Hdrdpphi = %.16e dHdpphi = %.16e\n",
1238 d2Hdr2, d2Hdrdpphi, dHdpphi);
1239
1240 if (d2Hdr2 != 0.0) {
1241 /* We will need to calculate the Hamiltonian to get the flux */
1242 s1Vec.length = s2Vec.length = s1VecNorm.length = s2VecNorm.length = sKerr.length = sStar.length = 3;
1243 s1Vec.data = tmpS1;
1244 s2Vec.data = tmpS2;
1245 s1VecNorm.data = tmpS1Norm;
1246 s2VecNorm.data = tmpS2Norm;
1247 sKerr.data = sKerrData;
1248 sStar.data = sStarData;
1249
1250 qCartVec.length = pCartVec.length = 3;
1251 qCartVec.data = qCart;
1252 pCartVec.data = pCart;
1253
1254 //chi1 = tmpS1[0] * LnHat[0] + tmpS1[1] * LnHat[1] + tmpS1[2] * LnHat[2];
1255 //chi2 = tmpS2[0] * LnHat[0] + tmpS2[1] * LnHat[1] + tmpS2[2] * LnHat[2];
1256
1257 //if (debugPK)
1258 //XLAL_PRINT_INFO("magS1 = %.16e, magS2 = %.16e\n", chi1, chi2);
1259
1260 //chiS = 0.5 * (chi1 / (mass1 * mass1) + chi2 / (mass2 * mass2));
1261 //chiA = 0.5 * (chi1 / (mass1 * mass1) - chi2 / (mass2 * mass2));
1262
1263 XLALSimIMRSpinEOBCalculateSigmaKerr(&sKerr, mass1, mass2, &s1Vec, &s2Vec);
1264 XLALSimIMRSpinEOBCalculateSigmaStar(&sStar, mass1, mass2, &s1Vec, &s2Vec);
1265
1266 /*
1267 * The a in the flux has been set to zero, but not in the
1268 * Hamiltonian
1269 */
1270 a = sqrt(sKerr.data[0] * sKerr.data[0] + sKerr.data[1] * sKerr.data[1] + sKerr.data[2] * sKerr.data[2]);
1271 //XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients(params->eobParams->hCoeffs, mass1, mass2, eta, /* a */ 0.0, chiS, chiA);
1272 //XLALSimIMRCalculateSpinPrecEOBHCoeffs(params->seobCoeffs, eta, a);
1273 ham = XLALSimIMRSpinPrecEOBHamiltonian(eta, &qCartVec, &pCartVec, &s1VecNorm, &s2VecNorm, &sKerr, &sStar, params->tortoise, params->seobCoeffs);
1274
1275 if (printPK)
1276 XLAL_PRINT_INFO("Stas: hamiltonian in ICs at this point is %.16e\n", ham);
1277
1278 /* And now, finally, the flux */
1279 REAL8Vector polarDynamics, cartDynamics;
1280 REAL8 polarData[4], cartData[12];
1281
1282 polarDynamics.length = 4;
1283 polarDynamics.data = polarData;
1284
1285 polarData[0] = qSph[0];
1286 polarData[1] = 0.;
1287 polarData[2] = pSph[0];
1288 polarData[3] = pSph[2];
1289
1290 cartDynamics.length = 12;
1291 cartDynamics.data = cartData;
1292
1293 memcpy(cartData, qCart, 3 * sizeof(REAL8));
1294 memcpy(cartData + 3, pCart, 3 * sizeof(REAL8));
1295 memcpy(cartData + 6, tmpS1Norm, 3 * sizeof(REAL8));
1296 memcpy(cartData + 9, tmpS2Norm, 3 * sizeof(REAL8));
1297
1298 //XLAL_PRINT_INFO("Stas: starting FLux calculations\n");
1299 if(use_optimized) {
1300 flux = XLALInspiralPrecSpinFactorizedFlux_exact(&polarDynamics, &cartDynamics, nqcCoeffs, omega, params, ham, lMax, SpinAlignedEOBversion);
1301 } else {
1302 flux = XLALInspiralPrecSpinFactorizedFlux(&polarDynamics, &cartDynamics, nqcCoeffs, omega, params, ham, lMax, SpinAlignedEOBversion);
1303 }
1304 /*
1305 * flux = XLALInspiralSpinFactorizedFlux( &polarDynamics,
1306 * nqcCoeffs, omega, params, ham, lMax, SpinAlignedEOBversion
1307 * );
1308 */
1309 //XLAL_PRINT_INFO("Stas flux = %.16e \n", flux);
1310 //exit(0);
1311 flux = flux / eta;
1312
1313 rDot = -flux / dEdr;
1314 if (debugPK) {
1315 XLAL_PRINT_INFO("Stas here I am 2 \n");
1316 }
1317 /*
1318 * We now need dHdpr - we take it that it is safely linear up
1319 * to a pr of 1.0e-3 PK: Ideally, the pr should be of the
1320 * order of other momenta, in order for its contribution to
1321 * the Hamiltonian to not get buried in the numerical noise
1322 * in the numerically larger momenta components
1323 */
1324 cartValues[3] = 1.0e-3;
1325 for (i = 0; i < 3; i++) {
1326 cartValues[i + 6] /= mTotal * mTotal;
1327 cartValues[i + 9] /= mTotal * mTotal;
1328 }
1329 oldignoreflux = params->ignoreflux;
1330 params->ignoreflux = 1;
1331 params->seobCoeffs->updateHCoeffs = 1;
1332 status = XLALSpinPrecHcapNumericalDerivative(0, cartValues, tmpDValues, params);
1333 params->ignoreflux = oldignoreflux;
1334 for (i = 0; i < 3; i++) {
1335 cartValues[i + 6] *= mTotal * mTotal;
1336 cartValues[i + 9] *= mTotal * mTotal;
1337 }
1338 REAL8 csi = sqrt(XLALSimIMRSpinPrecEOBHamiltonianDeltaT(params->seobCoeffs, qSph[0], eta, a)*XLALSimIMRSpinPrecEOBHamiltonianDeltaR(params->seobCoeffs, qSph[0], eta, a)) / (qSph[0] * qSph[0] + a * a);
1339
1340 dHdpr = csi*tmpDValues[0];
1341 //XLALSpinPrecHcapNumDerivWRTParam(3, cartValues, params);
1342
1343 if (debugPK) {
1344 XLAL_PRINT_INFO("Ingredients going into prDot:\n");
1345 XLAL_PRINT_INFO("flux = %.16e, dEdr = %.16e, dHdpr = %.16e, dHdpr/pr = %.16e\n", flux, dEdr, dHdpr, dHdpr / cartValues[3]);
1346 }
1347 /*
1348 * We can now calculate what pr should be taking into account
1349 * the flux
1350 */
1351 pSph[0] = rDot / (dHdpr / cartValues[3]);
1352 } else {
1353 /*
1354 * Since d2Hdr2 has evaluated to zero, we cannot do the
1355 * above. Just set pr to zero
1356 */
1357 //XLAL_PRINT_INFO("d2Hdr2 is zero!\n");
1358 pSph[0] = 0;
1359 }
1360
1361 /* Now we are done - convert back to cartesian coordinates ) */
1362 SphericalToCartesianPrec(qCart, pCart, qSph, pSph);
1363
1364 /*
1365 * STEP 5) Rotate back to the original inertial frame by inverting
1366 * the rotation of STEP 3 and then inverting the rotation of STEP 1.
1367 */
1368
1369 /* Undo rotations to get back to the original basis */
1370 /* Second rotation */
1371 ApplyRotationMatrixPrec(invMatrix2, rHat);
1372 ApplyRotationMatrixPrec(invMatrix2, vHat);
1373 ApplyRotationMatrixPrec(invMatrix2, LnHat);
1374 ApplyRotationMatrixPrec(invMatrix2, tmpS1);
1375 ApplyRotationMatrixPrec(invMatrix2, tmpS2);
1376 ApplyRotationMatrixPrec(invMatrix2, tmpS1Norm);
1377 ApplyRotationMatrixPrec(invMatrix2, tmpS2Norm);
1378 ApplyRotationMatrixPrec(invMatrix2, qCart);
1379 ApplyRotationMatrixPrec(invMatrix2, pCart);
1380
1381 /* First rotation */
1382 ApplyRotationMatrixPrec(invMatrix, rHat);
1383 ApplyRotationMatrixPrec(invMatrix, vHat);
1384 ApplyRotationMatrixPrec(invMatrix, LnHat);
1385 ApplyRotationMatrixPrec(invMatrix, tmpS1);
1386 ApplyRotationMatrixPrec(invMatrix, tmpS2);
1387 ApplyRotationMatrixPrec(invMatrix, tmpS1Norm);
1388 ApplyRotationMatrixPrec(invMatrix, tmpS2Norm);
1389 ApplyRotationMatrixPrec(invMatrix, qCart);
1390 ApplyRotationMatrixPrec(invMatrix, pCart);
1391
1392 gsl_matrix_free(invMatrix);
1393 gsl_matrix_free(invMatrix2);
1394
1395 /* If required, apply the tortoise transform */
1396 if (tmpTortoise) {
1397 REAL8 r = sqrt(qCart[0] * qCart[0] + qCart[1] * qCart[1] + qCart[2] * qCart[2]);
1400 REAL8 csi = sqrt(deltaT * deltaR) / (r * r + a * a);
1401
1402 REAL8 pr = (qCart[0] * pCart[0] + qCart[1] * pCart[1] + qCart[2] * pCart[2]) / r;
1403
1404 params->tortoise = tmpTortoise;
1405
1406 if (debugPK) {
1407 XLAL_PRINT_INFO("Applying the tortoise to p (csi = %.26e)\n", csi);
1408 XLAL_PRINT_INFO("pCart = %3.10f %3.10f %3.10f\n", pCart[0], pCart[1], pCart[2]);
1409 }
1410 for (i = 0; i < 3; i++) {
1411 pCart[i] = pCart[i] + qCart[i] * pr * (csi - 1.) / r;
1412 }
1413 }
1414
1415
1416 /* Now copy the initial conditions back to the return vector */
1417 memcpy(initConds->data, qCart, sizeof(qCart));
1418 memcpy(initConds->data + 3, pCart, sizeof(pCart));
1419 memcpy(initConds->data + 6, tmpS1Norm, sizeof(tmpS1Norm));
1420 memcpy(initConds->data + 9, tmpS2Norm, sizeof(tmpS2Norm));
1421
1422 for (i=0; i<12; i++) {
1423 if (fabs(initConds->data[i]) <=1.0e-15) {
1424 initConds->data[i] = 0.;
1425 }
1426 }
1427
1428 if (debugPK) {
1429 XLAL_PRINT_INFO("THE FINAL INITIAL CONDITIONS:\n");
1430 XLAL_PRINT_INFO(" %.16e %.16e %.16e\n%.16e %.16e %.16e\n%.16e %.16e %.16e\n%.16e %.16e %.16e\n", initConds->data[0], initConds->data[1], initConds->data[2],
1431 initConds->data[3], initConds->data[4], initConds->data[5], initConds->data[6], initConds->data[7], initConds->data[8],
1432 initConds->data[9], initConds->data[10], initConds->data[11]);
1433 }
1434 return XLAL_SUCCESS;
1435}
1436
1437#endif /* _LALSIMIMRSPINPRECEOBINITIALCONDITIONS_C */
static int XLALSimIMRSpinEOBCalculateSigmaKerr(REAL8Vector *sigmaKerr, REAL8 mass1, REAL8 mass2, REAL8Vector *s1, REAL8Vector *s2)
Function to calculate normalized spin of the deformed-Kerr background in SEOBNRv1.
static int XLALSimIMRSpinEOBCalculateSigmaStar(REAL8Vector *sigmaStar, REAL8 mass1, REAL8 mass2, REAL8Vector *s1, REAL8Vector *s2)
Function to calculate normalized spin of the test particle in SEOBNRv1.
static REAL8 XLALInspiralPrecSpinFactorizedFlux(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 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 REAL8 XLALSimIMRSpinPrecEOBHamiltonian(const REAL8 eta, REAL8Vector *restrict x, REAL8Vector *restrict p, REAL8Vector *restrict s1Vec, REAL8Vector *restrict s2Vec, REAL8Vector *restrict sigmaKerr, REAL8Vector *restrict sigmaStar, int tortoise, SpinEOBHCoeffs *coeffs)
static REAL8 XLALSimIMRSpinPrecEOBHamiltonianDeltaT(SpinEOBHCoeffs *coeffs, const REAL8 r, const REAL8 eta, const REAL8 a)
This function calculates the function which appears in the spinning EOB potential function.
static UNUSED REAL8 XLALSimIMRSpinPrecEOBHamiltonianDeltaR(SpinEOBHCoeffs *coeffs, const REAL8 r, const REAL8 eta, const REAL8 a)
This function calculates the DeltaR potential function in the spin EOB Hamiltonian.
static REAL8 XLALSpinPrecHcapNumDerivWRTParam(const INT4 paramIdx, const REAL8 values[], SpinEOBParams *funcParams)
Calculate the derivative of the Hamiltonian w.r.t.
static int XLALSpinPrecHcapNumericalDerivative(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 int XLALFindSphericalOrbitPrec(const gsl_vector *x, void *params, gsl_vector *f)
Root function for gsl root finders.
static double GSLSpinHamiltonianDerivWrapperPrec(double x, void *params)
Wrapper for calculating specific derivative of the Hamiltonian in spherical co-ordinates,...
static int NormalizeVectorPrec(REAL8 a[])
Normalizes the given vector.
static REAL8 CalculateCrossProductPrec(const INT4 i, const REAL8 a[], const REAL8 b[])
Calculates the ith component of the cross product of two vectors, where i is in the range 0-2.
static int XLALSimIMRSpinEOBInitialConditionsPrec(REAL8Vector *initConds, const REAL8 mass1, const REAL8 mass2, const REAL8 fMin, const REAL8 inc, const REAL8 spin1[], const REAL8 spin2[], SpinEOBParams *params, INT4 use_optimized)
Main function for calculating the spinning EOB initial conditions, following the quasi-spherical init...
static int CalculateRotationMatrixPrec(gsl_matrix *rotMatrix, gsl_matrix *rotInverse, REAL8 r[], REAL8 v[], REAL8 L[])
Calculate the rotation matrix and its inverse.
static int ApplyRotationMatrixPrec(gsl_matrix *rotMatrix, REAL8 a[])
Applies a rotation matrix to a given vector.
static int SphericalToCartesianPrec(REAL8 qCart[], REAL8 pCart[], const REAL8 qSph[], const REAL8 pSph[])
Performs a co-ordinate transformation from spherical to Cartesian co-ordinates.
static REAL8 XLALCalculateSphHamiltonianDeriv2Prec(const int idx1, const int idx2, const REAL8 values[], SpinEOBParams *params, INT4 use_optimized)
Function to calculate the second derivative of the Hamiltonian.
static int XLALRobustDerivative(const gsl_function *F, double x, double h, double *result, double *absErr)
static REAL8 CalculateDotProductPrec(const REAL8 a[], const REAL8 b[])
Calculate the dot product of two vectors.
static int CartesianToSphericalPrec(REAL8 qSph[], REAL8 pSph[], const REAL8 qCart[], const REAL8 pCart[])
Perform a co-ordinate transformation from Cartesian to spherical co-ordinates.
#define fprintf
#define XLAL_CALLGSL(statement)
double i
Definition: bh_ringdown.c:118
double e
Definition: bh_ringdown.c:117
const double pr
const double deltaR
const double csi
#define LAL_PI_2
#define LAL_PI
#define LAL_MTSUN_SI
double REAL8
uint32_t UINT4
int32_t INT4
static const INT4 r
static const INT4 a
#define XLAL_ERROR_REAL8(...)
#define XLAL_PRINT_INFO(...)
#define XLAL_ERROR(...)
#define XLAL_PRINT_WARNING(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_PRINT_ERROR(...)
#define XLAL_IS_REAL8_FAIL_NAN(val)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EDOM
XLAL_EINVAL
XLAL_FAILURE
list y
string status
x
The coefficients which are used in calculating the non-quasicircular correction to the EOBNRv2 model.
const REAL8 * sphValues
SpinEOBParams * params
REAL8 * data
Structure consisting SEOBNR parameters that can be used by gsl root finders.
REAL8 omega
< Orbital frequency
REAL8 values[12]
< Dynamical variables, x, y, z, px, py, pz, S1x, S1y, S1z, S2x, S2y and S2z
SpinEOBParams * params
< Spin EOB parameters – physical, pre-computed, etc.
Parameters for the spinning EOB model.
EOBParams * eobParams
Definition: burst.c:245
double deltaT
Definition: unicorn.c:24