Add IMRPhenomXO4a

Description

Creating the draft merge request to add the IMRPhenomXO4a model currently undergoing review here: https://git.ligo.org/waveforms/reviews/imrphenomxo4

API Changes and Justification

Backwards Compatible Changes

• This change does not modify any class/function/struct/type definitions in a public C header file or any Python class/function definitions
• This change adds new classes/functions/structs/types to a public C header file or Python module

Backwards Incompatible Changes

• This change modifies an existing class/function/struct/type definition in a public C header file or Python module
• This change removes an existing class/function/struct/type from a public C header file or Python module

If any of the Backwards Incompatible check boxes are ticked please provide a justification why this change is necessary and why it needs to be done in a backwards incompatible way.

Review Status

Review completed. For details, see https://git.ligo.org/waveforms/reviews/imrphenomxo4

lalsimulation/lib/LALSimIMRPhenomXPHM.c

@@ -71,6 +71,7 @@ static int IMRPhenomXPHMTwistUp(
 static int IMRPhenomXPHMTwistUpOneMode(
   const REAL8 Mf,                          /**< Frequency (Mf geometric units) */
   const COMPLEX16 hlmprime,                /**< Underlying aligned-spin IMRPhenomXHM waveform. The loop is with mprime positive, but the mode has to be the negative one for positive frequencies.*/
+  const COMPLEX16 hlmprime_antisym,         /**< antisymmetric waveform in the co-precessing frame */
   IMRPhenomXWaveformStruct *pWF,           /**< IMRPhenomX Waveform Struct */
   IMRPhenomXPrecessionStruct *pPrec,       /**< IMRPhenomXP Precession Struct */
   UINT4  l,                                /**< l index of the (l,m) (non-)precessing mode */
@@ -977,9 +978,6 @@ static int IMRPhenomXPHM_hplushcross(
       {
         if (AntisymmetricWaveform && IMRPhenomXPNRUseTunedAngles)
         {
-        //   XLAL_CHECK((IMRPhenomXPNRUseTunedAngles==1),XLAL_EFUNC,
-        //  "Error: Antisymmetric waveform called with PNR angles turned off\n");
-
           IMRPhenomX_PNR_GenerateAntisymmetricWaveform(antiSym_amp,antiSym_phi,freqs,pWF,pPrec,lalParams);
         }
       }
@@ -1480,6 +1478,8 @@ static int IMRPhenomXPHM_hplushcross(
   XLALDestroyValue(ModeArray);
   XLALDestroyREAL8Sequence(freqs);
 
+  
+
   #if DEBUG == 1
   printf("\n******Leaving IMRPhenomXPHM_hplushcross*****\n");
   #endif
@@ -2567,6 +2567,7 @@ static int IMRPhenomXPHM_OneMode(
 
   /* Set up code for using PNR tuned angles */
   int IMRPhenomXPNRUseTunedAngles = pPrec->IMRPhenomXPNRUseTunedAngles;
+  int AntisymmetricWaveform = pPrec->IMRPhenomXAntisymmetricWaveform;
 
   IMRPhenomX_PNR_angle_spline *hm_angle_spline = NULL;
   REAL8 Mf_RD_22 = pWF->fRING;
@@ -2601,6 +2602,16 @@ static int IMRPhenomXPHM_OneMode(
 
   }
 
+  REAL8Sequence *antiSym_amp = NULL;
+  REAL8Sequence *antiSym_phi = NULL;
+
+  if (AntisymmetricWaveform && IMRPhenomXPNRUseTunedAngles)
+  {
+    antiSym_amp = XLALCreateREAL8Sequence(freqs->length);
+    antiSym_phi = XLALCreateREAL8Sequence(freqs->length);
+  }
+
+
   /***** Loop over non-precessing modes ******/
   for (UINT4 emmprime = 1; emmprime <= ell; emmprime++)
   {
@@ -2665,6 +2676,15 @@ static int IMRPhenomXPHM_OneMode(
       }
     }
 
+    if(ell==2 && emmprime==2)
+      {
+        if (AntisymmetricWaveform && IMRPhenomXPNRUseTunedAngles)
+        {
+          IMRPhenomX_PNR_GenerateAntisymmetricWaveform(antiSym_amp,antiSym_phi,freqs,pWF,pPrec,lalParams);
+        }
+      }
+
+
     if (!(htildelm)){ XLAL_ERROR(XLAL_EFUNC);}
 
 
@@ -2689,6 +2709,7 @@ static int IMRPhenomXPHM_OneMode(
 
      /* Variable to store the non-precessing waveform in one frequency point. */
      COMPLEX16 hlmcoprec;
+     COMPLEX16 hlmcoprec_antiSym =0.0;
 
      if(XLALSimInspiralWaveformParamsLookupPhenomXPHMPrecModes(lalParams) == 1)
      {
@@ -2725,13 +2746,20 @@ static int IMRPhenomXPHM_OneMode(
           XLAL_CHECK(XLAL_SUCCESS == status, XLAL_EFUNC, "Error: IMRPhenomX_PNR_LinearFrequencyMapTransitionFrequencies failed.\n");
         }
       }
-
+      
       for (UINT4 idx = 0; idx < freqs->length; idx++)
       {
           REAL8 Mf = pWF->M_sec * freqs->data[idx];
           hlmcoprec  = htildelm->data->data[idx + offset];  /* Co-precessing waveform */
-          COMPLEX16Sequence *hlm;
-          hlm = XLALCreateCOMPLEX16Sequence(2);
+
+          /***** construct h(2,2) or h(2,-2) in co-precessing frame with antisymmetric contribution *****/
+        if(ell == 2 && emmprime == 2 && AntisymmetricWaveform && IMRPhenomXPNRUseTunedAngles)
+        {
+          hlmcoprec_antiSym = antiSym_amp->data[idx]*cexp(I*antiSym_phi->data[idx]);
+        }
+
+        COMPLEX16Sequence *hlm;
+        hlm = XLALCreateCOMPLEX16Sequence(2);
 
         if(IMRPhenomXPNRUseTunedAngles)
         {
@@ -2743,8 +2771,11 @@ static int IMRPhenomXPHM_OneMode(
           pPrec->gammaPNR = gsl_spline_eval(hm_angle_spline->gamma_spline, f_mapped, hm_angle_spline->gamma_acc);
         }
 
+
+        
         // Twist up
-        IMRPhenomXPHMTwistUpOneMode(Mf, hlmcoprec, pWF, pPrec, ell, emmprime, m, hlm);
+        IMRPhenomXPHMTwistUpOneMode(Mf, hlmcoprec, hlmcoprec_antiSym, pWF, pPrec, ell, emmprime, m, hlm);
+        
 
         (*hlmpos)->data->data[idx + offset] += hlm->data[0];     // Positive frequencies. Freqs do 0, df, 2df, ...., fmax
         (*hlmneg)->data->data[idx + offset] += hlm->data[1];     // Negative frequencies. Freqs do 0, -df, -2df, ...., -fmax
@@ -2787,6 +2818,12 @@ XLALDestroyREAL8Sequence(freqs);
 XLALDestroyCOMPLEX16FrequencySeries(htilde22);
 XLALDestroyValue(ModeArray);
 
+if (AntisymmetricWaveform && IMRPhenomXPNRUseTunedAngles)
+{
+  XLALDestroyREAL8Sequence(antiSym_amp);
+  XLALDestroyREAL8Sequence(antiSym_phi);
+}
+
 
 #if DEBUG == 1
 printf("\n******Leaving IMRPhenomXPHM_OneMode*****\n");
@@ -2805,6 +2842,7 @@ return XLAL_SUCCESS;
 static int IMRPhenomXPHMTwistUpOneMode(
   const REAL8 Mf,                          /**< Frequency (Mf geometric units) */
   const COMPLEX16 hlmprime,                /**< Underlying aligned-spin IMRPhenomXHM waveform. The loop is with mprime positive, but the mode has to be the negative one for positive frequencies.*/
+  const COMPLEX16 hlmprime_antisym,        /**< antisymmetric waveform in the co-precessing frame */
   IMRPhenomXWaveformStruct *pWF,           /**< IMRPhenomX Waveform Struct */
   IMRPhenomXPrecessionStruct *pPrec,       /**< IMRPhenomXP Precession Struct */
   UINT4  l,                                /**< l index of the (l,m) (non-)precessing mode */
@@ -2972,11 +3010,15 @@ static int IMRPhenomXPHMTwistUpOneMode(
     hlmneg += cexp_im_alpha * d44[m+4];
   }
 
+  COMPLEX16 eps_phase_hP_lmprime;
+  COMPLEX16 eps_phase_hP_lmprime_neg;
+
   /* See eqs. E3-E4 in Precessing paper. */
-  COMPLEX16 exp_imprime_epsilon = cexp(mprime*I*epsilon);
-  COMPLEX16 eps_phase_hP_lmprime = 1./exp_imprime_epsilon * hlmprime;
-  COMPLEX16 eps_phase_hP_lmprime_neg = exp_imprime_epsilon * minus1l * conj(hlmprime);
+  COMPLEX16 exp_imprime_epsilon = cexp(mprime*I*epsilon);   
 
+  eps_phase_hP_lmprime = 1./exp_imprime_epsilon * (hlmprime + hlmprime_antisym);
+  eps_phase_hP_lmprime_neg = exp_imprime_epsilon * minus1l * conj(hlmprime - hlmprime_antisym);
+ 
  /* Return h_lminertail */
  (hlminertial)->data[0] = eps_phase_hP_lmprime * hlm;
  (hlminertial)->data[1] = eps_phase_hP_lmprime_neg * hlmneg;