LALInferenceTemplate.c

/*
 *  LALInferenceTemplate.c: Bayesian Followup, template calls to LAL
 *  template functions. Temporary GeneratePPN
 *
 *  Copyright (C) 2009 Ilya Mandel, Vivien Raymond, Christian Roever,
 *  Marc van der Sluys, John Veitch, Will M. Farr and Salvatore Vitale
 *
 *
 *  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., 51 Franklin Street, Fifth Floor, Boston,
 *  MA  02110-1301  USA
 */

#include <stdio.h>
#include <stdlib.h>
#include <lal/LALInspiral.h>
#include <lal/SeqFactories.h>
#include <lal/TimeSeries.h>
#include <lal/FrequencySeries.h>
#include <lal/Date.h>
#include <lal/VectorOps.h>
#include <lal/TimeFreqFFT.h>
#include <lal/GenerateInspiral.h>
#include <lal/GenerateInspRing.h>
#include <lal/LALStatusMacros.h>
#include <lal/LALInference.h>
#include <lal/XLALError.h>
#include <lal/LIGOMetadataRingdownUtils.h>
#include <lal/LALSimInspiral.h>
#include <lal/LALInferenceTemplate.h>
#include <lal/LALInferenceMultibanding.h>
#include <lal/LALSimNeutronStar.h>

/* LIB imports*/
#include <lal/LALInferenceBurstRoutines.h>


#define PROGRAM_NAME "LALInferenceTemplate.c"
#define CVS_ID_STRING "$Id$"
#define CVS_REVISION "$Revision$"
#define CVS_SOURCE "$Source$"
#define CVS_DATE "$Date$"
#define CVS_NAME_STRING "$Name$"

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

/* Max amplitude orders found in LALSimulation (not accessible from outside of LALSim) */
#define MAX_NONPRECESSING_AMP_PN_ORDER 6
#define MAX_PRECESSING_AMP_PN_ORDER 3

#define Pi_p2 9.8696044010893586188344909998761511
#define Pi_p2by3 2.1450293971110256000774441009412356
#define log4 1.3862943611198906188344642429163531

static void q2masses(double mc, double q, double *m1, double *m2);

/* list of testing GR parameters to be passed to the waveform */
/* the first batch of parameters dchis through dsigmas refer to the parameterised tests for generation (TIGER) while the parameters log10lambda_eff through LIV_A_sign are testing coefficients for the parameterised tests for propagation using a deformed dispersion relation (LIV); new parameters may be added at the end for  readability although the order of parameters in this list does not matter */


const char list_extra_parameters[42][16] = {"dchi0","dchi1","dchi2","dchi3","dchi4","dchi5","dchi5l","dchi6","dchi6l","dchi7","aPPE","alphaPPE","bPPE","betaPPE","betaStep","fStep","dxi1","dxi2","dxi3","dxi4","dxi5","dxi6","dalpha1","dalpha2","dalpha3","dalpha4","dalpha5","dbeta1","dbeta2","dbeta3","dsigma1","dsigma2","dsigma3","dsigma4","log10lambda_eff","lambda_eff","nonGR_alpha","LIV_A_sign","dQuadMon1","dQuadMon2","dQuadMonS","dQuadMonA"};


const UINT4 N_extra_params = 42;


/* Return the quadrupole moment of a neutron star given its lambda
 * We use the relations defined here. https://arxiv.org/pdf/1302.4499.pdf.
 * Note that the convention we use is that:
 * \f$\mathrm{dquadmon} = \bar{Q} - 1.\f$
 * Where \f$\bar{Q}\f$ (dimensionless) is the reduced quadrupole moment.
 */
static REAL8 dquadmon_from_lambda(REAL8 lambdav);
static REAL8 dquadmon_from_lambda(REAL8 lambdav)
{
    
    double ll = log(lambdav);
    double ai = .194, bi = .0936, ci = 0.0474, di = -4.21e-3, ei = 1.23e-4;
    double ln_quad_moment = ai + bi*ll + ci*ll*ll + di*pow(ll,3.0) + ei*pow(ll,4.0);
    return(exp(ln_quad_moment) - 1.0);
}

static int InterpolateWaveform(REAL8Vector *freqs, COMPLEX16FrequencySeries *src, COMPLEX16FrequencySeries *dest);
static int InterpolateWaveform(REAL8Vector *freqs, COMPLEX16FrequencySeries *src, COMPLEX16FrequencySeries *dest)
{
  REAL8 deltaF = dest->deltaF;
  UINT4 j=ceil(freqs->data[0] / deltaF);
  COMPLEX16 *d=dest->data->data;
  memset(d, 0, sizeof(*(d))*j);

  /* Loop over reduced frequency set */
  for(UINT4 i=0;i<freqs->length-1;i++)
  {
    double startpsi = carg(src->data->data[i]);
    double startamp = cabs(src->data->data[i]);
    double endpsi = carg(src->data->data[i+1]);
    double endamp = cabs(src->data->data[i+1]);

    double startf=freqs->data[i];
    double endf=freqs->data[i+1];

    double df = endf - startf; /* Big freq step */

    /* linear interpolation setup */
    double dpsi = (endpsi - startpsi);

    /* Catch case where phase wraps around */
    /* NOTE: If changing this check that waveforms are not corrupted
     * at high frequencies when dpsi/df can go slightly -ve without
     * the phase wrapping around (e.g. TF2 1.4-1.4 srate=4096)
     */
    if (dpsi/df<-LAL_PI ) {dpsi+=LAL_TWOPI;}

    double dpsidf = dpsi/df;
    double dampdf = (endamp - startamp)/df;

    double damp = dampdf *deltaF;

    const double dim = sin(dpsidf*deltaF);
    const double dre = 2.0*sin(dpsidf*deltaF*0.5)*sin(dpsidf*deltaF*0.5);

    /* Loop variables */
    double newRe,newIm,f,re,im,a;
    for(f=j*deltaF,
	re = cos(startpsi), im = sin(startpsi),
        a = startamp;

        f<endf;

        j++, f+=deltaF,
        newRe = re - dre*re-dim*im,
        newIm = im + re*dim-dre*im,
        re=newRe, im = newIm,
        a += damp )
    {
      d[j] = a * (re + I*im);
    }
  }
  memset(&(d[j]), 0, sizeof(d[j])*(dest->data->length - j));
  return 0;
}


void LALInferenceTemplateNullFreqdomain(LALInferenceModel *model)
/**********************************************/
/* returns a frequency-domain 'null' template */
/* (all zeroes, implying no signal present).  */
/**********************************************/
{
  UINT4 i;
  if ((model->freqhPlus==NULL) || (model->freqhCross==NULL)) {
    XLALPrintError(" ERROR in templateNullFreqdomain(): encountered unallocated 'freqModelhPlus/-Cross'.\n");
    XLAL_ERROR_VOID(XLAL_EFAULT);
  }
  for (i=0; i<model->freqhPlus->data->length; ++i){
    model->freqhPlus->data->data[i] = 0.0;
    model->freqhCross->data->data[i] = 0.0;
  }
  model->domain = LAL_SIM_DOMAIN_FREQUENCY;
  return;
}



void LALInferenceTemplateNullTimedomain(LALInferenceModel *model)
/*********************************************/
/* returns a time-domain 'null' template     */
/* (all zeroes, implying no signal present). */
/*********************************************/
{
  UINT4 i;
  if ((model->timehPlus==NULL) || (model->timehCross==NULL)) {
    XLALPrintError(" ERROR in templateNullTimedomain(): encountered unallocated 'timeModelhPlus/-Cross'.\n");
    XLAL_ERROR_VOID(XLAL_EFAULT);
  }
  for (i=0; i<model->timehPlus->data->length; ++i){
    model->timehPlus->data->data[i]  = 0.0;
    model->timehCross->data->data[i] = 0.0;
  }
  model->domain = LAL_SIM_DOMAIN_TIME;
  return;
}



/* ============ LAL template wrapper function: ========== */


static void mc2masses(double mc, double eta, double *m1, double *m2);

static void mc2masses(double mc, double eta, double *m1, double *m2)
/*  Compute individual companion masses (m1, m2)   */
/*  for given chirp mass (m_c) & mass ratio (eta)  */
/*  (note: m1 >= m2).                              */
{
  double root = sqrt(0.25-eta);
  double fraction = (0.5+root) / (0.5-root);
  *m2 = mc * (pow(1+fraction,0.2) / pow(fraction,0.6));
  *m1 = mc * (pow(1+1.0/fraction,0.2) / pow(1.0/fraction,0.6));
  return;
}

static void q2masses(double mc, double q, double *m1, double *m2)
/*  Compute individual companion masses (m1, m2)   */
/*  for given chirp mass (m_c) & asymmetric mass   */
/*  ratio (q).  note: q = m2/m1, where m1 >= m2    */
{
  *m1 = mc * pow(q, -3.0/5.0) * pow(q+1, 1.0/5.0);
  *m2 = (*m1) * q;
  return;
}

void LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence(LALInferenceModel *model){
/*************************************************************************************************************************/
  Approximant approximant = (Approximant) 0;

  int ret=0;
  INT4 errnum=0;

  model->roq->hptildeLinear=NULL, model->roq->hctildeLinear=NULL;
  model->roq->hptildeQuadratic=NULL, model->roq->hctildeQuadratic=NULL;
  REAL8 mc;
  REAL8 phi0, m1, m2, distance, inclination;

  REAL8 *m1_p,*m2_p;

  if (LALInferenceCheckVariable(model->params, "LAL_APPROXIMANT"))
    approximant = *(Approximant*) LALInferenceGetVariable(model->params, "LAL_APPROXIMANT");
  else {
    XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_APPROXIMANT\" parameter not provided!\n");
    XLAL_ERROR_VOID(XLAL_EDATA);
  }

  if (LALInferenceCheckVariable(model->params, "LAL_PNORDER"))
    XLALSimInspiralWaveformParamsInsertPNPhaseOrder(model->LALpars, *(INT4 *) LALInferenceGetVariable(model->params, "LAL_PNORDER"));
  else {
    XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_PNORDER\" parameter not provided!\n");
    XLAL_ERROR_VOID(XLAL_EDATA);
  }

  /* Explicitly set the default amplitude order if one is not specified.
   *   This serves two purposes:
   *     1) The default behavior of the code won't change unexpectedly due to changes in LALSimulation.
   *     2) We need to know the amplitude order in order to set the starting frequency of the waveform properly. */
  if (LALInferenceCheckVariable(model->params, "LAL_AMPORDER"))
    XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars, *(INT4 *) LALInferenceGetVariable(model->params, "LAL_AMPORDER"));
  else
    if (!XLALSimInspiralWaveformParamsPNAmplitudeOrderIsDefault(model->LALpars))
      XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,-1);

  REAL8 f_ref = 100.0;
  if (LALInferenceCheckVariable(model->params, "f_ref")) f_ref = *(REAL8 *)LALInferenceGetVariable(model->params, "f_ref");

  REAL8 fTemp = f_ref;
  if(LALInferenceCheckVariable(model->params,"chirpmass"))
    {
      mc  = *(REAL8*) LALInferenceGetVariable(model->params, "chirpmass");
      if (LALInferenceCheckVariable(model->params,"q")) {
	REAL8 q = *(REAL8 *)LALInferenceGetVariable(model->params,"q");
	q2masses(mc, q, &m1, &m2);
      } else {
	REAL8 eta = *(REAL8*) LALInferenceGetVariable(model->params, "eta");
	mc2masses(mc, eta, &m1, &m2);
      }
    }
  else if((m1_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass1")) && (m2_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass2")))
    {
      m1=*m1_p;
      m2=*m2_p;
    }
  else
    {
      fprintf(stderr,"No mass parameters found!");
      exit(0);
    }
  if(LALInferenceCheckVariable(model->params,"logdistance"))
  {
      distance = exp(LALInferenceGetREAL8Variable(model->params, "logdistance"))* LAL_PC_SI * 1.0e6;        /* distance (1 Mpc) in units of metres */
  }
  else
  {
      distance = LAL_PC_SI * 1.0e6; /* 1Mpc default */
  }

  phi0 = LALInferenceGetREAL8Variable(model->params, "phase"); /* START phase as per lalsimulation convention, radians*/
  /* Zenith angle between J and N in radians. Also known as inclination angle when spins are aligned */
  REAL8 thetaJN = acos(LALInferenceGetREAL8Variable(model->params, "costheta_jn"));     /* zenith angle between J and N in radians */

  /* ==== SPINS ==== */
  /* We will default to spinless signal and then add in the spin components if required */
  /* If there are non-aligned spins, we must convert between the System Frame coordinates
   * and the cartestian coordinates */

  /* The cartesian spin coordinates (default 0), as passed to LALSimulation */
  REAL8 spin1x = 0.0;
  REAL8 spin1y = 0.0;
  REAL8 spin1z = 0.0;
  REAL8 spin2x = 0.0;
  REAL8 spin2y = 0.0;
  REAL8 spin2z = 0.0;

  /* System frame coordinates as used for jump proposals */
  REAL8 a_spin1 = 0.0;  /* Magnitude of spin1 */
  REAL8 a_spin2 = 0.0;  /* Magnitude of spin2 */
  REAL8 phiJL  = 0.0;  /* azimuthal angle of L_N on its cone about J radians */
  REAL8 tilt1   = 0.0;  /* zenith angle between S1 and LNhat in radians */
  REAL8 tilt2   = 0.0;  /* zenith angle between S2 and LNhat in radians */
  REAL8 phi12   = 0.0;  /* difference in azimuthal angle btwn S1, S2 in radians */

  /* Now check if we have spin amplitudes */
  if(LALInferenceCheckVariable(model->params, "a_spin1"))    a_spin1   = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin1");
  if(LALInferenceCheckVariable(model->params, "a_spin2"))    a_spin2   = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin2");

  /* Check if we have spin angles too */
  if(LALInferenceCheckVariable(model->params, "phi_jl"))
      phiJL = LALInferenceGetREAL8Variable(model->params, "phi_jl");
  if(LALInferenceCheckVariable(model->params, "tilt_spin1"))
      tilt1 = LALInferenceGetREAL8Variable(model->params, "tilt_spin1");
  if(LALInferenceCheckVariable(model->params, "tilt_spin2"))
      tilt2 = LALInferenceGetREAL8Variable(model->params, "tilt_spin2");
  if(LALInferenceCheckVariable(model->params, "phi12"))
      phi12 = LALInferenceGetREAL8Variable(model->params, "phi12");

  /* If we have tilt angles zero, then the spins are aligned and we just set the z component */
  /* However, if the waveform supports precession then we still need to get the right coordinate components */
  if(tilt1==0.0 && tilt2==0.0)
  {
      spin1z=a_spin1;
      spin2z=a_spin2;
      inclination = thetaJN; /* Inclination angle is just thetaJN */
  }
  else
  {   /* Template is not aligned-spin only. */
      /* Set all the other spin components according to the angles we received above */
      /* The transformation function doesn't know fLow, so f_ref==0 isn't interpretted as a request to use the starting frequency for reference. */
      XLAL_TRY(ret=XLALSimInspiralTransformPrecessingNewInitialConditions(
                    &inclination, &spin1x, &spin1y, &spin1z, &spin2x, &spin2y, &spin2z,
                    thetaJN, phiJL, tilt1, tilt2, phi12, a_spin1, a_spin2, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, fTemp, phi0), errnum);
      if (ret == XLAL_FAILURE)
      {
        XLALPrintError(" ERROR in XLALSimInspiralTransformPrecessingNewInitialConditions(): error converting angles. errnum=%d\n",errnum );
        return;
      }
  }


/* ==== Spin induced quadrupole moment PARAMETERS ==== */
 
 if(LALInferenceCheckVariable(model->params, "dQuadMonS")&&LALInferenceCheckVariable(model->params, "dQuadMonA")){
    REAL8 dQuadMon1=0.;
    REAL8 dQuadMon2=0.;
    REAL8 dQuadMonS = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonS");
    REAL8 dQuadMonA = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonA");
    LALInferencedQuadMonSdQuadMonA(dQuadMonS,dQuadMonA,&dQuadMon1,&dQuadMon2);
    XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
    XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
        fprintf(stdout,"Both dQuadMonS and dQaudMonA are  sampled");
    fprintf(stdout,"dQM1: %e, dQM2: %e, dQMS: %e, dQMA: %e \n",dQuadMon1,dQuadMon2,dQuadMonS,dQuadMonA);
  }
 else  if(LALInferenceCheckVariable(model->params, "dQuadMon1")&&LALInferenceCheckVariable(model->params, "dQuadMon2")){
    REAL8 dQuadMon1 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon1");
    REAL8 dQuadMon2 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon2");
    XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
    XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
        fprintf(stdout,"Both dQuadMon1 and dQaudMon2 are  sampled");
    fprintf(stdout,"dQM1: %e, dQM2: %e \n",dQuadMon1,dQuadMon2);
  }



  /* ==== TIDAL PARAMETERS ==== */
  if(LALInferenceCheckVariable(model->params, "lambda1"))
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda1"));
  if(LALInferenceCheckVariable(model->params, "lambda2"))
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda2"));
  REAL8 lambdaT = 0.;
  REAL8 dLambdaT = 0.;
  REAL8 sym_mass_ratio_eta = 0.;

  if(LALInferenceCheckVariable(model->params, "lambdaT")&&LALInferenceCheckVariable(model->params, "dLambdaT")){
    REAL8 lambda1=0.;
    REAL8 lambda2=0.;
    lambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "lambdaT");
    dLambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "dLambdaT");
    sym_mass_ratio_eta = m1*m2/((m1+m2)*(m1+m2));
    LALInferenceLambdaTsEta2Lambdas(lambdaT,dLambdaT,sym_mass_ratio_eta,&lambda1,&lambda2);
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars,lambda1);
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars,lambda2);
  }

  /* ==== 4-PIECE POLYTROPE EOS PARAMETERS ==== */
  REAL8 logp1 = 0.;
  REAL8 gamma1 = 0.;
  REAL8 gamma2 = 0.;
  REAL8 gamma3 = 0.;
  if(LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "gamma1")&&LALInferenceCheckVariable(model->params, "gamma2")&&LALInferenceCheckVariable(model->params, "gamma3")){
    REAL8 lambda1 = 0.;
    REAL8 lambda2 = 0.;
    logp1 = *(REAL8*) LALInferenceGetVariable(model->params, "logp1");
    gamma1 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma1");
    gamma2 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma2");
    gamma3 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma3");
    // Find lambda1,2(m1,2|eos)
    LALInferenceLogp1GammasMasses2Lambdas(logp1,gamma1,gamma2,gamma3,m1,m2,&lambda1,&lambda2);
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
  }

  /* ==== SPECTRAL DECOMPOSITION PARAMETERS ==== */
  REAL8 SDgamma0 = 0.;
  REAL8 SDgamma1 = 0.;
  REAL8 SDgamma2 = 0.;
  REAL8 SDgamma3 = 0.;
  /* Checks for 4 spectral parameters */
  if(!LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "SDgamma0")&&LALInferenceCheckVariable(model->params, "SDgamma1")&&LALInferenceCheckVariable(model->params, "SDgamma2")&&LALInferenceCheckVariable(model->params,"SDgamma3")){
    REAL8 lambda1 = 0.;
    REAL8 lambda2 = 0.;
    SDgamma0 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma0");
    SDgamma1 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma1");
    SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
    SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
    REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
    LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
  }

  /* ==== BINARY_LOVE PARAMETERS ==== */
  if(LALInferenceCheckVariable(model->params, "lambdaS")){
    REAL8 lambda1=0.;
    REAL8 lambda2=0.;
    LALInferenceBinaryLove(model->params, &lambda1, &lambda2);
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
    LALInferenceAddVariable(model->params, "lambda1", &lambda1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
    LALInferenceAddVariable(model->params, "lambda2", &lambda2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
  }

  /* Only use GR templates */
  /* Fill in the extra parameters for testing GR, if necessary */
  for (UINT4 k=0; k<N_extra_params; k++)
  {
    if(LALInferenceCheckVariable(model->params,list_extra_parameters[k]))
    {
      XLALDictInsert(model->LALpars, list_extra_parameters[k], (void *)LALInferenceGetVariable(model->params,list_extra_parameters[k]), sizeof(double), LAL_D_TYPE_CODE);

      //XLALSimInspiralAddTestGRParam(&nonGRparams,list_extra_parameters[k],*(REAL8 *)LALInferenceGetVariable(model->params,list_extra_parameters[k]));
    }
  }
  /* Fill in PPE params if they are available */
  char PPEparam[64]="";
  const char *PPEnames[]={"aPPE","alphaPPE","bPPE","betaPPE",NULL};
  for(UINT4 idx=0;PPEnames[idx];idx++)
  {
    for(UINT4 ppeidx=0;;ppeidx++)
    {
      sprintf(PPEparam, "%s%d",PPEnames[idx],ppeidx);
      if(LALInferenceCheckVariable(model->params,PPEparam))
	XLALDictInsert(model->LALpars, PPEparam, (void *)LALInferenceGetVariable(model->params,PPEparam), sizeof(double), LAL_D_TYPE_CODE);
      else
        break;
    }
  }

  /* ==== Call the waveform generator ==== */
    /* Correct distance to account for renormalisation of data due to window RMS */
    double corrected_distance = distance * sqrt(model->window->sumofsquares/model->window->data->length);
    XLAL_TRY(ret=XLALSimInspiralChooseFDWaveformSequence (&(model->roq->hptildeLinear), &(model->roq->hctildeLinear), phi0, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI,
                spin1x, spin1y, spin1z, spin2x, spin2y, spin2z, f_ref, corrected_distance, inclination, model->LALpars, approximant, (model->roq->frequencyNodesLinear)), errnum);

    XLAL_TRY(ret=XLALSimInspiralChooseFDWaveformSequence (&(model->roq->hptildeQuadratic), &(model->roq->hctildeQuadratic), phi0, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI,
							spin1x, spin1y, spin1z, spin2x, spin2y, spin2z, f_ref, corrected_distance, inclination, model->LALpars, approximant, (model->roq->frequencyNodesQuadratic)), errnum);

    REAL8 instant = model->freqhPlus->epoch.gpsSeconds + 1e-9*model->freqhPlus->epoch.gpsNanoSeconds;
    LALInferenceSetVariable(model->params, "time", &instant);

        return;
}

void LALInferenceTemplateSineGaussian(LALInferenceModel *model)
/*****************************************************/
/* Sine-Gaussian (burst) template.                   */
/* Signal is (by now?) linearly polarised,           */
/* i.e., the cross-waveform remains zero.            */
/* * * * * * * * * * * * * * * * * * * * * * * * * * */
/* The (plus-) waveform is:                          */
/*   a * exp(-((t-mu)/sigma)^2) * sin(2*pi*f*t-phi)  */
/* Note that by setting f=0, phi=pi/2 you also get   */
/* a `pure' Gaussian template.                       */
/*                                                   */
/* * * * * * * * * * * * * * * * * * * * * * * * * * ************************************/
/* Required (`model->params') parameters are:                                    */
/*   - "time"       (the "mu" parameter of the Gaussian part; REAL8, GPS sec.)          */
/*   - "sigma"      (width, the "sigma" parameter of the Gaussian part; REAL8, seconds) */
/*   - "frequency"  (frequency of the sine part; REAL8, Hertz)                          */
/*   - "phase"      (phase (at above "mu"); REAL8, radians)                             */
/*   - "amplitude"  (amplitude, REAL8)                                                  */
/****************************************************************************************/
{
  double endtime  = *(REAL8*) LALInferenceGetVariable(model->params, "time");       /* time parameter ("mu"), GPS sec.  */
  double sigma = *(REAL8*) LALInferenceGetVariable(model->params, "sigma");      /* width parameter, seconds         */
  double f     = *(REAL8*) LALInferenceGetVariable(model->params, "frequency");  /* frequency, Hz                    */
  double phi   = *(REAL8*) LALInferenceGetVariable(model->params, "phase");      /* phase, rad                       */
  double a     = *(REAL8*) LALInferenceGetVariable(model->params, "amplitude");  /* amplitude                        */
  double t, tsigma, twopif = LAL_TWOPI*f;
  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
  unsigned long i;
  if (sigma <= 0.0) {
    fprintf(stderr, " ERROR in templateSineGaussian(): zero or negative \"sigma\" parameter (sigma=%e).\n", sigma);
    exit(1);
  }
  if (f < 0.0)
    fprintf(stderr, " WARNING in templateSineGaussian(): negative \"frequency\" parameter (f=%e).\n", f);
  if (a < 0.0)
    fprintf(stderr, " WARNING in templateSineGaussian(): negative \"amplitude\" parameter (a=%e).\n", a);
  for (i=0; i<model->timehPlus->data->length; ++i){
    t = ((double)i)*model->deltaT + (epochGPS-endtime);  /* t-mu         */
    tsigma = t/sigma;                                             /* (t-mu)/sigma */
    if (fabs(tsigma) < 5.0)   /*  (only do computations within a 10 sigma range)  */
      model->timehPlus->data->data[i] = a * exp(-0.5*tsigma*tsigma) * sin(twopif*t+phi);
    else
      model->timehPlus->data->data[i] = 0.0;
    model->timehCross->data->data[i] = 0.0;
  }
  model->domain = LAL_SIM_DOMAIN_TIME;
  return;
}

void LALInferenceTemplateDampedSinusoid(LALInferenceModel *model)
/*****************************************************/
/* Damped Sinusoid (burst) template.                 */
/* Signal is linearly polarized,                     */
/* i.e., cross term is zero.                         */
/* * * * * * * * * * * * * * * * * * * * * * * * * * */
/* The (plus-) waveform is an exponentially decaying */
/* sine wave:                                        */
/*   a * exp((t-time)/tau) * sin(2*pi*f*(t-time))    */
/* where "time" is the time parameter denoting the   */
/* instant where the signal starts.                  */
/* * * * * * * * * * * * * * * * * * * * * * * * * * **************************/
/* Required (`model->params') parameters are:                          */
/*   - "time"       (the instant at which the signal starts; REAL8, GPS sec.) */
/*   - "tau"        (width parameter; REAL8, seconds)                         */
/*   - "frequency"  (frequency of the sine part; REAL8, Hertz)                */
/*   - "amplitude"  (amplitude, REAL8)                                        */
/******************************************************************************/
{
  double endtime  = *(REAL8*) LALInferenceGetVariable(model->params, "time");       /* time parameter ("mu"), GPS sec.  */
  double tau   = *(REAL8*) LALInferenceGetVariable(model->params, "tau");        /* width parameter, seconds         */
  double f     = *(REAL8*) LALInferenceGetVariable(model->params, "frequency");  /* frequency, Hz                    */
  double a     = *(REAL8*) LALInferenceGetVariable(model->params, "amplitude");  /* amplitude                        */
  double t, ttau, twopif = LAL_TWOPI*f;
  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
  unsigned long i;
  if (tau <= 0.0) {
    fprintf(stderr, " ERROR in templateDampedSinusoid(): zero or negative \"tau\" parameter (tau=%e).\n", tau);
    exit(1);
  }
  if (f < 0.0)
    fprintf(stderr, " WARNING in templateDampedSinusoid(): negative \"frequency\" parameter (f=%e).\n", f);
  for (i=0; i<model->timehPlus->data->length; ++i){
    t = ((double)i)*model->deltaT + (epochGPS-endtime);  /* t-mu       */
    if ((t>0.0) && ((ttau=t/tau) < 10.0)) /*  (only do computations within a 10 tau range)  */
      model->timehPlus->data->data[i] = a * exp(-ttau) * sin(twopif*t);
    else
      model->timehPlus->data->data[i] = 0.0;
    model->timehCross->data->data[i] = 0.0;
  }
  model->domain = LAL_SIM_DOMAIN_TIME;
  return;
}



void LALInferenceTemplateSinc(LALInferenceModel *model)
/*****************************************************/
/* Sinc function (burst) template.                   */
/* Signal is linearly polarized,                     */
/* i.e., cross term is zero.                         */
/* * * * * * * * * * * * * * * * * * * * * * * * * * */
/* The (plus-) waveform is a sinc function of given  */
/* frequency:                                        */
/*   a * sinc(2*pi*f*(t-time))                       */
/*   = a * sin(2*pi*f*(t-time)) / (2*pi*f*(t-time))  */
/* where "time" is the time parameter denoting the   */
/* signal's central peak location.                   */
/* * * * * * * * * * * * * * * * * * * * * * * * * * *************************/
/* Required (`model->params') parameters are:                         */
/*   - "time"       (the instant at which the signal peaks; REAL8, GPS sec.) */
/*   - "frequency"  (frequency of the sine part; REAL8, Hertz)               */
/*   - "amplitude"  (amplitude, REAL8)                                       */
/*****************************************************************************/
{
  double endtime  = *(REAL8*) LALInferenceGetVariable(model->params, "time");       /* time parameter ("mu"), GPS sec.  */
  double f     = *(REAL8*) LALInferenceGetVariable(model->params, "frequency");  /* frequency, Hz                    */
  double a     = *(REAL8*) LALInferenceGetVariable(model->params, "amplitude");  /* amplitude                        */
  double t, sinArg, sinc, twopif = LAL_TWOPI*f;
  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
  unsigned long i;
  if (f < 0.0)
    fprintf(stderr, " WARNING in templateSinc(): negative \"frequency\" parameter (f=%e).\n", f);
  for (i=0; i<model->timehPlus->data->length; ++i){
    t = ((double)i)*model->deltaT + (epochGPS-endtime);  /* t-mu       */
    sinArg = twopif*t;
    sinc = (sinArg==0.0) ? 1.0 : sin(sinArg)/sinArg;
    model->timehPlus->data->data[i] = a * sinc;
    model->timehCross->data->data[i] = 0.0;
  }
  model->domain = LAL_SIM_DOMAIN_TIME;
  return;
}


void LALInferenceTemplateASinOmegaT(LALInferenceModel *model)
/************************************************************/
/* Trivial h(t)=A*sin(Omega*t) template						*/
/*  Required (`model->params') parameters are:       */
/*   - "A"       (dimensionless amplitude, REAL8)			*/
/*   - "Omega"   (frequency; REAL8, radians/sec)            */
/************************************************************/
{
  double A		= *(REAL8*) LALInferenceGetVariable(model->params, "A");				/* dim-less	   */
  double Omega	= *(REAL8*) LALInferenceGetVariable(model->params, "Omega");			/* rad/sec     */
  double t;
  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));

  unsigned long i;
  for (i=0; i<model->timehPlus->data->length; ++i){
    t = ((double)i)*model->deltaT + (epochGPS);  /* t-mu       */
    model->timehPlus->data->data[i] = A * sin(Omega*t);
    model->timehCross->data->data[i] = 0.0;
  }
  model->domain = LAL_SIM_DOMAIN_TIME;
  return;
}

void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
/*************************************************************************************************************************/
/* Wrapper for LALSimulation waveforms:						                                                             */
/* XLALSimInspiralChooseFDWaveform() and XLALSimInspiralChooseTDWaveform().                                              */
/*                                                                                                                       */
/*  model->params parameters are:										                                                 */
/*  - "name" description; type OPTIONAL (default value)										                             */
/*										                                                                                 */
/*   MODEL PARAMETERS										                                                             */
/*   - "LAL_APPROXIMANT     Approximant;        Approximant                                                              */
/*   - "LAL_PNORDER"        Phase PN order;     INT4                                                                     */
/*   - "LAL_AMPORDER"       Amplitude PN order; INT4 OPTIONAL (-1)                                                       */
/*   - "spinO"              Spin order;         LALSimInspiralSpinOrder OPTIONAL (LAL_SIM_INSPIRAL_SPIN_ORDER_DEFAULT)   */
/*   - "tideO"              Tidal order;        LALSimInspiralTidalOrder OPTIONAL (LAL_SIM_INSPIRAL_TIDAL_ORDER_DEFAULT) */
/*   - "f_ref"               frequency at which the (frequency dependent) parameters are defined; REAL8 OPTIONAL (0.0)    */
/*   - "fLow"               lower frequency bound; REAL8 OPTIONAL (model->fLow)                                          */
/*                                                                                                                       */
/*   MASS PARAMETERS; either:                                                                                            */
/*      - "mass1"           mass of object 1 in solar mass; REAL8								                         */
/*      - "mass2"		        mass of object 1 in solar mass; REAL8								                     */
/*      OR                                                                                                               */
/*      - "chirpmass"       chirpmass in solar mass; REAL8                                                               */
/*      - "q"  asymmetric mass ratio m2/m1, 0<q<1; REAL8                                      */
/*      OR                                                                                                               */
/*      - "chirpmass"       chirpmass in solar mass; REAL8                                                               */
/*      - "eta"             symmetric mass ratio (m1*m2)/(m1+m2)^2; REAL8                                                */
/*                                                                                                                       */
/*   ORIENTATION AND SPIN PARAMETERS                                                                                     */
/*   - "phi0"               reference phase as per LALSimulation convention; REAL8                                       */
/*   - "logdistance"           distance in Mpc                                                                              */
/*   - "costheta_jn");      cos of zenith angle between J and N in radians;            REAL8                                    */
/*   - "phi_jl");        azimuthal angle of L_N on its cone about J radians; REAL8                                    */
/*   - "tilt_spin1");    zenith angle between S1 and LNhat in radians;       REAL8                                    */
/*   - "tilt_spin2");    zenith angle between S2 and LNhat in radians;       REAL8                                    */
/*   - "phi12");         difference in azimuthal angle between S1, S2 in radians;   REAL8                             */
/*   - "a_spin1"            magnitude of spin 1 in general configuration, -1<a_spin1<1; REAL8 OPTIONAL (0.0)              */
/*   - "a_spin2"            magnitude of spin 2 in general configuration, -1<a_spin1<1; REAL8 OPTIONAL (0.0)              */
/*                                                                                                                       */
/*   OTHER PARAMETERS                                                                                                    */
/*   - "lambda1"            tidal parameter of object 1; REAL8  OPTIONAL (0.0)                                           */
/*   - "lambda2"            tidal parameter of object 1; REAL8  OPTIONAL (0.0)                                           */
/*                                                                                                                       */
/*   - "time"               used as an OUTPUT only; REAL8								                                 */
/*                                                                                                                       */
/*                                                                                                                       */
/*   model needs to also contain:                                                                                        */
/*   - model->fLow Unless  - "fLow" OPTIONAL                                                                             */
/*   - model->deltaT                                                                                                     */
/*   - if model->domain == LAL_SIM_DOMAIN_FREQUENCY                                                                      */
/*      - model->deltaF                                                                                                  */
/*      - model->freqhCross                                                                                              */
/*      - model->freqhPlus                                                                                               */
/*   - else                                                                                                              */
/*      - model->timehPlus                                                                                               */
/*      - model->timehCross                                                                                              */
/*************************************************************************************************************************/
{

  Approximant approximant = (Approximant) 0;

  static int sizeWarning = 0;
  int ret=0;
  INT4 errnum=0;

  REAL8TimeSeries *hplus=NULL;  /**< +-polarization waveform [returned] */
  REAL8TimeSeries *hcross=NULL; /**< x-polarization waveform [returned] */
  COMPLEX16FrequencySeries *hptilde=NULL, *hctilde=NULL;
  REAL8 mc;
  REAL8 phi0, deltaT, m1, m2, f_low, f_start, distance, inclination;

  REAL8 *m1_p,*m2_p;
  REAL8 deltaF, f_max;

  /* Sampling rate for time domain models */
  deltaT = model->deltaT;

  if (LALInferenceCheckVariable(model->params, "LAL_APPROXIMANT"))
    approximant = *(Approximant*) LALInferenceGetVariable(model->params, "LAL_APPROXIMANT");
  else {
    XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_APPROXIMANT\" parameter not provided!\n");
    XLAL_ERROR_VOID(XLAL_EDATA);
  }

  if (LALInferenceCheckVariable(model->params, "LAL_PNORDER"))
    XLALSimInspiralWaveformParamsInsertPNPhaseOrder(model->LALpars, *(INT4*) LALInferenceGetVariable(model->params, "LAL_PNORDER"));
  else {
    XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_PNORDER\" parameter not provided!\n");
    XLAL_ERROR_VOID(XLAL_EDATA);
  }

  /* Explicitly set the default amplitude order if one is not specified.
   *   This serves two purposes:
   *     1) The default behavior of the code won't change unexpectedly due to changes in LALSimulation.
   *     2) We need to know the amplitude order in order to set the starting frequency of the waveform properly. */
  if (LALInferenceCheckVariable(model->params, "LAL_AMPORDER"))
    XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,*(INT4*) LALInferenceGetVariable(model->params, "LAL_AMPORDER"));
  else
    if (!XLALSimInspiralWaveformParamsPNAmplitudeOrderIsDefault(model->LALpars))
      XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,-1);

  REAL8 f_ref = 100.0;
  if (LALInferenceCheckVariable(model->params, "f_ref")) f_ref = *(REAL8 *)LALInferenceGetVariable(model->params, "f_ref");

  REAL8 fTemp = f_ref;

  if(LALInferenceCheckVariable(model->params,"chirpmass"))
    {
      mc  = *(REAL8*) LALInferenceGetVariable(model->params, "chirpmass");
      if (LALInferenceCheckVariable(model->params,"q")) {
	REAL8 q = *(REAL8 *)LALInferenceGetVariable(model->params,"q");
	q2masses(mc, q, &m1, &m2);
      } else {
	REAL8 eta = *(REAL8*) LALInferenceGetVariable(model->params, "eta");
	mc2masses(mc, eta, &m1, &m2);
      }
    }
  else if((m1_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass1")) && (m2_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass2")))
    {
      m1=*m1_p;
      m2=*m2_p;
    }
  else
    {
      fprintf(stderr,"No mass parameters found!");
      exit(0);
    }

    if(!LALInferenceCheckVariable(model->params,"logdistance")) distance=LAL_PC_SI * 1e6; /* If distance not given, 1Mpc used */
    else
    {
        distance	= exp(LALInferenceGetREAL8Variable(model->params, "logdistance"))* LAL_PC_SI * 1.0e6;        /* distance (1 Mpc) in units of metres */
    }
  phi0		= LALInferenceGetREAL8Variable(model->params, "phase"); /* START phase as per lalsimulation convention, radians*/

  /* Zenith angle between J and N in radians. Also known as inclination angle when spins are aligned */
  REAL8 thetaJN = acos(LALInferenceGetREAL8Variable(model->params, "costheta_jn"));     /* zenith angle between J and N in radians */

  /* Check if fLow is a model parameter, otherwise use data structure definition */
  if(LALInferenceCheckVariable(model->params, "flow"))
    f_low = *(REAL8*) LALInferenceGetVariable(model->params, "flow");
  else
    f_low = model->fLow;

  f_start = XLALSimInspiralfLow2fStart(f_low, XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(model->LALpars), approximant);

  /* Don't let TaylorF2 generate unphysical inspiral up to Nyquist */
  if (approximant == TaylorF2)
      f_max = 0.0; /* this will stop at ISCO */
  else
      f_max = model->fHigh; /* this will be the highest frequency used across the network */

  /* ==== SPINS ==== */
  /* We will default to spinless signal and then add in the spin components if required */
  /* If there are non-aligned spins, we must convert between the System Frame coordinates
   * and the cartestian coordinates */

  /* The cartesian spin coordinates (default 0), as passed to LALSimulation */
  REAL8 spin1x = 0.0;
  REAL8 spin1y = 0.0;
  REAL8 spin1z = 0.0;
  REAL8 spin2x = 0.0;
  REAL8 spin2y = 0.0;
  REAL8 spin2z = 0.0;

  /* System frame coordinates as used for jump proposals */
  REAL8 a_spin1 = 0.0;  /* Magnitude of spin1 */
  REAL8 a_spin2 = 0.0;  /* Magnitude of spin2 */
  REAL8 phiJL  = 0.0;  /* azimuthal angle of L_N on its cone about J radians */
  REAL8 tilt1   = 0.0;  /* zenith angle between S1 and LNhat in radians */
  REAL8 tilt2   = 0.0;  /* zenith angle between S2 and LNhat in radians */
  REAL8 phi12   = 0.0;  /* difference in azimuthal angle btwn S1, S2 in radians */

  /* Now check if we have spin amplitudes */
  if(LALInferenceCheckVariable(model->params, "a_spin1"))    a_spin1   = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin1");
  if(LALInferenceCheckVariable(model->params, "a_spin2"))    a_spin2   = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin2");

  /* Check if we have spin angles too */
  if(LALInferenceCheckVariable(model->params, "phi_jl"))
      phiJL = LALInferenceGetREAL8Variable(model->params, "phi_jl");
  if(LALInferenceCheckVariable(model->params, "tilt_spin1"))
      tilt1 = LALInferenceGetREAL8Variable(model->params, "tilt_spin1");
  if(LALInferenceCheckVariable(model->params, "tilt_spin2"))
      tilt2 = LALInferenceGetREAL8Variable(model->params, "tilt_spin2");
  if(LALInferenceCheckVariable(model->params, "phi12"))
      phi12 = LALInferenceGetREAL8Variable(model->params, "phi12");

  /* If we have tilt angles zero, then the spins are aligned and we just set the z component */
  /* However, if the waveform supports precession then we still need to get the right coordinate components */
  if(tilt1==0.0 && tilt2==0.0)
  {
      spin1z=a_spin1;
      spin2z=a_spin2;
      inclination = thetaJN; /* Inclination angle is just thetaJN */
  }
  else
  {   /* Template is not aligned-spin only. */
      /* Set all the other spin components according to the angles we received above */
      /* The transformation function doesn't know fLow, so f_ref==0 isn't interpretted as a request to use the starting frequency for reference. */
      if(fTemp==0.0)
        fTemp = f_start;
      // If we have a *time-domain* AND *precessing*  EOB template, set fTemp to f_start.
      // This is done since EOB only uses f_start and ignores f_ref.
      if(model->domain == LAL_SIM_DOMAIN_TIME && (strstr(XLALSimInspiralGetStringFromApproximant(approximant),"EOB") != NULL)){
	fTemp = f_start;
      }
      XLAL_TRY(ret=XLALSimInspiralTransformPrecessingNewInitialConditions(
                    &inclination, &spin1x, &spin1y, &spin1z, &spin2x, &spin2y, &spin2z,
                    thetaJN, phiJL, tilt1, tilt2, phi12, a_spin1, a_spin2, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, fTemp, phi0), errnum);
      if (ret == XLAL_FAILURE)
      {
        XLALPrintError(" ERROR in XLALSimInspiralTransformPrecessingNewInitialConditions(): error converting angles. errnum=%d: %s\n",errnum, XLALErrorString(errnum) );
        return;
      }
  }


  
/* ==== Spin induced quadrupole moment PARAMETERS ==== */

 if(LALInferenceCheckVariable(model->params, "dQuadMonS")&&LALInferenceCheckVariable(model->params, "dQuadMonA")){
    REAL8 dQuadMon1=0.;
    REAL8 dQuadMon2=0.;
    REAL8 dQuadMonS = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonS");
    REAL8 dQuadMonA = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonA");
    LALInferencedQuadMonSdQuadMonA(dQuadMonS,dQuadMonA,&dQuadMon1,&dQuadMon2);
    XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
    XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
        fprintf(stdout,"Both dQuadMonS and dQaudMonA are  sampled");
    fprintf(stdout,"dQM1: %e, dQM2: %e, dQMS: %e, dQMA: %e \n",dQuadMon1,dQuadMon2,dQuadMonS,dQuadMonA);
  }
 else  if(LALInferenceCheckVariable(model->params, "dQuadMon1")&&LALInferenceCheckVariable(model->params, "dQuadMon2")){
    REAL8 dQuadMon1 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon1");
    REAL8 dQuadMon2 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon2");
    XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
    XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
        fprintf(stdout,"Both dQuadMon1 and dQaudMon2 are  sampled");
    fprintf(stdout,"dQM1: %e, dQM2: %e \n",dQuadMon1,dQuadMon2);
  }


  



/* ==== TIDAL PARAMETERS ==== */
  if(LALInferenceCheckVariable(model->params, "lambda1"))
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda1"));
  if(LALInferenceCheckVariable(model->params, "lambda2"))
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda2"));
  REAL8 lambdaT = 0.;
  REAL8 dLambdaT = 0.;
  REAL8 sym_mass_ratio_eta = 0.;
  if(LALInferenceCheckVariable(model->params, "lambdaT")&&LALInferenceCheckVariable(model->params, "dLambdaT")){
    REAL8 lambda1=0.;
    REAL8 lambda2=0.;
    lambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "lambdaT");
    dLambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "dLambdaT");
    sym_mass_ratio_eta = m1*m2/((m1+m2)*(m1+m2));
    LALInferenceLambdaTsEta2Lambdas(lambdaT,dLambdaT,sym_mass_ratio_eta,&lambda1,&lambda2);
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
  }
  if(model->eos_fam)
  {
      LALSimNeutronStarFamily *eos_fam = model->eos_fam;
      REAL8 r1=0, r2=0, k2_1=0, k2_2=0, lambda1=0, lambda2=0;
      REAL8 mass_max = XLALSimNeutronStarMaximumMass(eos_fam) / LAL_MSUN_SI;
      REAL8 mass_min = XLALSimNeutronStarFamMinimumMass(eos_fam) / LAL_MSUN_SI;

      if(m1<mass_max && m1>mass_min)
      {
        /* Compute l1, l2 from mass and EOS */
        r1 = XLALSimNeutronStarRadius(m1*LAL_MSUN_SI, eos_fam);
        k2_1 = XLALSimNeutronStarLoveNumberK2(m1*LAL_MSUN_SI, eos_fam);
        lambda1 = (2./3.)*k2_1 * pow(r1/(m1*LAL_MRSUN_SI), 5.0);
      }
      if(m2<mass_max && m2>mass_min)
      {
         r2 = XLALSimNeutronStarRadius(m2*LAL_MSUN_SI, eos_fam);
         k2_2 = XLALSimNeutronStarLoveNumberK2(m2*LAL_MSUN_SI, eos_fam);
         lambda2 = (2./3.)*k2_2 * pow(r2/(m2*LAL_MRSUN_SI), 5.0);          
      }
      /* Set waveform params */
      XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
      XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
      REAL8 dQuadMon1 = dquadmon_from_lambda(lambda1);
      REAL8 dQuadMon2 = dquadmon_from_lambda(lambda2);
      XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars, dQuadMon1);
      XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars, dQuadMon2);

      /* Calculate maximum frequency */
      /* If both lambdas are non-zero compute EOS-dependent f_max */
      if((lambda1 > 0) && (lambda2 > 0) && (approximant == TaylorF2))
      {
        /* Start with ISCO */
        f_max = 1. / (pow(6,1.5) * LAL_PI * (m1*LAL_MTSUN_SI + m2*LAL_MTSUN_SI));
        REAL8 log_lambda1 = log(lambda1);
        REAL8 log_lambda2 = log(lambda2);
        REAL8 compactness1 = 0.371 - 0.0391 * log_lambda1 + 0.001056 * log_lambda1 * log_lambda1;
        REAL8 compactness2 = 0.371 - 0.0391 * log_lambda2 + 0.001056 * log_lambda2 * log_lambda2;
        REAL8 rad1 = m1*LAL_MTSUN_SI / compactness1;
        REAL8 rad2 = m2*LAL_MTSUN_SI / compactness2;
        /* Use the smaller of ISCO and the EOS-dependent termination */
        REAL8 fmax_eos = 1. / LAL_PI * pow((m1*LAL_MTSUN_SI + m2*LAL_MTSUN_SI) / pow((rad1 + rad2),3.0),0.5);
        if (fmax_eos < f_max)
        {
          f_max = fmax_eos;
        }
      }