From d20558ad2c91479ff583b9f56999c8b81f5ddf3d Mon Sep 17 00:00:00 2001
From: shrobana <shrobana.ghosh@gmail.com>
Date: Fri, 31 Mar 2023 20:41:16 +0200
Subject: [PATCH] added non-uniform frequency sampling for antisymmetric
waveform and turned antisymmetric wf flag on as default in
XLALSimInspiralChooseFDWaveformSequence; address all comments on code review
---
lalsimulation/lib/LALSimIMRPhenomX.c | 2 +-
.../LALSimIMRPhenomX_AntisymmetricWaveform.c | 107 +++++++++---------
lalsimulation/lib/LALSimInspiral.c | 4 +-
.../lib/LALSimInspiralWaveformCache.c | 2 +-
4 files changed, 59 insertions(+), 56 deletions(-)
diff --git a/lalsimulation/lib/LALSimIMRPhenomX.c b/lalsimulation/lib/LALSimIMRPhenomX.c
index 34cc61a97d..7703665195 100644
--- a/lalsimulation/lib/LALSimIMRPhenomX.c
+++ b/lalsimulation/lib/LALSimIMRPhenomX.c
@@ -1932,7 +1932,7 @@ int IMRPhenomXPGenerateFD(
{
if (!PNRUseTunedAngles)
{
- XLAL_ERROR(XLAL_EFUNC, "antisymmetric waveform cannot be generated without PNR angles");
+ XLAL_ERROR(XLAL_EFUNC, "Error: Antisymmetric waveform generation not supported without PNR angles, please turn on PNR angles to produce waveform with asymmetries in the (2,2) and (2,-2) modes\n");
}
else
{
diff --git a/lalsimulation/lib/LALSimIMRPhenomX_AntisymmetricWaveform.c b/lalsimulation/lib/LALSimIMRPhenomX_AntisymmetricWaveform.c
index bc60695f13..8f6b178b3a 100644
--- a/lalsimulation/lib/LALSimIMRPhenomX_AntisymmetricWaveform.c
+++ b/lalsimulation/lib/LALSimIMRPhenomX_AntisymmetricWaveform.c
@@ -63,8 +63,9 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
*/
/**
* EXTERNAL GENERATE antisymmetric waveform
- * This is an external wrapper to generate the (2,2) antisymmetric waveform,
- * amplitude given the standard inputs given to generate FD waveforms.
+ * This is an external wrapper to generate the (2,2) and (2,-2) antisymmetric waveform,
+ * given the standard inputs given to generate FD waveforms.
+ * Note that at present this is only compatible with the PNR angles (refer arxiv XXXX.YYYYY)
*/
int XLALSimIMRPhenomX_PNR_GenerateAntisymmetricWaveform(
REAL8Sequence **antisymamp, /**< [out] Amplitude of antisymmetric (2,2) waveform */
@@ -97,6 +98,9 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
/* Ensure we have a dictionary */
+ status = IMRPhenomX_Initialize_Powers(&powers_of_lalpi, LAL_PI);
+ XLAL_CHECK(XLAL_SUCCESS == status, status, "Failed to initialize useful powers of LAL_PI.\n");
+
LALDict *lalParams_aux;
if (lalParams == NULL)
{
@@ -107,15 +111,11 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
lalParams_aux = XLALDictDuplicate(lalParams);
}
- int UseAntisymmetric = XLALSimInspiralWaveformParamsLookupPhenomXAntisymmetricWaveform(lalParams_aux);
- int UseTunedAngles = XLALSimInspiralWaveformParamsLookupPhenomXPNRUseTunedAngles(lalParams_aux);
-
- if (!UseAntisymmetric || !UseTunedAngles)
+ XLALSimInspiralWaveformParamsInsertPhenomXAntisymmetricWaveform(lalParams_aux, 1);
+ if(!XLALSimInspiralWaveformParamsLookupPhenomXPNRUseTunedAngles(lalParams_aux))
{
- UseAntisymmetric = 1;
- UseTunedAngles = 1;
- XLALSimInspiralWaveformParamsInsertPhenomXAntisymmetricWaveform(lalParams_aux, UseAntisymmetric);
- XLALSimInspiralWaveformParamsInsertPhenomXPNRUseTunedAngles(lalParams_aux, UseTunedAngles);
+ XLAL_PRINT_WARNING("Warning:Antisymmetric waveform generation currently not supported without PNR angles. Turning on PNR angles ... \n");
+ XLALSimInspiralWaveformParamsInsertPhenomXPNRUseTunedAngles(lalParams_aux,1);
}
/* Map fRef to the start frequency if need be, then make sure it's within the frequency range */
@@ -142,14 +142,23 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
freqs->data[i - iStart] = i * deltaF;
}
- REAL8 f_min_eval = freqs->data[0];
- REAL8 f_max_eval = freqs->data[freqs->length - 1];
-
+ // REAL8 f_min_eval = freqs->data[0];
+ // REAL8 f_max_eval = freqs->data[freqs->length - 1];
+
IMRPhenomXWaveformStruct *pWF;
pWF = XLALMalloc(sizeof(IMRPhenomXWaveformStruct));
- status = IMRPhenomXSetWaveformVariables(pWF, m1_SI, m2_SI, chi1z, chi2z, deltaF, fRef, phiRef, f_min_eval, f_max_eval, distance, inclination, lalParams_aux, DEBUG);
+ status = IMRPhenomXSetWaveformVariables(pWF, m1_SI, m2_SI, chi1z, chi2z, deltaF, fRef, phiRef, f_min, f_max, distance, inclination, lalParams_aux, PHENOMXDEBUG);
XLAL_CHECK(XLAL_SUCCESS == status, XLAL_EFUNC, "Error: IMRPhenomXSetWaveformVariables failed.\n");
+ // IMRPhenomXWaveformStruct *pWF;
+ // pWF = XLALMalloc(sizeof(IMRPhenomXWaveformStruct));
+ // status = IMRPhenomXSetWaveformVariables(pWF, m1_SI, m2_SI, chi1z, chi2z, deltaF, fRef, phiRef, f_min, f_max, distance, inclination, lalParams_aux, PHENOMXDEBUG);
+ // XLAL_CHECK(XLAL_SUCCESS == status, XLAL_EFUNC, "Error: IMRPhenomXSetWaveformVariables failed.\n");
+
+ FILE *fpp = fopen ("/Users/shrobana/lalsuite_forks/lalsuite_review/lalsuite/test_data_second.txt", "w");
+ fprintf(fpp,"%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f \n", m1_SI, m2_SI, chi1z, chi2z, deltaF, fRef, phiRef, f_min, f_max, distance, inclination);
+ fclose(fpp);
+
IMRPhenomXPrecessionStruct *pPrec;
pPrec = XLALMalloc(sizeof(IMRPhenomXPrecessionStruct));
status = IMRPhenomXGetAndSetPrecessionVariables(pWF, pPrec, m1_SI, m2_SI, chi1x, chi1y, chi1z, chi2x, chi2y, chi2z, lalParams_aux, DEBUG);
@@ -163,6 +172,8 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
}
LALFree(pPrec);
XLALDestroyDict(lalParams_aux);
+ XLALDestroyREAL8Sequence(freqs);
+
return XLAL_SUCCESS;
}
@@ -276,6 +287,7 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
/****** antisymmetric amplitude ******/
amp_AS = kappa->data[idx] * amp;
+
if(Mf < MfT)
{
@@ -283,7 +295,7 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
}
else
{
- phi_AS = phi + phi_B0;
+ phi_AS = phi + phi_B0-2;
}
REAL8 Amp0 = pWF->amp0 * pWF->ampNorm;
@@ -292,6 +304,7 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
antisymphase->data[idx] = phi_AS;
}
+
/* Clean up memory allocation */
LALFree(pPhase22);
@@ -299,8 +312,6 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
XLALDestroyREAL8Sequence(kappa);
XLALDestroyREAL8Sequence(alphaPNR);
- // XLALDestroyREAL8Sequence(betaPNR);
- // XLALDestroyREAL8Sequence(gammaPNR);
return XLAL_SUCCESS;
}
@@ -335,7 +346,7 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
XLAL_CHECK(
XLAL_SUCCESS == status,
XLAL_EFUNC,
- "Error: XLALIMRPhenomXPCheckMassesAndSpins failed in XLALSimIMRPhenomX_PNR_GenerateAntisymmetricWaveform.\n");
+ "Error: XLALIMRPhenomXPCheckMassesAndSpins failed in XLALSimIMRPhenomX_PNR_GenerateAntisymmetricAmpRatio.\n");
/* Ensure we have a dictionary */
LALDict *lalParams_aux;
@@ -348,11 +359,11 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
lalParams_aux = XLALDictDuplicate(lalParams);
}
- int UseAntisymmetric = XLALSimInspiralWaveformParamsLookupPhenomXAntisymmetricWaveform(lalParams_aux);
- if (!UseAntisymmetric)
+ XLALSimInspiralWaveformParamsInsertPhenomXAntisymmetricWaveform(lalParams_aux, 1);
+ if(!XLALSimInspiralWaveformParamsLookupPhenomXPNRUseTunedAngles(lalParams_aux))
{
- UseAntisymmetric = 1;
- XLALSimInspiralWaveformParamsInsertPhenomXAntisymmetricWaveform(lalParams_aux, UseAntisymmetric);
+ XLAL_PRINT_WARNING("Warning:Antisymmetric waveform generation currently not supported without PNR angles. Turning on PNR angles ... \n");
+ XLALSimInspiralWaveformParamsInsertPhenomXPNRUseTunedAngles(lalParams_aux,1);
}
/* Map fRef to the start frequency if need be, then make sure it's within the frequency range */
@@ -425,14 +436,6 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
)
{
- int UseAntisymmetric = pPrec->IMRPhenomXAntisymmetricWaveform;
- XLAL_CHECK(
- UseAntisymmetric,
- XLAL_EFUNC,
- "Error: Antisymmetric waveform called without being activated!\n");
-
- REAL8 deltaF = pWF->deltaF;
-
/**** populating parameters of binary from pWF ******/
const double m2 = pWF->m2;
const double M = pWF->Mtot;
@@ -451,42 +454,42 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
double b = b0 + b1*eta + b2*theta + b3*eta*theta;
REAL8 vRD = cbrt (LAL_PI * MfRD );
- // const double kappaPNRD = (21 * vRD*vRD * (1 + delta) * Chi * sin(theta))/( 2 * (42 + vRD * vRD * (55 * eta - 107)) + vRD * vRD * vRD *(84 * LAL_PI - 28 * (1 + delta - eta) * Chi * cos ( theta )));
- const double kappaRD = GetKappa_at_frequency(vRD,delta,Chi,theta,eta,b);/*kappaPNRD * (1 + b * vRD * vRD * vRD * vRD * vRD);*/
-
- if (deltaF > 0.0)
+ const double kappaRD = GetKappa_at_frequency(vRD,delta,Chi,theta,eta,b);
+
+ for (size_t i = 0; i < freqs->length; i++)
{
- for (size_t i = 0; i < freqs->length; i++)
+ REAL8 Mf = XLALSimPhenomUtilsHztoMf(freqs->data[i], M);
+ REAL8 v = cbrt (LAL_PI * Mf );
+ if (Mf < MfRD)
{
- REAL8 Mf = XLALSimPhenomUtilsHztoMf(freqs->data[i], M);
- REAL8 v = cbrt (LAL_PI * Mf * (2.0 / 2.0) );
- if (Mf < MfRD)
- {
- // double kappaPN = (21 * v * v * (1 + delta) * Chi * sin(theta))/(2 * ((42 + v * v * (55 * eta - 107)) + v * v * v *(84 * LAL_PI - 28 * (1 + delta - eta) * Chi * cos ( theta ))));
- kappa->data[i] = GetKappa_at_frequency(v,delta,Chi,theta,eta,b);/*kappaPN * (1 + b * v * v * v * v * v);*/
- }
- else
- {
- kappa->data[i] = kappaRD;
- }
+ kappa->data[i] = GetKappa_at_frequency(v,delta,Chi,theta,eta,b);
+ }
+ else
+ {
+ kappa->data[i] = kappaRD;
}
}
-
+
size_t width = 80;
+ double df = 0.0;
if(width > kappa->length - 1)
{
width = (size_t)floor((double)(kappa->length) / 2.0);
}
size_t half_width = (size_t)floor((double)(width / 2.0));
+
for (size_t id = 0; id < kappa->length-width; id++)
{
double smoothed_ratio = 0.0;
+ double frequency_width = 0.0;
for (size_t j = 0; j < width+1; j++)
{
- smoothed_ratio += kappa->data[id+j];
+ df = freqs->data[id+j+1] - freqs->data[id+j];
+ smoothed_ratio += (kappa->data[id+j]*df);
+ frequency_width += df;
}
- kappa->data[id+half_width] = smoothed_ratio *1.0/((double)(width+1));
+ kappa->data[id+half_width] = smoothed_ratio *1.0/frequency_width;
}
return XLAL_SUCCESS;
@@ -495,8 +498,8 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
double GetKappa_at_frequency(REAL8 v,REAL8 delta,REAL8 Chi,REAL8 theta,REAL8 eta,double b)
{
REAL8 v2 = v * v;
- REAL8 v3 = v * v * v;
- REAL8 v5 = v * v * v * v * v;
+ REAL8 v3 = v2 * v;
+ REAL8 v5 = v3 * v2;
REAL8 kappaPNnum = (21 * v2 * (1 + delta) * Chi * sin(theta));
REAL8 kappaPNden = 2 * (42 + 84 * LAL_PI * v3 + v2 * (55 * eta - 107) - 28 * v3 * (1 + delta - eta) * Chi * cos ( theta ));
REAL8 amp_ratio = kappaPNnum/kappaPNden * (1 + b * v5 );
@@ -579,7 +582,7 @@ IMRPhenomX_UsefulPowers powers_of_lalpi;
*A0 = phi_der_MfT/2 - alpha_der_MfT;
*phi_A0 = pPrec-> alpha_offset;
- *phi_B0 = alpha_MfT - phi_MfT/2 + *A0 * MfT + *phi_A0;
+ *phi_B0 = alpha_MfT - phi_MfT/2 + *A0 * MfT + *phi_A0;
LALFree(alphaParams);
diff --git a/lalsimulation/lib/LALSimInspiral.c b/lalsimulation/lib/LALSimInspiral.c
index 1f2d998578..54499834b9 100644
--- a/lalsimulation/lib/LALSimInspiral.c
+++ b/lalsimulation/lib/LALSimInspiral.c
@@ -2405,7 +2405,7 @@ int XLALSimInspiralChooseFDWaveform(
{
if(!XLALSimInspiralWaveformParamsLookupPhenomXPNRUseTunedAngles(LALparams_aux))
{
- XLAL_ERROR(XLAL_EFUNC,"Error: Antisymmetric waveform generation not supported without PNR angles \n");
+ XLAL_ERROR(XLAL_EFUNC,"Error: Antisymmetric waveform generation not supported without PNR angles, please turn on PNR angles to produce waveform with asymmetries in the (2,2) and (2,-2) modes \n");
}
}
@@ -3751,7 +3751,7 @@ SphHarmFrequencySeries *XLALSimInspiralChooseFDModes(
{
if(!XLALSimInspiralWaveformParamsLookupPhenomXPNRUseTunedAngles(LALparams_aux))
{
- XLAL_ERROR_NULL(XLAL_EFUNC,"Error: Antisymmetric waveform generation not supported without PNR angles \n");
+ XLAL_ERROR_NULL(XLAL_EFUNC,"Error: Antisymmetric waveform generation not supported without PNR angles, please turn on PNR angles to produce waveform with asymmetries in the (2,2) and (2,-2) modes \n");
}
}
diff --git a/lalsimulation/lib/LALSimInspiralWaveformCache.c b/lalsimulation/lib/LALSimInspiralWaveformCache.c
index 44cc8e83a6..8cd2b4dae0 100644
--- a/lalsimulation/lib/LALSimInspiralWaveformCache.c
+++ b/lalsimulation/lib/LALSimInspiralWaveformCache.c
@@ -1520,7 +1520,7 @@ int XLALSimInspiralChooseFDWaveformSequence(
/* Toggle on antisymmetric contributions */
if(!XLALDictContains(LALparams_aux, "AntisymmetricWaveform")){
- XLALSimInspiralWaveformParamsInsertPhenomXAntisymmetricWaveform(LALparams_aux, 0);
+ XLALSimInspiralWaveformParamsInsertPhenomXAntisymmetricWaveform(LALparams_aux, 1);
}
ret = XLALSimIMRPhenomXPHMFrequencySequence(
--
GitLab