1 #include <lal/LALSimInspiral.h>
2 #include <lal/LALSimIMR.h>
4 #include <gsl/gsl_vector.h>
5 #include <gsl/gsl_matrix.h>
6 #include <gsl/gsl_blas.h>
8 #include <gsl/gsl_multiroots.h>
9 #include <gsl/gsl_deriv.h>
20 #ifndef _LALSIMIMRSPINEOBINITIALCONDITIONS_C
21 #define _LALSIMIMRSPINEOBINITIALCONDITIONS_C
27 struct tagSEOBRootParams
43 return a[0]*b[0] +
a[1]*b[1] +
a[2]*b[2];
55 return a[(
i+1)%3]*b[(
i+2)%3] -
a[(
i+2)%3]*b[(
i+1)%3];
66 REAL8 norm = sqrt(
a[0]*
a[0] +
a[1]*
a[1] +
a[2]*
a[2] );
82 gsl_matrix *rotMatrix,
83 gsl_matrix *rotInverse,
93 REAL8 cosa, sina, cosb, sinb, cosg, sing;
104 a = atan2( L[0], -L[1] );
105 b = atan2( sqrt( L[0]*L[0] + L[1]*L[1] ), L[2] );
106 g = atan2(
r[2], v[2] );
109 if ( ( cosa = cos(
a ) ) < 1.0e-16 ) cosa = 0.0;
110 if ( ( sina = sin(
a ) ) < 1.0e-16 ) sina = 0.0;
111 if ( ( cosb = cos( b ) ) < 1.0e-16 ) cosb = 0.0;
112 if ( ( sinb = sin( b ) ) < 1.0e-16 ) sinb = 0.0;
113 if ( ( cosg = cos( g ) ) < 1.0e-16 ) cosg = 0.0;
114 if ( ( sing = sin( g ) ) < 1.0e-16 ) sing = 0.0;
123 gsl_matrix_set( rotMatrix, 0, 0, cosg*cosa - cosb*sina*sing );
124 gsl_matrix_set( rotMatrix, 0, 1, cosg*sina + cosb*cosa*sing );
125 gsl_matrix_set( rotMatrix, 0, 2, sing*sinb );
126 gsl_matrix_set( rotMatrix, 1, 0, -sing*cosa - cosb*sina*cosg );
127 gsl_matrix_set( rotMatrix, 1, 1, -sing*sina + cosb*cosa*cosg );
128 gsl_matrix_set( rotMatrix, 1, 2, cosg*sinb );
129 gsl_matrix_set( rotMatrix, 2, 0, sinb*sina );
130 gsl_matrix_set( rotMatrix, 2, 1, -sinb*cosa );
131 gsl_matrix_set( rotMatrix, 2, 2, cosb );
134 gsl_matrix_transpose_memcpy( rotInverse, rotMatrix );
171 gsl_matrix *rotMatrix,
176 gsl_vector *tmpVec1 = gsl_vector_alloc( 3 );
177 gsl_vector *tmpVec2 = gsl_vector_calloc( 3 );
179 gsl_vector_set( tmpVec1, 0,
a[0] );
180 gsl_vector_set( tmpVec1, 1,
a[1] );
181 gsl_vector_set( tmpVec1, 2,
a[2] );
183 gsl_blas_dgemv( CblasNoTrans, 1.0, rotMatrix, tmpVec1, 0.0, tmpVec2 );
185 a[0] = gsl_vector_get( tmpVec2, 0 );
186 a[1] = gsl_vector_get( tmpVec2, 1 );
187 a[2] = gsl_vector_get( tmpVec2, 2 );
189 gsl_vector_free( tmpVec1 );
190 gsl_vector_free( tmpVec2 );
298 REAL8 py, pz,
r, ptheta, pphi;
301 REAL8 dHdx, dHdpy, dHdpz;
302 REAL8 dHdr, dHdptheta, dHdpphi;
306 rootParams->
values[0] =
r = gsl_vector_get(
x, 0 );
307 rootParams->
values[4] = py = gsl_vector_get(
x, 1 );
308 rootParams->
values[5] = pz = gsl_vector_get(
x, 2 );
340 dHdr = dHdx - dHdpy * pphi / (
r*
r) + dHdpz * ptheta / (
r*
r);
341 dHdptheta = - dHdpz /
r;
345 gsl_vector_set( f, 0, dHdr );
346 gsl_vector_set( f, 1, dHdptheta );
347 gsl_vector_set( f, 2, dHdpphi - rootParams->
omega );
366 REAL8 cartValues[12];
368 REAL8 dHdr, dHdx, dHdpy, dHdpz;
371 memcpy( sphValues, dParams->
sphValues,
sizeof( sphValues ) );
375 memcpy( cartValues+6, sphValues+6, 6*
sizeof(
REAL8) );
378 ptheta = sphValues[4];
391 dHdr = dHdx - dHdpy * pphi / (
r*
r) + dHdpz * ptheta / (
r*
r);
407 XLALPrintError(
"This option is not supported in the second derivative function!\n" );
426 REAL8 cartValues[12];
428 REAL8 dHdr, dHdx, dHdpy, dHdpz;
431 memcpy( sphValues, dParams->
sphValues,
sizeof( sphValues ) );
435 memcpy( cartValues+6, sphValues+6, 6*
sizeof(
REAL8) );
438 ptheta = sphValues[4];
450 dHdr = dHdx - dHdpy * pphi / (
r*
r) + dHdpz * ptheta / (
r*
r);
466 XLALPrintError(
"This option is not supported in the second derivative function!\n" );
474 const REAL8 values[],
478 static const REAL8 STEP_SIZE = 1.0e-4;
486 INT4 UNUSED gslStatus;
495 XLAL_CALLGSL( gslStatus = gsl_deriv_central( &F, values[idx1],
496 STEP_SIZE, &result, &absErr ) );
498 if ( gslStatus != GSL_SUCCESS )
500 XLALPrintError(
"XLAL Error %s - Failure in GSL function\n", __func__ );
512 const REAL8 values[],
517 static const REAL8 STEP_SIZE = 1.0e-4;
525 INT4 UNUSED gslStatus;
550 XLAL_CALLGSL( gslStatus = gsl_deriv_central( &F, values[idx1],
551 STEP_SIZE, &result, &absErr ) );
553 if ( gslStatus != GSL_SUCCESS )
555 XLALPrintError(
"XLAL Error %s - Failure in GSL function\n", __func__ );
604 INT4 use_optimized_v2
614 static const int UNUSED lMax = 8;
621 UINT4 SpinAlignedEOBversion;
637 REAL8 qCart[3], pCart[3];
638 REAL8 qSph[3], pSph[3];
650 REAL8 sKerrData[3], sStarData[3];
656 REAL8 cartValues[12];
660 gsl_matrix *rotMatrix = NULL;
661 gsl_matrix *invMatrix = NULL;
662 gsl_matrix *rotMatrix2 = NULL;
663 gsl_matrix *invMatrix2 = NULL;
667 const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids;
668 gsl_multiroot_fsolver *rootSolver = NULL;
670 gsl_multiroot_function rootFunction;
671 gsl_vector *initValues = NULL;
672 gsl_vector *finalValues = NULL;
674 const int maxIter = 100;
676 memset( &rootParams, 0,
sizeof( rootParams ) );
678 mTotal = mass1 + mass2;
679 eta = mass1 * mass2 / (mTotal * mTotal);
680 memcpy( tmpS1, spin1,
sizeof(tmpS1) );
681 memcpy( tmpS2, spin2,
sizeof(tmpS2) );
682 memcpy( tmpS1Norm, spin1,
sizeof(tmpS1Norm) );
683 memcpy( tmpS2Norm, spin2,
sizeof(tmpS2Norm) );
684 for (
i = 0;
i < 3;
i++ )
686 tmpS1Norm[
i] /= mTotal * mTotal;
687 tmpS2Norm[
i] /= mTotal * mTotal;
689 SpinAlignedEOBversion =
params->seobCoeffs->SpinAlignedEOBversion;
691 tmpTortoise =
params->tortoise;
695 nqcCoeffs =
params->nqcCoeffs;
707 if ( LnHat[2] > 0.9999 )
710 rHat[1] = rHat[2] = 0.;
714 REAL8 theta0 = atan( - LnHat[2] / LnHat[0] );
715 rHat[0] = sin( theta0 );
717 rHat[2] = cos( theta0 );
742 if ( !rotMatrix || !invMatrix )
744 if ( rotMatrix ) gsl_matrix_free( rotMatrix );
745 if ( invMatrix ) gsl_matrix_free( invMatrix );
751 gsl_matrix_free( rotMatrix );
752 gsl_matrix_free( invMatrix );
789 rootParams.
omega = omega;
793 rootParams.
values[0] = 1./(v0*v0);
794 rootParams.
values[4] = v0;
795 memcpy( rootParams.
values+6, tmpS1,
sizeof( tmpS1 ) );
796 memcpy( rootParams.
values+9, tmpS2,
sizeof( tmpS2 ) );
801 XLAL_CALLGSL( rootSolver = gsl_multiroot_fsolver_alloc( T, 3 ) );
804 gsl_matrix_free( rotMatrix );
805 gsl_matrix_free( invMatrix );
812 gsl_multiroot_fsolver_free( rootSolver );
813 gsl_matrix_free( rotMatrix );
814 gsl_matrix_free( invMatrix );
818 gsl_vector_set( initValues, 0, rootParams.
values[0] );
819 gsl_vector_set( initValues, 1, rootParams.
values[4] );
823 rootFunction.params = &rootParams;
825 gsl_multiroot_fsolver_set( rootSolver, &rootFunction, initValues );
832 XLAL_CALLGSL( gslStatus = gsl_multiroot_fsolver_iterate( rootSolver ) );
833 if ( gslStatus != GSL_SUCCESS )
836 gsl_multiroot_fsolver_free( rootSolver );
837 gsl_vector_free( initValues );
838 gsl_matrix_free( rotMatrix );
839 gsl_matrix_free( invMatrix );
843 XLAL_CALLGSL( gslStatus = gsl_multiroot_test_residual( rootSolver->f, 1.0e-10 ) );
846 while ( gslStatus == GSL_CONTINUE &&
i <= maxIter );
848 if (
i > maxIter && gslStatus != GSL_SUCCESS )
850 gsl_multiroot_fsolver_free( rootSolver );
851 gsl_vector_free( initValues );
852 gsl_matrix_free( rotMatrix );
853 gsl_matrix_free( invMatrix );
857 finalValues = gsl_multiroot_fsolver_root( rootSolver );
863 memset( qCart, 0,
sizeof(qCart) );
864 memset( pCart, 0,
sizeof(pCart) );
866 qCart[0] = gsl_vector_get( finalValues, 0 );
867 pCart[1] = gsl_vector_get( finalValues, 1 );
868 pCart[2] = gsl_vector_get( finalValues, 2 );
871 gsl_multiroot_fsolver_free( rootSolver );
872 gsl_vector_free( initValues );
879 memcpy( qHat, qCart,
sizeof(qCart) );
880 memcpy( pHat, pCart,
sizeof(pCart) );
896 gsl_matrix_free( rotMatrix );
897 gsl_matrix_free( invMatrix );
917 memcpy( sphValues, qSph,
sizeof( qSph ) );
918 memcpy( sphValues+3, pSph,
sizeof( pSph ) );
919 memcpy( sphValues+6, tmpS1,
sizeof(tmpS1) );
920 memcpy( sphValues+9, tmpS2,
sizeof(tmpS2) );
922 memcpy( cartValues, qCart,
sizeof(qCart) );
923 memcpy( cartValues+3, pCart,
sizeof(pCart) );
924 memcpy( cartValues+6, tmpS1,
sizeof(tmpS1) );
925 memcpy( cartValues+9, tmpS2,
sizeof(tmpS2) );
927 REAL8 dHdpphi, d2Hdr2, d2Hdrdpphi;
928 REAL8 rDot, dHdpr, flux, dEdr;
930 if(use_optimized_v2) {
942 dEdr = - dHdpphi * d2Hdr2 / d2Hdrdpphi;
952 s1VecNorm.
data = tmpS1Norm;
953 s2VecNorm.
data = tmpS2Norm;
954 sKerr.
data = sKerrData;
955 sStar.
data = sStarData;
958 qCartVec.
data = qCart;
959 pCartVec.
data = pCart;
976 if(use_optimized_v2) {
988 polarDynamics.
data = polarData;
990 polarData[0] = qSph[0];
992 polarData[2] = pSph[0];
993 polarData[3] = pSph[2];
998 rDot = - flux / dEdr;
1001 cartValues[3] = 1.0e-3;
1002 if(use_optimized_v2) {
1011 pSph[0] = rDot / (dHdpr / cartValues[3] );
1053 REAL8 r = sqrt(qCart[0]*qCart[0] + qCart[1]*qCart[1] + qCart[2]*qCart[2] );
1058 REAL8 pr = (qCart[0]*pCart[0] + qCart[1]*pCart[1] + qCart[2]*pCart[2])/
r;
1060 params->tortoise = tmpTortoise;
1064 for (
i = 0;
i < 3;
i++ )
1066 pCart[
i] = pCart[
i] + qCart[
i] *
pr * (
csi - 1.) /
r;
1071 memcpy( initConds->
data, qCart,
sizeof(qCart) );
1072 memcpy( initConds->
data+3, pCart,
sizeof(pCart) );
1073 memcpy( initConds->
data+6, tmpS1Norm,
sizeof(tmpS1Norm) );
1074 memcpy( initConds->
data+9, tmpS2Norm,
sizeof(tmpS2Norm) );
1076 gsl_matrix_free(rotMatrix2);
1077 gsl_matrix_free(invMatrix2);
1079 gsl_matrix_free(rotMatrix);
1080 gsl_matrix_free(invMatrix);
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 XLALInspiralSpinFactorizedFlux(REAL8Vector *values, EOBNonQCCoeffs *nqcCoeffs, const REAL8 omega, SpinEOBParams *ak, const REAL8 H, const UINT4 lMax, const UINT4 SpinAlignedEOBversion)
static REAL8 XLALSimIMRSpinEOBHamiltonianDeltaR(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 XLALSimIMRSpinEOBHamiltonian(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 XLALSimIMRSpinEOBHamiltonianDeltaT(SpinEOBHCoeffs *coeffs, const REAL8 r, const REAL8 eta, const REAL8 a)
This function calculates the function which appears in the spinning EOB potential function.
UNUSED REAL8 XLALSimIMRSpinEOBHamiltonianOptimized(const REAL8 eta, REAL8Vector *restrict x, REAL8Vector *restrict p, REAL8Vector *restrict s1Vec, REAL8Vector *restrict s2Vec, REAL8Vector *restrict sigmaKerr, REAL8Vector *restrict sigmaStar, INT4 tortoise, SpinEOBHCoeffs *coeffs)
Function to calculate the value of the spinning Hamiltonian for given values of the dynamical variabl...
static REAL8 XLALSpinHcapExactDerivWRTParam(const INT4 paramIdx, const REAL8 values[], SpinEOBParams *params)
Calculate the derivative of the Hamiltonian w.r.t.
static REAL8 XLALSpinHcapHybDerivWRTParam(const INT4 paramIdx, const REAL8 values[], SpinEOBParams *funcParams)
static REAL8 XLALSpinHcapNumDerivWRTParam(const INT4 paramIdx, const REAL8 values[], SpinEOBParams *params)
Calculate the derivative of the Hamiltonian w.r.t.
static int ApplyRotationMatrix(gsl_matrix *rotMatrix, REAL8 a[])
Applies a rotation matrix to a given vector.
static REAL8 XLALCalculateSphHamiltonianDeriv2(const int idx1, const int idx2, const REAL8 values[], SpinEOBParams *params)
static double GSLSpinHamiltonianDerivWrapperHybrid(double x, void *params)
Wrapper for calculating specific derivative of the Hamiltonian in spherical co-ordinates,...
static int CartesianToSpherical(REAL8 qSph[], REAL8 pSph[], const REAL8 qCart[], const REAL8 pCart[])
Perform a co-ordinate transformation from Cartesian to spherical co-ordinates.
static REAL8 CalculateCrossProduct(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 XLALSimIMRSpinEOBInitialConditions(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_v2)
Main function for calculating the spinning EOB initial conditions, following the quasi-spherical init...
static int NormalizeVector(REAL8 a[])
Normalizes the given vector.
static REAL8 XLALCalculateSphHamiltonianDeriv2Hybrid(const int idx1, const int idx2, const REAL8 values[], SpinEOBParams *params)
static REAL8 CalculateDotProduct(const REAL8 a[], const REAL8 b[])
Calculates the dot product of two vectors.
static int CalculateRotationMatrix(gsl_matrix *rotMatrix, gsl_matrix *rotInverse, REAL8 r[], REAL8 v[], REAL8 L[])
Calculate the rotation matrix and its inverse.
static int XLALFindSphericalOrbit(const gsl_vector *x, void *params, gsl_vector *f)
Root function for gsl root finders.
static double GSLSpinHamiltonianDerivWrapper(double x, void *params)
Wrapper for calculating specific derivative of the Hamiltonian in spherical co-ordinates,...
static int SphericalToCartesian(REAL8 qCart[], REAL8 pCart[], const REAL8 qSph[], const REAL8 pSph[])
Performs a co-ordinate transformation from spherical to Cartesian co-ordinates.
#define XLAL_CALLGSL(statement)
#define XLAL_ERROR_REAL8(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_IS_REAL8_FAIL_NAN(val)
The coefficients which are used in calculating the non-quasicircular correction to the EOBNRv2 model.
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.