Commit 3d3763b7 authored by HyungWon Lee's avatar HyungWon Lee

move LALSimInspiralTaylorF2Ecc.c file reference into Makefile.am for automatic build

add LALSimInspiralTaylorF2Ecc.c noinst_HEADERS list in Makefile.am

add corrections for switch fall thourgh problem for GNUC >= 7

add sanity check for f_ecc adn eccentricity

phase calculation for TaylorF2Ecc

add TaylorF2Ecc test

limit the f_max <= f_ISCO in TaylorF2Ecc wavefunction code

correct the sanity check error message for wrong eccentricity

print warning message for wrong eccentricity specified instead aborting the code

correct warning message

check for the case of f_ref = 0, in this case set f_ecc = f_min

use XLAL_PRINT_WARNING macro to print warning
parent ced0d1d3
......@@ -1079,6 +1079,14 @@ int XLALSimInspiralChooseFDWaveform(
ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
if( !LALSimInspiralEccentricityIsCorrect(eccentricity, LALparams) ){
/* we set f_ecc to be f_ref for wrong eccentricity specification after long discussion.
ABORT_WRONG_ECCENTRICITY(LALparams); */
REAL8 f_ecc = f_ref;
if( f_ecc == 0 ) f_ecc = f_min;
XLALSimInspiralWaveformParamsInsertEccentricityFreq(LALparams, f_ecc);
XLAL_PRINT_WARNING("Warning... The reference frequency for eccentricity was set as default value(%f). This might be not optimal case for you.\n", f_ecc);
}
/* Call the waveform driver routine */
ret = XLALSimInspiralTaylorF2Ecc(hptilde, phiRef, deltaF, m1, m2,
......
......@@ -637,14 +637,18 @@ int XLALHGimriGenerator(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8
/* TaylorF2 functions */
/* in module LALSimInspiralTaylorF2.c */
int XLALSimInspiralTaylorF2CoreEcc(COMPLEX16FrequencySeries **htilde, const REAL8Sequence *freqs, const REAL8 phi_ref, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 S1z, const REAL8 S2z, const REAL8 f_ref, const REAL8 shft, const REAL8 r, const REAL8 eccentricity, LALDict *LALparams);
int XLALSimInspiralTaylorF2Ecc(COMPLEX16FrequencySeries **htilde, const REAL8 phi_ref, const REAL8 deltaF, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 S1z, const REAL8 S2z, const REAL8 fStart, const REAL8 fEnd, const REAL8 f_ref, const REAL8 r, const REAL8 eccentricity, LALDict *LALparams);
int XLALSimInspiralTaylorF2AlignedPhasing(PNPhasingSeries **pfa, const REAL8 m1, const REAL8 m2, const REAL8 chi1, const REAL8 chi2, LALDict *extraPars);
int XLALSimInspiralTaylorF2AlignedPhasingArray(REAL8Vector **phasingvals, REAL8Vector mass1, REAL8Vector mass2, REAL8Vector chi1, REAL8Vector chi2, REAL8Vector lambda1, REAL8Vector lambda2, REAL8Vector dquadmon1, REAL8Vector dquadmon2);
int XLALSimInspiralTaylorF2Core(COMPLEX16FrequencySeries **htilde, const REAL8Sequence *freqs, const REAL8 phi_ref, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 f_ref, const REAL8 shft, const REAL8 r, LALDict *LALparams, PNPhasingSeries *pfaP);
int XLALSimInspiralTaylorF2(COMPLEX16FrequencySeries **htilde, const REAL8 phi_ref, const REAL8 deltaF, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 S1z, const REAL8 S2z, const REAL8 fStart, const REAL8 fEnd, const REAL8 f_ref, const REAL8 r, LALDict *LALpars);
/* TaylorF2Ecc functions */
/* in module LALSimInspiralTaylorF2Ecc.c */
int XLALSimInspiralTaylorF2CoreEcc(COMPLEX16FrequencySeries **htilde, const REAL8Sequence *freqs, const REAL8 phi_ref, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 f_ref, const REAL8 shft, const REAL8 r, const REAL8 eccentricity, LALDict *LALparams, PNPhasingSeries *pfaP);
int XLALSimInspiralTaylorF2Ecc(COMPLEX16FrequencySeries **htilde, const REAL8 phi_ref, const REAL8 deltaF, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 S1z, const REAL8 S2z, const REAL8 fStart, const REAL8 fEnd, const REAL8 f_ref, const REAL8 r, const REAL8 eccentricity, LALDict *LALparams);
int LALSimInspiralEccentricityIsCorrect(REAL8 eccentricity, LALDict *params);
/* TaylorF2NLPhase functions */
/* in module LALSimInspiralTaylorF2NLTides.c */
......
......@@ -64,13 +64,12 @@ int XLALSimInspiralTaylorF2CoreEcc(
const REAL8 phi_ref, /**< reference orbital phase (rad) */
const REAL8 m1_SI, /**< mass of companion 1 (kg) */
const REAL8 m2_SI, /**< mass of companion 2 (kg) */
const REAL8 S1z, /**< z component of the spin of companion 1 */
const REAL8 S2z, /**< z component of the spin of companion 2 */
const REAL8 f_ref, /**< Reference GW frequency (Hz) - if 0 reference point is coalescence */
const REAL8 shft, /**< time shift to be applied to frequency-domain phase (sec)*/
const REAL8 r, /**< distance of source (m) */
const REAL8 eccentricity, /**< eccentricity effect control < 0 : no eccentricity effect */
LALDict *p /**< Linked list containing the extra parameters >**/
LALDict *p, /**< Linked list containing the extra parameters >**/
PNPhasingSeries *pfaP /**< Phasing coefficients >**/
)
{
......@@ -83,8 +82,6 @@ int XLALSimInspiralTaylorF2CoreEcc(
const REAL8 m_sec = m * LAL_MTSUN_SI; /* total mass in seconds */
const REAL8 eta = m1 * m2 / (m * m);
const REAL8 piM = LAL_PI * m_sec;
const REAL8 m1OverM = m1 / m;
const REAL8 m2OverM = m2 / m;
REAL8 amp0;
size_t i;
COMPLEX16 *data = NULL;
......@@ -105,8 +102,7 @@ int XLALSimInspiralTaylorF2CoreEcc(
}
/* phasing coefficients */
PNPhasingSeries pfa;
XLALSimInspiralPNPhasing_F2(&pfa, m1, m2, S1z, S2z, S1z*S1z, S2z*S2z, S1z*S2z, p);
PNPhasingSeries pfa = *pfaP;
REAL8 pfaN = 0.; REAL8 pfa1 = 0.;
REAL8 pfa2 = 0.; REAL8 pfa3 = 0.; REAL8 pfa4 = 0.;
......@@ -120,20 +116,41 @@ int XLALSimInspiralTaylorF2CoreEcc(
case -1:
case 7:
pfa7 = pfa.v[7];
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 6:
pfa6 = pfa.v[6];
pfl6 = pfa.vlogv[6];
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 5:
pfa5 = pfa.v[5];
pfl5 = pfa.vlogv[5];
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 4:
pfa4 = pfa.v[4];
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 3:
pfa3 = pfa.v[3];
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 2:
pfa2 = pfa.v[2];
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 1:
pfa1 = pfa.v[1];
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 0:
pfaN = pfa.v[0];
break;
......@@ -169,15 +186,19 @@ int XLALSimInspiralTaylorF2CoreEcc(
*/
REAL8 pft10 = 0.;
REAL8 pft12 = 0.;
REAL8 lambda1=XLALSimInspiralWaveformParamsLookupTidalLambda1(p);
REAL8 lambda2=XLALSimInspiralWaveformParamsLookupTidalLambda2(p);
switch( XLALSimInspiralWaveformParamsLookupPNTidalOrder(p) )
{
case LAL_SIM_INSPIRAL_TIDAL_ORDER_ALL:
case LAL_SIM_INSPIRAL_TIDAL_ORDER_6PN:
pft12 = pfaN * (lambda1*XLALSimInspiralTaylorF2Phasing_12PNTidalCoeff(m1OverM) + lambda2*XLALSimInspiralTaylorF2Phasing_12PNTidalCoeff(m2OverM) );
pft12 = pfa.v[12];
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case LAL_SIM_INSPIRAL_TIDAL_ORDER_5PN:
pft10 = pfaN * ( lambda1*XLALSimInspiralTaylorF2Phasing_10PNTidalCoeff(m1OverM) + lambda2*XLALSimInspiralTaylorF2Phasing_10PNTidalCoeff(m2OverM) );
pft10 = pfa.v[10];
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case LAL_SIM_INSPIRAL_TIDAL_ORDER_0PN:
break;
default:
......@@ -314,19 +335,37 @@ int XLALSimInspiralTaylorF2CoreEcc(
{
case 7:
flux += FTa7 * v7;
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 6:
flux += (FTa6 + FTl6*logv) * v6;
dEnergy += dETa3 * v6;
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 5:
flux += FTa5 * v5;
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 4:
flux += FTa4 * v4;
dEnergy += dETa2 * v4;
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 3:
flux += FTa3 * v3;
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case 2:
flux += FTa2 * v2;
dEnergy += dETa1 * v2;
#if __GNUC__ >= 7
__attribute__ ((fallthrough));
#endif
case -1: /* Default to no SPA amplitude corrections */
case 0:
flux += 1.;
......@@ -409,7 +448,7 @@ int XLALSimInspiralTaylorF2Ecc(
if (r <= 0) XLAL_ERROR(XLAL_EDOM);
/* allocate htilde */
if ( fEnd == 0. ) // End at ISCO
if ( fEnd == 0. || fEnd > fISCO) // End at ISCO
f_max = fISCO;
else // End at user-specified freq.
f_max = fEnd;
......@@ -435,8 +474,12 @@ int XLALSimInspiralTaylorF2Ecc(
for (i = iStart; i < n; i++) {
freqs->data[i-iStart] = i * deltaF;
}
/* phasing coefficients */
PNPhasingSeries pfa;
XLALSimInspiralPNPhasing_F2(&pfa, m1, m2, S1z, S2z, S1z*S1z, S2z*S2z, S1z*S2z, p);
ret = XLALSimInspiralTaylorF2CoreEcc(&htilde, freqs, phi_ref, m1_SI, m2_SI,
S1z, S2z, f_ref, shft, r, eccentricity, p);
f_ref, shft, r, eccentricity, p, &pfa);
XLALDestroyREAL8Sequence(freqs);
......@@ -445,5 +488,22 @@ int XLALSimInspiralTaylorF2Ecc(
return ret;
}
/**
* Returns true if f_ecc and eccentricity were set correctly, check for the case of non-zero eccentricity without setting f_ecc
* value; returns false otherwise.
* Pointed out by Riccaro Sturani
*/
int LALSimInspiralEccentricityIsCorrect(REAL8 eccentricity, LALDict *params)
{
/**
* returns true for the case
* 1) eccentricity = 0.0, no eccentricity effect is required
* 2) eccentricy > 0.0 and f_ecc > 0.0, f_ecc and eccentricity were set intentionally.Default value of f_ecc = -1.0
*/
REAL8 f_ecc = 0;
f_ecc = XLALSimInspiralWaveformParamsLookupEccentricityFreq(params); /** get f_ecc */
return ( eccentricity == 0.0 || (eccentricity > 0.0 && f_ecc > 0.0));
}
/** @} */
/** @} */
......@@ -81,7 +81,7 @@ DEFINE_INSERT_FUNC(TidalOctupolarFMode2, REAL8, "TidalOctupolarFMode2", 0)
DEFINE_INSERT_FUNC(dQuadMon1, REAL8, "dQuadMon1", 0)
DEFINE_INSERT_FUNC(dQuadMon2, REAL8, "dQuadMon2", 0)
DEFINE_INSERT_FUNC(Redshift, REAL8, "redshift", 0)
DEFINE_INSERT_FUNC(EccentricityFreq, REAL8, "f_ecc", 10)
DEFINE_INSERT_FUNC(EccentricityFreq, REAL8, "f_ecc", -1.0)
DEFINE_INSERT_FUNC(NonGRPhi1, REAL8, "phi1", 0)
DEFINE_INSERT_FUNC(NonGRPhi2, REAL8, "phi2", 0)
......@@ -181,7 +181,7 @@ DEFINE_LOOKUP_FUNC(TidalOctupolarFMode2, REAL8, "TidalOctupolarFMode2", 0)
DEFINE_LOOKUP_FUNC(dQuadMon1, REAL8, "dQuadMon1", 0)
DEFINE_LOOKUP_FUNC(dQuadMon2, REAL8, "dQuadMon2", 0)
DEFINE_LOOKUP_FUNC(Redshift, REAL8, "redshift", 0)
DEFINE_LOOKUP_FUNC(EccentricityFreq, REAL8, "f_ecc", 10)
DEFINE_LOOKUP_FUNC(EccentricityFreq, REAL8, "f_ecc", -1.0)
DEFINE_LOOKUP_FUNC(NonGRPhi1, REAL8, "phi1", 0)
DEFINE_LOOKUP_FUNC(NonGRPhi2, REAL8, "phi2", 0)
......@@ -274,7 +274,7 @@ DEFINE_ISDEFAULT_FUNC(TidalOctupolarFMode2, REAL8, "TidalOctupolarFMode2", 0)
DEFINE_ISDEFAULT_FUNC(dQuadMon1, REAL8, "dQuadMon1", 0)
DEFINE_ISDEFAULT_FUNC(dQuadMon2, REAL8, "dQuadMon2", 0)
DEFINE_ISDEFAULT_FUNC(Redshift, REAL8, "redshift", 0)
DEFINE_ISDEFAULT_FUNC(EccentricityFreq, REAL8, "f_ecc", 10)
DEFINE_ISDEFAULT_FUNC(EccentricityFreq, REAL8, "f_ecc", -1.0)
DEFINE_ISDEFAULT_FUNC(NonGRPhi1, REAL8, "phi1", 0)
DEFINE_ISDEFAULT_FUNC(NonGRPhi2, REAL8, "phi2", 0)
......
......@@ -185,6 +185,7 @@ noinst_HEADERS = \
LALSimIMRSpinEOBInitialConditions.c \
LALSimIMRSpinEOBInitialConditionsPrec.c \
LALSimInspiralPNCoefficients.c \
LALSimInspiralTaylorF2Ecc.c \
LALSimInspiraldEnergyFlux.c \
LALSimNeutronStarEOSPiecewisePolytrope.c \
LALSimNeutronStarEOSTabular.c \
......@@ -220,8 +221,8 @@ liblalsimulation_la_SOURCES = \
LALSimIMRPhenomD.c \
LALSimIMRPhenomD_NRTidal.c \
LALSimIMRPhenomP.c \
LALSimIMRPhenomInternalUtils.c \
LALSimIMRPhenomUtils.c \
LALSimIMRPhenomInternalUtils.c \
LALSimIMRPhenomUtils.c \
LALSimIMRPSpinInspiralRD.c \
LALSimPhenSpinRingDown.c \
LALSimInspiralWaveformFlags.c \
......
......@@ -279,6 +279,18 @@ XLAL_ERROR_NULL(XLAL_EINVAL);\
XLAL_ERROR(XLAL_EINVAL);\
} while (0)
/*
* Macro procedure for aborting if not correct eccentricity parameters are given
* if eccentricity > 0 ,then f_ecc > 0 also.
* default f_ecc = -1.0, when do not specify
*/
#define ABORT_WRONG_ECCENTRICITY(LALparams)\
do {\
XLALDestroyDict(LALparams);\
XLALPrintError("XLAL Error - %s: Not correct eccentricity poarameters are given, eccentricity is greater than zero without specifying f_ecc. This approximant does not support this case.\n", __func__);\
XLAL_ERROR(XLAL_EINVAL);\
} while (0)
/* Internal utility macro to check all spin components are zero
returns 1 if all spins zero, otherwise returns 0 */
#define checkSpinsZero(s1x, s1y, s1z, s2x, s2y, s2z) \
......
......@@ -131,7 +131,7 @@ do
done
# Test FD aligned-spin, tidal approximants
approx="TaylorF2 TaylorF2RedSpinTidal"
approx="TaylorF2 TaylorF2RedSpinTidal TaylorF2Ecc"
for i in $approx;
do
......
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