Merge branch 'phenompnrtidal-review' into 'master'

Phenompnrtidal review

See merge request !660

lalsimulation/src/LALSimIMR.h

@@ -63,7 +63,8 @@ extern "C" {
 
 typedef enum tagIMRPhenomP_version_type {
  IMRPhenomPv1_V, /**< version 1: based on IMRPhenomC */
- IMRPhenomPv2_V  /**< version 2: based on IMRPhenomD */
+ IMRPhenomPv2_V,  /**< version 2: based on IMRPhenomD */
+ IMRPhenomPv2NRTidal_V /**< version Pv2_NRTidal: based on IMRPhenomPv2; NRTides (https://arxiv.org/pdf/1706.02969.pdf) added before precession */
 } IMRPhenomP_version_type;
 
 /** @} */

lalsimulation/src/LALSimIMRPhenom.c

@@ -1205,6 +1205,10 @@ static gsl_matrix *XLALSimIMRPhenomBProjectExtrinsicParam(
  *    based on IMRPhenomD
  *    (should be used instead of IMRPhenomP,
  *    unless there are good reasons not to).
+ *  * IMRPhenomPv2_NRTidal models precessing binaries,
+ *    adds NRTides (https://arxiv.org/pdf/1706.02969.pdf) 
+ *    to IMRPhenomD phasing and twists the waveform 
+ *    to generate the corresponding precessing waveform.
  *
  * @review IMRPhenomB routines reviewed by Frank Ohme, P. Ajith, Alex Nitz
  * and Riccardo Sturani. The review concluded with git hash
@@ -1217,7 +1221,10 @@ static gsl_matrix *XLALSimIMRPhenomBProjectExtrinsicParam(
  * db16d17013531cd10451c7d0c6906972ce731866 (Oct/Nov 2015).
  *
  * @review original IMRPhenomP not reviewed, nor going to be.
- * IMRPhenomPv2 currently under review (Dec 2015).
+ * IMRPhenomPv2 reviewed by Capano, Pürrer, Bohe et al. Conludeded
+ * with git hash 1354291cf6a897995a04cd12dce42b7acaca7b34 (May 2016)
+ *
+ * @review IMRPhenomPv2_NRTidal currently under review (2018).
  * @{
  */
 

lalsimulation/src/LALSimIMRPhenomD_internals.c

@@ -1135,6 +1135,7 @@ static double PhiInsAnsatzInt(double Mf, UsefulPowers * powers_of_Mf, PhiInsPref
   phasing += prefactors->minus_third / powers_of_Mf->third;
   phasing += prefactors->minus_two_thirds / powers_of_Mf->two_thirds;
   phasing += prefactors->minus_one / Mf;
+  phasing += prefactors->minus_four_thirds / powers_of_Mf->four_thirds;
   phasing += prefactors->minus_five_thirds / powers_of_Mf->five_thirds; // * v^0
 
   // Now add higher order terms that were calibrated for PhenomD
@@ -1166,6 +1167,7 @@ static int init_phi_ins_prefactors(PhiInsPrefactors * prefactors, IMRPhenomDPhas
 	prefactors->minus_third = pn->v[4] / powers_of_pi.third;
 	prefactors->minus_two_thirds = pn->v[3] / powers_of_pi.two_thirds;
 	prefactors->minus_one = pn->v[2] / Pi;
+        prefactors->minus_four_thirds = pn->v[1] / powers_of_pi.four_thirds;
 	prefactors->minus_five_thirds = pn->v[0] / powers_of_pi.five_thirds; // * v^0
 
   // higher order terms that were calibrated for PhenomD

lalsimulation/src/LALSimIMRPhenomD_internals.h

@@ -301,6 +301,7 @@ typedef struct tagPhiInsPrefactors
 	double minus_third;
 	double minus_two_thirds;
 	double minus_one;
+        double minus_four_thirds;
 	double minus_five_thirds;
 } PhiInsPrefactors;
 

lalsimulation/src/LALSimIMRPhenomP.h

@@ -97,7 +97,7 @@ static int PhenomPCore(
    * If deltaF > 0, the frequency points given in freqs are uniformly spaced with
    * spacing deltaF. Otherwise, the frequency points are spaced non-uniformly.
    * Then we will use deltaF = 0 to create the frequency series we return. */
-  IMRPhenomP_version_type IMRPhenomP_version, /**< IMRPhenomPv1 uses IMRPhenomC, IMRPhenomPv2 uses IMRPhenomD */
+  IMRPhenomP_version_type IMRPhenomP_version, /**< IMRPhenomPv1 uses IMRPhenomC, IMRPhenomPv2 uses IMRPhenomD, IMRPhenomPv2_NRTidal is a tidal version of IMRPhenomPv2 */
   LALDict *extraParams /**< linked list containing the extra testing GR parameters */
 );
 
@@ -105,9 +105,6 @@ static int PhenomPCore(
 static int PhenomPCoreOneFrequency(
   const REAL8 fHz,                        /**< Frequency (Hz) */
   const REAL8 eta,                        /**< Symmetric mass ratio */
-  const REAL8 chi1_l,                     /**< Dimensionless aligned spin on companion 1 */
-  const REAL8 chi2_l,                     /**< Dimensionless aligned spin on companion 2 */
-  const REAL8 chip,                       /**< Dimensionless spin in the orbital plane */
   const REAL8 distance,                   /**< Distance of source (m) */
   const REAL8 M,                          /**< Total mass (Solar masses) */
   const REAL8 phic,                       /**< Orbital phase at the peak of the underlying non precessing model (rad) */
@@ -115,18 +112,30 @@ static int PhenomPCoreOneFrequency(
   IMRPhenomDPhaseCoefficients *pPhi,      /**< Internal IMRPhenomD phase coefficients */
   BBHPhenomCParams *PCparams,             /**< Internal PhenomC parameters */
   PNPhasingSeries *PNparams,              /**< PN inspiral phase coefficients */
-  NNLOanglecoeffs *angcoeffs,             /**< Struct with PN coeffs for the NNLO angles */
-  SpinWeightedSphericalHarmonic_l2 *Y2m,  /**< Struct of l=2 spherical harmonics of spin weight -2 */
-  const REAL8 alphaoffset,                /**< f_ref dependent offset for alpha angle (azimuthal precession angle) */
-  const REAL8 epsilonoffset,              /**< f_ref dependent offset for epsilon angle */
-  COMPLEX16 *hp,                          /**< Output: tilde h_+ */
-  COMPLEX16 *hc,                          /**< Output: tilde h_+ */
+  COMPLEX16 *hPhenom,                     /**< Output: IMRPhenom waveform (before precession) */
   REAL8 *phasing,                         /**< Output: overall phasing */
-  const UINT4 IMRPhenomP_version,         /**< Version number: 1 uses IMRPhenomC, 2 uses IMRPhenomD */
+  const UINT4 IMRPhenomP_version,         /**< Version number: 1 uses IMRPhenomC, 2 uses IMRPhenomD, NRTidal uses IMRPhenomPv2 with the NRTidal framework */
   AmpInsPrefactors *amp_prefactors,       /**< pre-calculated (cached for saving runtime) coefficients for amplitude. See LALSimIMRPhenomD_internals.c*/
   PhiInsPrefactors *phi_prefactors        /**< pre-calculated (cached for saving runtime) coefficients for phase. See LALSimIMRPhenomD_internals.*/
 );
 
+static int PhenomPCoreTwistUp(
+  const REAL8 fHz,                            /**< Frequency (Hz) */
+  COMPLEX16 hPhenom,                    /**< [in] IMRPhenom waveform (before precession) */
+  const REAL8 eta,                            /**< Symmetric mass ratio */
+  const REAL8 chi1_l,                         /**< Dimensionless aligned spin on companion 1 */
+  const REAL8 chi2_l,                         /**< Dimensionless aligned spin on companion 2 */
+  const REAL8 chip,                           /**< Dimensionless spin in the orbital plane */
+  const REAL8 M,                              /**< Total mass (Solar masses) */
+  NNLOanglecoeffs *angcoeffs,                 /**< Struct with PN coeffs for the NNLO angles */
+  SpinWeightedSphericalHarmonic_l2 *Y2m,      /**< Struct of l=2 spherical harmonics of spin weight -2 */
+  const REAL8 alphaoffset,                    /**< f_ref dependent offset for alpha angle (azimuthal precession angle) */
+  const REAL8 epsilonoffset,                  /**< f_ref dependent offset for epsilon angle */
+  COMPLEX16 *hp,                              /**< [out] plus polarization \f$\tilde h_+\f$ */
+  COMPLEX16 *hc,                              /**< [out] cross polarization \f$\tilde h_x\f$ */
+  IMRPhenomP_version_type IMRPhenomP_version  /**< IMRPhenomP(v1) uses IMRPhenomC, IMRPhenomPv2 uses IMRPhenomD, IMRPhenomPv2_NRTidal uses NRTidal framework with IMRPhenomPv2 */
+);
+
 /* Simple 2PN version of L, without any spin terms expressed as a function of v */
 static REAL8 L2PNR(
   const REAL8 v,   /**< Cubic root of (Pi * Frequency (geometric)) */
@@ -190,4 +199,20 @@ static REAL8 FinalSpinBarausse2009(  /* Barausse & Rezzolla, Astrophys.J.Lett.70
 static bool approximately_equal(REAL8 x, REAL8 y, REAL8 epsilon);
 static void nudge(REAL8 *x, REAL8 X, REAL8 epsilon);
 
+static int PhenomPCoreOneFrequency_withTides(
+  const REAL8 fHz,                            /**< Frequency (Hz) */
+  const REAL8 ampTidal,                      /**< tidal amplitude at a frequency sample; planck window */
+  COMPLEX16 phaseTidal,                      /**< tidal phasing at a frequency sample from NRTidal infrastructure*/
+  const REAL8 distance,                       /**< Distance of source (m) */
+  const REAL8 M,                              /**< Total mass (Solar masses) */
+  const REAL8 phic,                           /**< Orbital phase at the peak of the underlying non precessing model (rad) */
+  IMRPhenomDAmplitudeCoefficients *pAmp,      /**< Internal IMRPhenomD amplitude coefficients */
+  IMRPhenomDPhaseCoefficients *pPhi,          /**< Internal IMRPhenomD phase coefficients */
+  PNPhasingSeries *PNparams,                  /**< PN inspiral phase coefficients */
+  COMPLEX16 *hPhenom,                         /**< [out] IMRPhenom waveform (before precession) */
+  REAL8 *phasing,                             /**< [out] overall phasing */
+  AmpInsPrefactors *amp_prefactors,           /**< pre-calculated (cached for saving runtime) coefficients for amplitude. See LALSimIMRPhenomD_internals.c*/
+  PhiInsPrefactors *phi_prefactors            /**< pre-calculated (cached for saving runtime) coefficients for phase. See LALSimIMRPhenomD_internals.*/
+);
+
 #endif	// of #ifndef _LALSIM_IMR_PHENOMP_H

lalsimulation/src/LALSimInspiral.c

@@ -27,6 +27,7 @@
 
 #include <lal/SphericalHarmonics.h>
 #include <lal/LALSimInspiral.h>
+#include <lal/LALSimInspiralEOS.h>
 #include <lal/LALSimIMR.h>
 #include <lal/LALSimSphHarmMode.h>
 #include <lal/LALConstants.h>
@@ -151,6 +152,7 @@ static const char *lalSimulationApproximantNames[] = {
     INITIALIZE_NAME(IMRPhenomD_NRTidal),
     INITIALIZE_NAME(IMRPhenomP),
     INITIALIZE_NAME(IMRPhenomPv2),
+    INITIALIZE_NAME(IMRPhenomPv2_NRTidal),
     INITIALIZE_NAME(IMRPhenomFC),
     INITIALIZE_NAME(TaylorEt),
     INITIALIZE_NAME(TaylorT4),
@@ -271,6 +273,8 @@ static double fixReferenceFrequency(const double f_ref, const double f_min, cons
         case SpinTaylorF2:
         case IMRPhenomP:
         case IMRPhenomPv2:
+        case IMRPhenomPv2_NRTidal:
+            return f_min;
         case NRSur4d2s:
             return f_min;
         default:
@@ -758,6 +762,10 @@ int XLALSimInspiralChooseTDWaveform(
 	    ret = XLALSimInspiralTDFromFD(hplus, hcross, m1, m2, S1x, S1y, S1z, S2x, S2y, S2z, distance, inclination, phiRef, longAscNodes, eccentricity, meanPerAno, deltaT, f_min, f_ref, LALparams, approximant);
 	    break;
 
+        case IMRPhenomPv2_NRTidal:
+            ret = XLALSimInspiralTDFromFD(hplus, hcross, m1, m2, S1x, S1y, S1z, S2x, S2y, S2z, distance, inclination, phiRef, longAscNodes, eccentricity, meanPerAno, deltaT, f_min, f_ref, LALparams, approximant);
+           break;
+
         case PhenSpinTaylorRD:
             /* Waveform-specific sanity checks */
             if( !checkTidesZero(lambda1, lambda2) )
@@ -963,12 +971,13 @@ int XLALSimInspiralChooseFDWaveform(
     unsigned int j;
     REAL8 pfac, cfac;
     INT4 phiRefAtEnd;
+    int amplitudeO = XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(LALparams);
+    int phaseO = XLALSimInspiralWaveformParamsLookupPNPhaseOrder(LALparams);
+
     REAL8 quadparam1 = 1.+XLALSimInspiralWaveformParamsLookupdQuadMon1(LALparams);
     REAL8 quadparam2 = 1.+XLALSimInspiralWaveformParamsLookupdQuadMon2(LALparams);
     REAL8 lambda1 = XLALSimInspiralWaveformParamsLookupTidalLambda1(LALparams);
     REAL8 lambda2 = XLALSimInspiralWaveformParamsLookupTidalLambda2(LALparams);
-    int amplitudeO = XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(LALparams);
-    int phaseO =XLALSimInspiralWaveformParamsLookupPNPhaseOrder(LALparams);
 
     /* Support variables for precessing wfs*/
     REAL8 spin1x,spin1y,spin1z;
@@ -1518,6 +1527,33 @@ int XLALSimInspiralChooseFDWaveform(
             }
             break;
 
+        case IMRPhenomPv2_NRTidal:
+            /* Waveform-specific sanity checks */
+            if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
+                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);/* Default is LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L : z-axis along direction of orbital angular momentum. */
+            if(!XLALSimInspiralWaveformParamsModesChoiceIsDefault(          /* Default is (2,2) or l=2 modes. */LALparams) )
+                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
+            /* Tranform to model parameters */
+            if(f_ref==0.0)
+                f_ref = f_min; /* Default reference frequency is minimum frequency */
+            XLALSimIMRPhenomPCalculateModelParametersFromSourceFrame(
+                &chi1_l, &chi2_l, &chip, &thetaJN, &alpha0, &phi_aligned, &zeta_polariz,
+                m1, m2, f_ref, phiRef, inclination,
+                S1x, S1y, S1z,
+                S2x, S2y, S2z, IMRPhenomPv2NRTidal_V);
+            /* Call the waveform driver routine */
+            ret = XLALSimIMRPhenomP(hptilde, hctilde,
+              chi1_l, chi2_l, chip, thetaJN,
+              m1, m2, distance, alpha0, phi_aligned, deltaF, f_min, f_max, f_ref, IMRPhenomPv2NRTidal_V, LALparams);
+            if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
+            for (UINT4 idx=0;idx<(*hptilde)->data->length;idx++) {
+                PhPpolp=(*hptilde)->data->data[idx];
+                PhPpolc=(*hctilde)->data->data[idx];
+                (*hptilde)->data->data[idx] =cos(2.*zeta_polariz)*PhPpolp+sin(2.*zeta_polariz)*PhPpolc;
+                (*hctilde)->data->data[idx]=cos(2.*zeta_polariz)*PhPpolc-sin(2.*zeta_polariz)*PhPpolp;
+            }
+            break;
+
         case SpinTaylorT4Fourier:
             /* Waveform-specific sanity checks */
             if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
@@ -4626,6 +4662,7 @@ int XLALSimInspiralImplementedTDApproximants(
         case IMRPhenomC:
 	case IMRPhenomD:
 	case IMRPhenomPv2:
+        case IMRPhenomPv2_NRTidal:
         case PhenSpinTaylorRD:
         case SEOBNRv1:
         case SpinDominatedWf:
@@ -4667,6 +4704,7 @@ int XLALSimInspiralImplementedFDApproximants(
         case IMRPhenomD_NRTidal:
         case IMRPhenomP:
         case IMRPhenomPv2:
+        case IMRPhenomPv2_NRTidal:
         case EOBNRv2_ROM:
         case EOBNRv2HM_ROM:
         case SEOBNRv1_ROM_EffectiveSpin:
@@ -5070,6 +5108,7 @@ int XLALSimInspiralGetSpinSupportFromApproximant(Approximant approx){
     case SpinTaylorT3:
     case IMRPhenomP:
     case IMRPhenomPv2:
+    case IMRPhenomPv2_NRTidal:
     case SpinTaylorT2Fourier:
     case SpinTaylorT4Fourier:
     case SpinDominatedWf:
@@ -5231,6 +5270,7 @@ int XLALSimInspiralApproximantAcceptTestGRParams(Approximant approx){
     case IMRPhenomD:
     case IMRPhenomP:
     case IMRPhenomPv2:
+    case IMRPhenomPv2_NRTidal:
       testGR_accept=LAL_SIM_INSPIRAL_TESTGR_PARAMS;
       break;
     default:
@@ -7048,6 +7088,34 @@ int XLALSimInspiralChooseFDWaveformOLD(
     return ret;
 }
 
+/** 
+ * if you do NOT provide a quadparam[1,2] term and you DO provide
+ * lamdba[1,2] then we calculate quad-mono term using universal relations
+ * quadparam[1,2]_UR: Quadrupole-Monopole parameter computed using
+ * universal relations (UR) 
+ */
+
+int XLALSimInspiralSetQuadMonParamsFromLambdas(
+       LALDict *LALparams /**< LAL dictionary containing accessory parameters */
+       )
+{
+    REAL8 quadparam1 = XLALSimInspiralWaveformParamsLookupdQuadMon1(LALparams);
+    REAL8 quadparam2 = XLALSimInspiralWaveformParamsLookupdQuadMon2(LALparams);
+    REAL8 lambda1 = XLALSimInspiralWaveformParamsLookupTidalLambda1(LALparams);
+    REAL8 lambda2 = XLALSimInspiralWaveformParamsLookupTidalLambda2(LALparams);
+
+    if ((lambda1 > 0) && (quadparam1 == 0)) {
+        REAL8 quadparam1_UR = XLALSimInspiralEOSQfromLambda(lambda1);
+        XLALSimInspiralWaveformParamsInsertdQuadMon1(LALparams, quadparam1_UR - 1.);
+    }
+
+    if ((lambda2 > 0) && (quadparam2 == 0)) {
+        REAL8 quadparam2_UR = XLALSimInspiralEOSQfromLambda(lambda2);
+        XLALSimInspiralWaveformParamsInsertdQuadMon2(LALparams, quadparam2_UR - 1.);
+    }
+    return XLAL_SUCCESS;
+}
+
 /** @} */
 
 /** @} */

lalsimulation/src/LALSimInspiral.h

@@ -374,6 +374,7 @@ typedef enum tagApproximant {
                          * @remarks Implemented in lalsimulation (frequency domain).  */
    IMRPhenomPv2,		/**< Frequency domain (generic spins) inspiral-merger-ringdown templates of Hannam et al., arXiv:1308.3271 [gr-qc]. Based on IMRPhenomD, arXiv:1508.07250 and arXiv:1508.07253.
                          * @remarks Implemented in lalsimulation (frequency domain).  */
+   IMRPhenomPv2_NRTidal, /**< Frequency domain tidal version of IMRPhenomPv2, using NRTidal framework from arXiv:1706.02969 */
    IMRPhenomFC,		/**< Frequency domain (non-precessing spins) inspiral-merger-ringdown templates of Santamaria et al [Santamaria:2010yb] with phenomenological coefficients defined in the Table II of [Santamaria:2010yb]
                          * @attention Not implemented in lalsimulation.*/
    TaylorEt,		/**< UNDOCUMENTED
@@ -463,7 +464,7 @@ typedef enum tagTestGRaccept {
  * Structure for passing around PN phasing coefficients.
  * For use with the TaylorF2 waveform.
  */
-#define PN_PHASING_SERIES_MAX_ORDER 12
+#define PN_PHASING_SERIES_MAX_ORDER 15
 typedef struct tagPNPhasingSeries
 {
     REAL8 v[PN_PHASING_SERIES_MAX_ORDER+1];
@@ -820,6 +821,8 @@ int XLALSimInspiralREAL8WaveTaper(REAL8Vector *signalvec, LALSimInspiralApplyTap
 
 int XLALSimInspiralTEOBResumROM(REAL8TimeSeries **hPlus, REAL8TimeSeries **hCross, REAL8 phiRef, REAL8 deltaT, REAL8 fLow, REAL8 fRef, REAL8 distance, REAL8 inclination, REAL8 m1SI, REAL8 m2SI, REAL8 lambda1, REAL8 lambda2);
 
+int XLALSimInspiralSetQuadMonParamsFromLambdas(LALDict *LALpars);
+
 
 
 #if 0

lalsimulation/src/LALSimInspiralEOS.c

+/*
+ *  Copyright (C) 2012 Walter Del Pozzo, Tjonnie Li, Michalis Agathos
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with with program; see the file COPYING. If not, write to the
+ *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+ *  MA  02111-1307  USA
+ */
+ #include <stdlib.h> 
+ #include <stdio.h>
+ #include <string.h>
+ #include <lal/LALSimInspiralEOS.h>
+ #include <lal/LALSimInspiral.h>
+
+ LALEquationOfState XLALSimEOSfromString(char eos_name[])
+ {
+   LALEquationOfState eos;
+   if (!strcmp("MS1",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_MS1;
+   else if (!strcmp("H4",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_H4;
+   else if (!strcmp("SQM3",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_SQM3;
+   else if (!strcmp("MPA1",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_MPA1;
+   else if (!strcmp("GNH3",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_GNH3;
+   else if (!strcmp("A",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_A;
+   else if (!strcmp("AU",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_AU;
+   else if (!strcmp("FPS",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_FPS;
+   else if (!strcmp("APR",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_APR;
+   else if (!strcmp("UU",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_UU;
+   else if (!strcmp("L",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_L;
+   else if (!strcmp("PP",eos_name)) eos = LAL_SIM_INSPIRAL_EOS_NONE;
+   else
+   {
+    XLALPrintError( "XLAL Error - %s: Equation of state %s not recognized.", __func__, eos_name);
+    XLAL_ERROR(XLAL_EINVAL);
+   }
+   return eos;
+ }
+
+REAL8 XLALSimInspiralEOSLambda(LALEquationOfState eos_type, REAL8 m_intr_msun){/** this must be fed the INTRINSIC mass */
+
+    /* this is fed the intrinsic masses and then computes the value of \Lambda(m) See Hinderer et al ( http://arxiv.org/abs/0911.3535 ) for details of the EOSes*/
+    /* \Lambda(m) is in units of s^5 */
+    REAL8 lambda=0.;
+//  printf("EOS number: %d\n", eos_type);
+//  printf("mass: %e\n", m_intr_msun);
+    switch (eos_type)
+    {
+        case LAL_SIM_INSPIRAL_EOS_NONE:
+            lambda = 0.0;
+        break;
+    // MS1
+        case LAL_SIM_INSPIRAL_EOS_MS1:
+           // printf("Using EOS MS1\n");
+            lambda = 2.755956E-24*(2.19296 + 20.0273*m_intr_msun - 17.9443*m_intr_msun*m_intr_msun 
+            + 5.75129*m_intr_msun*m_intr_msun*m_intr_msun - 0.699095*m_intr_msun*m_intr_msun*m_intr_msun*m_intr_msun);
+        break;
+    // H4
+        case LAL_SIM_INSPIRAL_EOS_H4:
+            lambda = 2.755956E-24*(0.743351 + 15.8917*m_intr_msun - 14.7348*m_intr_msun*m_intr_msun 
+            + 5.32863*m_intr_msun*m_intr_msun*m_intr_msun - 0.942625*m_intr_msun*m_intr_msun*m_intr_msun*m_intr_msun);
+        break; 
+    // SQM3
+        case LAL_SIM_INSPIRAL_EOS_SQM3:
+            lambda = 2.755956E-24*(-5.55858 + 20.8977*m_intr_msun - 20.5583*m_intr_msun*m_intr_msun 
+            + 9.55465*m_intr_msun*m_intr_msun*m_intr_msun - 1.84933*m_intr_msun*m_intr_msun*m_intr_msun*m_intr_msun);
+        break;
+    // MPA1
+    case LAL_SIM_INSPIRAL_EOS_MPA1:
+        lambda = 2.755956E-24*(0.276761 + 7.26925*m_intr_msun - 5.72102*m_intr_msun*m_intr_msun
+        + 1.51347*m_intr_msun*m_intr_msun*m_intr_msun - 0.152181*m_intr_msun*m_intr_msun*m_intr_msun*m_intr_msun);
+        break;
+    // GNH3
+    case LAL_SIM_INSPIRAL_EOS_GNH3:
+        lambda = 2.755956E-24*(7.80715 + 0.683549*m_intr_msun + 1.21351*m_intr_msun*m_intr_msun
+        - 3.50234*m_intr_msun*m_intr_msun*m_intr_msun + 0.894662*m_intr_msun*m_intr_msun*m_intr_msun*m_intr_msun);
+        break;
+    default:
+        lambda = 0.0;
+        break;
+    }
+//  printf("calculated love number: %e\n", lambda);
+    if (lambda<0.0) return 0.0;
+    else return lambda;
+}
+
+REAL8 XLALLambdaQuadratic(REAL8 c0, REAL8 c1, REAL8 c2, REAL8 mass) {
+    mass = mass*LAL_MTSUN_SI;
+    // [LAMBDA0] = SEC^5; [LAMBDA1] = SEC^4; [LAMBDA2] = SEC^3
+    REAL8 lambda = 1.0E-23*c0 + 1.0E-18*(mass-1.4*LAL_MTSUN_SI)*c1 + 1.0E-13*(mass-1.4*LAL_MTSUN_SI)*(mass-1.4*LAL_MTSUN_SI)*c2;
+    lambda = (lambda > 0.0) ? lambda : 0.0;
+    return lambda;
+}
+ 
+
+REAL8 XLALSimInspiralEOSQfromLambda(REAL8 lambda) {
+    /* Quadrupole-monopole parameter calculated from love number;
+       see http://arxiv.org/abs/1303.1528 */
+    REAL8 q, loglam;
+    REAL8 tolerance = 5E-1;
+    if(lambda<tolerance) { //printf("Love number is (nearly) zero; cannot compute QM parameter. Setting to 1.0 (BH value).\n");
+                      q = 1.0; } 
+    else {
+    loglam = log(lambda);
+    q =  0.194 + 0.0936*loglam + 0.0474*loglam*loglam;
+    q -= 0.00421*loglam*loglam*loglam;
+    q += 0.000123*loglam*loglam*loglam*loglam;
+    q = exp(q);
+    }
+
+//  printf("%e %e\n", l, q); // Testing numerical results from these functions.
+
+    return q;
+
+}
+
+REAL8 XLALSimInspiralEOSqmparameter(LALEquationOfState eos_type, REAL8 m_intr_msun){
+  
+  REAL8 q = 0.0 ;
+  REAL8 m = m_intr_msun ;
+  REAL8 m2 = m*m ;
+  REAL8 m3 = m2*m ;
+  
+  switch (eos_type) {
+  /*  */
+  case LAL_SIM_INSPIRAL_EOS_A:
+    q = -6.41414141*m3 + 30.70779221*m2 - 53.37417027*m + 35.62253247 ;
+    break;
+  /*  */
+  case LAL_SIM_INSPIRAL_EOS_AU:
+    q = -6.18686869*m3 + 30.15909091*m2 - 52.87806638*m + 35.86616883 ;
+    break;
+  /*  */
+  case LAL_SIM_INSPIRAL_EOS_FPS:
+    q = -3.86363636*m3 + 21.03030303*m2 - 42.19448052*m + 32.83722944 ;
+    break;
+  /*  */
+  case LAL_SIM_INSPIRAL_EOS_APR:
+    q = -10.55555556*m3 + 49.52380952*m2 - 82.77063492*m + 53.02428571 ;
+    break;
+  /*  */
+  case LAL_SIM_INSPIRAL_EOS_UU:
+    q = -8.03030303*m3 + 37.61363636*m2 - 63.48733766*m + 41.75080087 ;
+    break;
+  /*  */
+  case LAL_SIM_INSPIRAL_EOS_L:
+    q = -6.59090909*m3 + 33.67424242*m2 - 63.77034632*m + 48.98073593 ;
+    break;
+  case LAL_SIM_INSPIRAL_EOS_NONE:
+    q = 1.0 ;
+    break;
+
+  default:
+    q = 1.0 ;
+    break ;
+  }
+  
+  if (q < 1.0) {
+    q = 1.0;
+  }
+  
+  return q ;
+}
+
+/**
+ * This function estimates the radius of a NS of a given mass and
+ * tidal deformability parameter, based on the "I-Love-Q forever" relation
+ * of Maselli et al, arXiv:1304.2052v1.
+ * To be used for masses within [1.2,2]M_sun, and preferably not for strange
+ * quark stars (since the relation is calibrated for this mass range and for 
+ * the EoS APR4, MS1, H4).
+ * For a BH, (lambda=0) it returns the Schwarzschild radius.
+ * The arguments are:
+ * m_intr_msun              the intrinsic mass in solar masses
+ * barlambda                the dimensionless tidal deformability (lambda/m^5)
+ * The return value is the radius in meters.
+ */
+
+REAL8 XLALSimInspiralNSRadiusOfLambdaM(REAL8 m_intr_msun, REAL8 barlambda){
+
+  REAL8 loglambda;
+  REAL8 compactness, radius ;
+  REAL8 tolerance = 1E-15;
+
+  /* Check for sign of lambda */
+  if ( barlambda <= tolerance && barlambda >= 0.0 ) {
+  /* This is a black hole */
+    compactness = 0.5;
+  }
+  else if ( barlambda > tolerance ) {
+  loglambda = log(barlambda);
+  /* Calculate compactness according to arXiv:1304.2052v1 */
+  compactness = 0.371 - 0.0391*loglambda + 0.001056*loglambda*loglambda;
+  }
+  else {
+    XLALPrintError( "XLAL Error - %s: Tidal deformability is negative. Only positive values are acceptable.", __func__);
+    XLAL_ERROR_REAL8(XLAL_EDOM);
+  }
+
+  /* Check that radius is larger than Schwarzschild radius */
+  if ( compactness > 0.5 ) {
+    XLALPrintWarning( "XLAL Warning - %s: Neutron Star is calculated to have compactness larger than a black hole (C = %f, lambda = %f, m = %f).\n Setting C=0.5 ...", __func__, compactness, barlambda, m_intr_msun);
+    compactness = 0.5;
+  }
+
+  if ( compactness < 0.0 ) {
+    XLALPrintError( "XLAL Error - %s: Neutron Star is calculated to have negative compactness (C = %f, lambda = %f, m = %f).", __func__, compactness, barlambda, m_intr_msun);
+    XLAL_ERROR_REAL8(XLAL_ERANGE);
+  }
+
+  radius = LAL_MRSUN_SI * m_intr_msun / compactness;
+
+ return radius;
+
+}
+
+
+/**
+ * This function estimates the radius for a binary of given masses and
+ * tidal deformability parameters.
+ * It uses XLALSimInspiralNSRadiusOfLambdaM() to calculate radii (see above).
+ * The arguments are:
+ * m1_intr, m2_intr                      the intrinsic masses in solar masses
+ * barlambda1, barlambda2                the dimensionless tidal deformabilities (lambda_i/m_i^5)
+ * The return value is the GW contact frequency in Hz.
+ */
+
+REAL8 XLALSimInspiralContactFrequency(REAL8 m1_intr, REAL8 barlambda1, REAL8 m2_intr, REAL8 barlambda2){
+
+  REAL8 r1, r2, rtot, mtot, f_gw_contact;
+
+  /* Calculate radii for the two components */
+  r1 = XLALSimInspiralNSRadiusOfLambdaM(m1_intr, barlambda1);
+  r2 = XLALSimInspiralNSRadiusOfLambdaM(m2_intr, barlambda2);
+
+  rtot = (r1 + r2)/LAL_C_SI;                             // Orbital distance in seconds
+  mtot = (m1_intr + m2_intr)*LAL_MTSUN_SI;               // Total mass in seconds
+
+  /* Calculate the GW contact frequency */
+  f_gw_contact = sqrt(mtot/(rtot*rtot*rtot))/LAL_PI;
+  if ( f_gw_contact < 0.0 ) {
+    XLALPrintError( "XLAL Error - %s: Contact frequency is calculated to be negative  (fcontact = %f)", __func__, f_gw_contact);
+    XLAL_ERROR_REAL8(XLAL_ERANGE);
+  }
+
+  return f_gw_contact;
+
+}

lalsimulation/src/LALSimInspiralEOS.h

+/*
+ *  Copyright (C) 2012 Walter Del Pozzo, Tjonnie Li, Michalis Agathos
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with with program; see the file COPYING. If not, write to the
+ *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+ *  MA  02111-1307  USA
+ */
+
+#ifndef _LALSIMINSPIRALEOS_H  /* Double-include protection. */
+#define _LALSIMINSPIRALEOS_H
+
+#ifdef  __cplusplus   /* C++ protection. */
+extern "C" {
+#endif
+