Commit 5d820966 authored by Roberto Cotesta's avatar Roberto Cotesta Committed by Riccardo Sturani

Change the behaviour of the final spin fit

Under some very rare circumstances, the final spin formula may give
a negative number inside a square root. Previously this was handleded
by having a tolerance and just taking the abs of the argument of the
square root if it was below tolerance and emitting XLAL_EDOM otherwise.
Now, we follow what is done in the final spin python code and simply set
the argument inside the square root to 0 if it is <0. A warning is always
printed when this is done.
parent 8619740b
......@@ -296,18 +296,8 @@ INT4 XLALSimIMREOBFinalMassSpinPrec(
// Guard against NANs inside the fit
if (inside_sqrt < 0)
{
// If the final spin was going to be small, just multiply the argument of the sqrt by -1.
if (sqrt(fabs(inside_sqrt)) * prefactor < 3.0e-2)
{
inside_sqrt *= -1;
}
else
{
// Final spin would have been huge and the argument of the sqrt is negative. Fail.
XLAL_PRINT_ERROR("Final spin fit had a large negative argument of sqrt! Aborting");
XLAL_ERROR(XLAL_EDOM);
}
printf("Warning: in the final spin fit, argument of sqrt was negative!\n");
XLAL_PRINT_WARNING("Warning: in the final spin fit, argument of sqrt was negative. Truncating the spin to 0.\n");
inside_sqrt = 0.0;
}
*finalSpin = prefactor * sqrt(inside_sqrt);
// If somehow we have gone above the Kerr limit, don't.
......
......@@ -2899,10 +2899,8 @@ static int SEOBCalculateSphHarmListNQCCoefficientsV4(
SpinEOBParams *seobParams, /**<< Input: SEOB params */
REAL8Vector *chi1_omegaPeak, /**<< Input: dimensionless spin 1 at peak of
omega in L_N frame */
REAL8Vector *chi2_omegaPeak, /**<< Input: dimensionless spin 2 at peak of
REAL8Vector *chi2_omegaPeak /**<< Input: dimensionless spin 2 at peak of
omega in L_N frame */
UINT4 flag0NQCForOddmModes /**<< Input: flag to put all NQC coeffs to 0 in
the case of a symmetric binary */
) {
/* Masses */
REAL8 m1 = seobParams->eobParams->m1;
......@@ -2958,7 +2956,8 @@ static int SEOBCalculateSphHarmListNQCCoefficientsV4(
REAL8 chi2dotZfinal = chi2_omegaPeak->data[2];
REAL8 chiSfinal = SEOBCalculateChiS(chi1dotZfinal, chi2dotZfinal);
REAL8 chiAfinal = SEOBCalculateChiA(chi1dotZfinal, chi2dotZfinal);
REAL8 q = m1/m2;
//printf("chiA = %.16f\n",chiAfinal);
/* Time elapsed from the start of the dynamics to tPeakOmega */
REAL8 tPeakOmegaFromStartDyn = tPeakOmega - seobdynamics->tVec[0];
......@@ -2977,9 +2976,8 @@ static int SEOBCalculateSphHarmListNQCCoefficientsV4(
EOBNonQCCoeffs *nqcCoeffs = XLALMalloc(sizeof(EOBNonQCCoeffs));
memset(nqcCoeffs, 0, sizeof(EOBNonQCCoeffs));
if (flag0NQCForOddmModes &&
(m % 2 != 0)) { /* In this case, set NQC coeffs to 0 for odd m */
/* In the equal mass equal spins case the odd-m modes are 0, so we set the NQCs to 0 */
if (q<1.005 && (m % 2 != 0) && (fabs(chiAfinal) < 0.15)) { /* In this case, set NQC coeffs to 0 for odd m */
nqcCoeffs->a1 = 0.;
nqcCoeffs->a2 = 0.;
nqcCoeffs->a3 = 0.;
......@@ -5271,16 +5269,9 @@ int XLALSimIMRSpinPrecEOBWaveformAll(
printf("STEP 5) Compute P-frame of modes amp/phase on HiS and compute NQC "
"coefficients\n");
/* Flag setting the NQCs to 0 for odd-m modes in the case of a symmetric
* binary (these modes are then 0 identically) */
INT4 flag0NQCForOddmModes = 0;
if ((eta == 0.25) && (chi1z == chi2z))
flag0NQCForOddmModes = 1;
if (SEOBCalculateSphHarmListNQCCoefficientsV4(
&nqcCoeffsList, modes, nmodes, tPeakOmega, seobdynamicsHiS,
&seobParams, chi1L_tPeakOmega, chi2L_tPeakOmega,
flag0NQCForOddmModes) == XLAL_FAILURE) {
&seobParams, chi1L_tPeakOmega, chi2L_tPeakOmega) == XLAL_FAILURE) {
FREE_ALL
XLALPrintError("XLAL Error - %s: NQC computation failed.\n", __func__);
PRINT_ALL_PARAMS
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment