LALInferenceTemplate.c 115 KB
Newer Older
1
/*
Will Farr's avatar
Will Farr committed
2
3
 *  LALInferenceTemplate.c: Bayesian Followup, template calls to LAL
 *  template functions. Temporary GeneratePPN
Christian Roever's avatar
Christian Roever committed
4
 *
Will Farr's avatar
Will Farr committed
5
 *  Copyright (C) 2009 Ilya Mandel, Vivien Raymond, Christian Roever,
6
 *  Marc van der Sluys, John Veitch, Will M. Farr and Salvatore Vitale
Christian Roever's avatar
Christian Roever committed
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 *
 *
 *  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
21
22
 *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 *  MA  02110-1301  USA
Christian Roever's avatar
Christian Roever committed
23
 */
24
25
26

#include <stdio.h>
#include <stdlib.h>
27
#include <lal/LALInspiral.h>
Christian Roever's avatar
Christian Roever committed
28
#include <lal/SeqFactories.h>
29
#include <lal/TimeSeries.h>
30
#include <lal/FrequencySeries.h>
31
#include <lal/Date.h>
32
33
#include <lal/VectorOps.h>
#include <lal/TimeFreqFFT.h>
Vivien Raymond's avatar
Vivien Raymond committed
34
#include <lal/GenerateInspiral.h>
35
#include <lal/GenerateInspRing.h>
John Douglas Veitch's avatar
John Douglas Veitch committed
36
#include <lal/LALStatusMacros.h>
Matthew Pitkin's avatar
Matthew Pitkin committed
37
#include <lal/LALInference.h>
38
#include <lal/XLALError.h>
39
#include <lal/LIGOMetadataRingdownUtils.h>
40
#include <lal/LALSimInspiral.h>
41
#include <lal/LALInferenceTemplate.h>
42
#include <lal/LALInferenceMultibanding.h>
John Douglas Veitch's avatar
John Douglas Veitch committed
43
#include <lal/LALSimNeutronStar.h>
44

45
/* LIB imports*/
46
#include <lal/LALInferenceBurstRoutines.h>
47
48


John Douglas Veitch's avatar
John Douglas Veitch committed
49
50
51
52
53
54
#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$"
55

56
57
58
59
60
61
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

62
63
64
65
/* 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

66
67
68
69
#define Pi_p2 9.8696044010893586188344909998761511
#define Pi_p2by3 2.1450293971110256000774441009412356
#define log4 1.3862943611198906188344642429163531

70
static void q2masses(double mc, double q, double *m1, double *m2);
John Douglas Veitch's avatar
John Douglas Veitch committed
71

72
/* list of testing GR parameters to be passed to the waveform */
Anuradha Samajdar's avatar
Anuradha Samajdar committed
73
/* 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 */
74
75


76
77
78
79
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;
80

NV KRISHNENDU's avatar
NV KRISHNENDU committed
81

Ian Harry's avatar
Ian Harry committed
82
83
84
85
86
87
/* 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.
 */
88
89
90
91
92
static REAL8 dquadmon_from_lambda(REAL8 lambdav);
static REAL8 dquadmon_from_lambda(REAL8 lambdav)
{
    
    double ll = log(lambdav);
John Douglas Veitch's avatar
John Douglas Veitch committed
93
94
95
    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);
96
}
97
98
99
100
101
102
103
104

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);
105

106
107
  /* Loop over reduced frequency set */
  for(UINT4 i=0;i<freqs->length-1;i++)
108
  {
109
110
111
112
113
114
115
116
117
    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 */
118

119
120
121
122
123
124
125
126
127
    /* 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;}
128

129
130
131
132
    double dpsidf = dpsi/df;
    double dampdf = (endamp - startamp)/df;

    double damp = dampdf *deltaF;
133

134
135
136
137
138
139
140
141
    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;
142

143
        f<endf;
144

145
146
147
148
149
150
151
152
153
154
155
156
157
158
        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;
}


159
void LALInferenceTemplateNullFreqdomain(LALInferenceModel *model)
160
161
162
163
164
/**********************************************/
/* returns a frequency-domain 'null' template */
/* (all zeroes, implying no signal present).  */
/**********************************************/
{
John Douglas Veitch's avatar
John Douglas Veitch committed
165
  UINT4 i;
166
  if ((model->freqhPlus==NULL) || (model->freqhCross==NULL)) {
167
    XLALPrintError(" ERROR in templateNullFreqdomain(): encountered unallocated 'freqModelhPlus/-Cross'.\n");
168
    XLAL_ERROR_VOID(XLAL_EFAULT);
169
  }
170
171
172
  for (i=0; i<model->freqhPlus->data->length; ++i){
    model->freqhPlus->data->data[i] = 0.0;
    model->freqhCross->data->data[i] = 0.0;
173
  }
174
  model->domain = LAL_SIM_DOMAIN_FREQUENCY;
175
176
177
178
179
  return;
}



180
void LALInferenceTemplateNullTimedomain(LALInferenceModel *model)
181
182
183
184
185
/*********************************************/
/* returns a time-domain 'null' template     */
/* (all zeroes, implying no signal present). */
/*********************************************/
{
John Douglas Veitch's avatar
John Douglas Veitch committed
186
  UINT4 i;
187
  if ((model->timehPlus==NULL) || (model->timehCross==NULL)) {
188
    XLALPrintError(" ERROR in templateNullTimedomain(): encountered unallocated 'timeModelhPlus/-Cross'.\n");
189
    XLAL_ERROR_VOID(XLAL_EFAULT);
190
  }
191
192
193
  for (i=0; i<model->timehPlus->data->length; ++i){
    model->timehPlus->data->data[i]  = 0.0;
    model->timehCross->data->data[i] = 0.0;
194
  }
195
  model->domain = LAL_SIM_DOMAIN_TIME;
196
197
198
199
200
201
202
  return;
}



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

203

Rory Smith's avatar
Rory Smith committed
204
static void mc2masses(double mc, double eta, double *m1, double *m2);
205

Rory Smith's avatar
Rory Smith committed
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
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;
}
227

228
229
230
void LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence(LALInferenceModel *model){
/*************************************************************************************************************************/
  Approximant approximant = (Approximant) 0;
231

232
233
  int ret=0;
  INT4 errnum=0;
234

235
236
237
238
  model->roq->hptildeLinear=NULL, model->roq->hctildeLinear=NULL;
  model->roq->hptildeQuadratic=NULL, model->roq->hctildeQuadratic=NULL;
  REAL8 mc;
  REAL8 phi0, m1, m2, distance, inclination;
239

240
  REAL8 *m1_p,*m2_p;
241

242
243
244
245
246
  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);
247
248
  }

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

256
257
258
259
260
  /* 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"))
261
    XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars, *(INT4 *) LALInferenceGetVariable(model->params, "LAL_AMPORDER"));
262
  else
263
264
265
    if (!XLALSimInspiralWaveformParamsPNAmplitudeOrderIsDefault(model->LALpars))
      XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,-1);

266
267
268
269
270
  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"))
271
    {
272
273
274
275
276
277
278
279
280
281
      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")))
282
    {
283
284
      m1=*m1_p;
      m2=*m2_p;
285
286
287
    }
  else
    {
288
289
290
      fprintf(stderr,"No mass parameters found!");
      exit(0);
    }
291
292
293
294
295
296
297
298
  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 */
  }
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340

  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 */
341
  if(tilt1==0.0 && tilt2==0.0)
342
343
344
345
346
347
348
349
350
351
352
  {
      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,
353
                    thetaJN, phiJL, tilt1, tilt2, phi12, a_spin1, a_spin2, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, fTemp, phi0), errnum);
354
355
356
357
      if (ret == XLAL_FAILURE)
      {
        XLALPrintError(" ERROR in XLALSimInspiralTransformPrecessingNewInitialConditions(): error converting angles. errnum=%d\n",errnum );
        return;
358
      }
359
  }
360

361

362
363
364
365
366
367
368
369
370
371
372
373
374
/* ==== 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);
  }
NV KRISHNENDU's avatar
NV KRISHNENDU committed
375
 else  if(LALInferenceCheckVariable(model->params, "dQuadMon1")&&LALInferenceCheckVariable(model->params, "dQuadMon2")){
376
377
    REAL8 dQuadMon1 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon1");
    REAL8 dQuadMon2 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon2");
378
379
    XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
    XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
NV KRISHNENDU's avatar
NV KRISHNENDU committed
380
381
        fprintf(stdout,"Both dQuadMon1 and dQaudMon2 are  sampled");
    fprintf(stdout,"dQM1: %e, dQM2: %e \n",dQuadMon1,dQuadMon2);
382
383
384
  }


NV KRISHNENDU's avatar
NV KRISHNENDU committed
385

386
  /* ==== TIDAL PARAMETERS ==== */
387
388
389
390
  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"));
391
392
393
  REAL8 lambdaT = 0.;
  REAL8 dLambdaT = 0.;
  REAL8 sym_mass_ratio_eta = 0.;
394

395
  if(LALInferenceCheckVariable(model->params, "lambdaT")&&LALInferenceCheckVariable(model->params, "dLambdaT")){
396
397
    REAL8 lambda1=0.;
    REAL8 lambda2=0.;
398
399
400
401
    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);
402
403
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars,lambda1);
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars,lambda2);
404
405
  }

406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
  /* ==== 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);
  }
423

424
425
426
427
428
429
430
431
432
433
434
435
436
437
  /* ==== 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};
438
    LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
439
440
441
442
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
  }

Carl-Johan Haster's avatar
Carl-Johan Haster committed
443
444
445
446
447
448
449
450
451
452
453
  /* ==== 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);
  }

454
  /* Only use GR templates */
455
456
457
458
459
  /* 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]))
    {
460
461
462
      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]));
463
464
465
466
467
468
469
470
471
472
473
    }
  }
  /* 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))
474
	XLALDictInsert(model->LALpars, PPEparam, (void *)LALInferenceGetVariable(model->params,PPEparam), sizeof(double), LAL_D_TYPE_CODE);
475
476
477
478
      else
        break;
    }
  }
479
480

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

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

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

        return;
493
494
}

495
void LALInferenceTemplateSineGaussian(LALInferenceModel *model)
496
497
498
499
500
501
502
/*****************************************************/
/* 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)  */
503
504
/* Note that by setting f=0, phi=pi/2 you also get   */
/* a `pure' Gaussian template.                       */
505
506
/*                                                   */
/* * * * * * * * * * * * * * * * * * * * * * * * * * ************************************/
507
/* Required (`model->params') parameters are:                                    */
508
509
510
511
512
513
514
/*   - "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)                                                  */
/****************************************************************************************/
{
515
516
517
518
519
  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                        */
520
  double t, tsigma, twopif = LAL_TWOPI*f;
521
  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
Matthew Pitkin's avatar
Matthew Pitkin committed
522
  unsigned long i;
523
524
525
526
527
528
529
530
  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);
531
532
  for (i=0; i<model->timehPlus->data->length; ++i){
    t = ((double)i)*model->deltaT + (epochGPS-endtime);  /* t-mu         */
533
534
    tsigma = t/sigma;                                             /* (t-mu)/sigma */
    if (fabs(tsigma) < 5.0)   /*  (only do computations within a 10 sigma range)  */
535
      model->timehPlus->data->data[i] = a * exp(-0.5*tsigma*tsigma) * sin(twopif*t+phi);
536
    else
537
538
      model->timehPlus->data->data[i] = 0.0;
    model->timehCross->data->data[i] = 0.0;
539
  }
540
  model->domain = LAL_SIM_DOMAIN_TIME;
541
  return;
542
}
543

544
void LALInferenceTemplateDampedSinusoid(LALInferenceModel *model)
545
546
547
548
549
550
551
552
553
554
555
/*****************************************************/
/* 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.                  */
/* * * * * * * * * * * * * * * * * * * * * * * * * * **************************/
556
/* Required (`model->params') parameters are:                          */
557
558
559
560
561
562
/*   - "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)                                        */
/******************************************************************************/
{
563
564
565
566
  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                        */
567
  double t, ttau, twopif = LAL_TWOPI*f;
568
  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
Matthew Pitkin's avatar
Matthew Pitkin committed
569
  unsigned long i;
570
571
572
573
574
575
  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);
576
577
  for (i=0; i<model->timehPlus->data->length; ++i){
    t = ((double)i)*model->deltaT + (epochGPS-endtime);  /* t-mu       */
578
    if ((t>0.0) && ((ttau=t/tau) < 10.0)) /*  (only do computations within a 10 tau range)  */
579
      model->timehPlus->data->data[i] = a * exp(-ttau) * sin(twopif*t);
580
    else
581
582
      model->timehPlus->data->data[i] = 0.0;
    model->timehCross->data->data[i] = 0.0;
583
  }
584
  model->domain = LAL_SIM_DOMAIN_TIME;
585
586
587
588
  return;
}


Christian Roever's avatar
Christian Roever committed
589

590
void LALInferenceTemplateSinc(LALInferenceModel *model)
591
592
593
594
595
596
597
/*****************************************************/
/* Sinc function (burst) template.                   */
/* Signal is linearly polarized,                     */
/* i.e., cross term is zero.                         */
/* * * * * * * * * * * * * * * * * * * * * * * * * * */
/* The (plus-) waveform is a sinc function of given  */
/* frequency:                                        */
Christian Roever's avatar
Christian Roever committed
598
599
/*   a * sinc(2*pi*f*(t-time))                       */
/*   = a * sin(2*pi*f*(t-time)) / (2*pi*f*(t-time))  */
600
601
602
/* where "time" is the time parameter denoting the   */
/* signal's central peak location.                   */
/* * * * * * * * * * * * * * * * * * * * * * * * * * *************************/
603
/* Required (`model->params') parameters are:                         */
604
605
606
607
608
/*   - "time"       (the instant at which the signal peaks; REAL8, GPS sec.) */
/*   - "frequency"  (frequency of the sine part; REAL8, Hertz)               */
/*   - "amplitude"  (amplitude, REAL8)                                       */
/*****************************************************************************/
{
609
610
611
  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                        */
612
  double t, sinArg, sinc, twopif = LAL_TWOPI*f;
613
  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
Matthew Pitkin's avatar
Matthew Pitkin committed
614
  unsigned long i;
615
616
  if (f < 0.0)
    fprintf(stderr, " WARNING in templateSinc(): negative \"frequency\" parameter (f=%e).\n", f);
617
618
  for (i=0; i<model->timehPlus->data->length; ++i){
    t = ((double)i)*model->deltaT + (epochGPS-endtime);  /* t-mu       */
619
    sinArg = twopif*t;
620
    sinc = (sinArg==0.0) ? 1.0 : sin(sinArg)/sinArg;
621
622
    model->timehPlus->data->data[i] = a * sinc;
    model->timehCross->data->data[i] = 0.0;
623
  }
624
  model->domain = LAL_SIM_DOMAIN_TIME;
625
626
  return;
}
627
628


629
void LALInferenceTemplateASinOmegaT(LALInferenceModel *model)
630
631
/************************************************************/
/* Trivial h(t)=A*sin(Omega*t) template						*/
632
/*  Required (`model->params') parameters are:       */
633
634
635
636
/*   - "A"       (dimensionless amplitude, REAL8)			*/
/*   - "Omega"   (frequency; REAL8, radians/sec)            */
/************************************************************/
{
637
638
  double A		= *(REAL8*) LALInferenceGetVariable(model->params, "A");				/* dim-less	   */
  double Omega	= *(REAL8*) LALInferenceGetVariable(model->params, "Omega");			/* rad/sec     */
639
  double t;
640
  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
641

Matthew Pitkin's avatar
Matthew Pitkin committed
642
  unsigned long i;
643
  for (i=0; i<model->timehPlus->data->length; ++i){
644
    t = ((double)i)*model->deltaT + (epochGPS);  /* t-mu       */
645
646
    model->timehPlus->data->data[i] = A * sin(Omega*t);
    model->timehCross->data->data[i] = 0.0;
647
  }
648
  model->domain = LAL_SIM_DOMAIN_TIME;
649
650
  return;
}
Vivien Raymond's avatar
Vivien Raymond committed
651

652
void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
653
/*************************************************************************************************************************/
654
/* Wrapper for LALSimulation waveforms:						                                                             */
655
656
/* XLALSimInspiralChooseFDWaveform() and XLALSimInspiralChooseTDWaveform().                                              */
/*                                                                                                                       */
657
/*  model->params parameters are:										                                                 */
658
659
660
/*  - "name" description; type OPTIONAL (default value)										                             */
/*										                                                                                 */
/*   MODEL PARAMETERS										                                                             */
661
/*   - "LAL_APPROXIMANT     Approximant;        Approximant                                                              */
662
663
664
665
/*   - "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) */
666
/*   - "f_ref"               frequency at which the (frequency dependent) parameters are defined; REAL8 OPTIONAL (0.0)    */
667
/*   - "fLow"               lower frequency bound; REAL8 OPTIONAL (model->fLow)                                          */
668
669
/*                                                                                                                       */
/*   MASS PARAMETERS; either:                                                                                            */
670
671
/*      - "mass1"           mass of object 1 in solar mass; REAL8								                         */
/*      - "mass2"		        mass of object 1 in solar mass; REAL8								                     */
672
673
/*      OR                                                                                                               */
/*      - "chirpmass"       chirpmass in solar mass; REAL8                                                               */
674
/*      - "q"  asymmetric mass ratio m2/m1, 0<q<1; REAL8                                      */
675
676
/*      OR                                                                                                               */
/*      - "chirpmass"       chirpmass in solar mass; REAL8                                                               */
677
/*      - "eta"             symmetric mass ratio (m1*m2)/(m1+m2)^2; REAL8                                                */
678
679
680
/*                                                                                                                       */
/*   ORIENTATION AND SPIN PARAMETERS                                                                                     */
/*   - "phi0"               reference phase as per LALSimulation convention; REAL8                                       */
Ben Farr's avatar
Ben Farr committed
681
/*   - "logdistance"           distance in Mpc                                                                              */
682
/*   - "costheta_jn");      cos of zenith angle between J and N in radians;            REAL8                                    */
683
684
685
686
/*   - "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                             */
687
688
/*   - "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)              */
689
690
691
692
693
/*                                                                                                                       */
/*   OTHER PARAMETERS                                                                                                    */
/*   - "lambda1"            tidal parameter of object 1; REAL8  OPTIONAL (0.0)                                           */
/*   - "lambda2"            tidal parameter of object 1; REAL8  OPTIONAL (0.0)                                           */
/*                                                                                                                       */
694
/*   - "time"               used as an OUTPUT only; REAL8								                                 */
695
696
/*                                                                                                                       */
/*                                                                                                                       */
697
698
699
700
701
702
703
/*   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                                                                                               */
704
/*   - else                                                                                                              */
705
706
/*      - model->timehPlus                                                                                               */
/*      - model->timehCross                                                                                              */
707
/*************************************************************************************************************************/
708
{
709

710
  Approximant approximant = (Approximant) 0;
711

Will Farr's avatar
Will Farr committed
712
  static int sizeWarning = 0;
713
  int ret=0;
714
  INT4 errnum=0;
715

716
717
  REAL8TimeSeries *hplus=NULL;  /**< +-polarization waveform [returned] */
  REAL8TimeSeries *hcross=NULL; /**< x-polarization waveform [returned] */
718
  COMPLEX16FrequencySeries *hptilde=NULL, *hctilde=NULL;
Will Farr's avatar
Will Farr committed
719
  REAL8 mc;
720
  REAL8 phi0, deltaT, m1, m2, f_low, f_start, distance, inclination;
721

722
  REAL8 *m1_p,*m2_p;
Will Farr's avatar
Will Farr committed
723
  REAL8 deltaF, f_max;
724

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

728
729
  if (LALInferenceCheckVariable(model->params, "LAL_APPROXIMANT"))
    approximant = *(Approximant*) LALInferenceGetVariable(model->params, "LAL_APPROXIMANT");
Will Farr's avatar
Will Farr committed
730
731
732
733
  else {
    XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_APPROXIMANT\" parameter not provided!\n");
    XLAL_ERROR_VOID(XLAL_EDATA);
  }
734

735
  if (LALInferenceCheckVariable(model->params, "LAL_PNORDER"))
736
    XLALSimInspiralWaveformParamsInsertPNPhaseOrder(model->LALpars, *(INT4*) LALInferenceGetVariable(model->params, "LAL_PNORDER"));
Will Farr's avatar
Will Farr committed
737
738
739
740
  else {
    XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_PNORDER\" parameter not provided!\n");
    XLAL_ERROR_VOID(XLAL_EDATA);
  }
741
742
743
744
745

  /* 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. */
746
  if (LALInferenceCheckVariable(model->params, "LAL_AMPORDER"))
747
    XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,*(INT4*) LALInferenceGetVariable(model->params, "LAL_AMPORDER"));
748
  else
749
750
    if (!XLALSimInspiralWaveformParamsPNAmplitudeOrderIsDefault(model->LALpars))
      XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,-1);
751

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

755
  REAL8 fTemp = f_ref;
756

757
  if(LALInferenceCheckVariable(model->params,"chirpmass"))
Will Farr's avatar
Will Farr committed
758
    {
759
      mc  = *(REAL8*) LALInferenceGetVariable(model->params, "chirpmass");
760
761
      if (LALInferenceCheckVariable(model->params,"q")) {
	REAL8 q = *(REAL8 *)LALInferenceGetVariable(model->params,"q");
Will Farr's avatar
Will Farr committed
762
763
	q2masses(mc, q, &m1, &m2);
      } else {
764
	REAL8 eta = *(REAL8*) LALInferenceGetVariable(model->params, "eta");
Will Farr's avatar
Will Farr committed
765
766
767
	mc2masses(mc, eta, &m1, &m2);
      }
    }
768
  else if((m1_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass1")) && (m2_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass2")))
Will Farr's avatar
Will Farr committed
769
770
771
772
773
774
775
776
777
    {
      m1=*m1_p;
      m2=*m2_p;
    }
  else
    {
      fprintf(stderr,"No mass parameters found!");
      exit(0);
    }
778

John Douglas Veitch's avatar
John Douglas Veitch committed
779
780
781
782
783
    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 */
    }
784
  phi0		= LALInferenceGetREAL8Variable(model->params, "phase"); /* START phase as per lalsimulation convention, radians*/
785

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

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

795
  f_start = XLALSimInspiralfLow2fStart(f_low, XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(model->LALpars), approximant);
796
797
798
799
800
801

  /* 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 */
802

803
804
805
806
807
808
809
810
811
812
813
814
  /* ==== 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;
815

816
817
818
  /* System frame coordinates as used for jump proposals */
  REAL8 a_spin1 = 0.0;  /* Magnitude of spin1 */
  REAL8 a_spin2 = 0.0;  /* Magnitude of spin2 */
819
  REAL8 phiJL  = 0.0;  /* azimuthal angle of L_N on its cone about J radians */
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
  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 */
840
  if(tilt1==0.0 && tilt2==0.0)
841
842
843
844
  {
      spin1z=a_spin1;
      spin2z=a_spin2;
      inclination = thetaJN; /* Inclination angle is just thetaJN */
845
  }
846
847
848
849
850
851
  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;
852
853
854
855
856
      // 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;
      }
Ben Farr's avatar
Ben Farr committed
857
      XLAL_TRY(ret=XLALSimInspiralTransformPrecessingNewInitialConditions(
858
                    &inclination, &spin1x, &spin1y, &spin1z, &spin2x, &spin2y, &spin2z,
859
                    thetaJN, phiJL, tilt1, tilt2, phi12, a_spin1, a_spin2, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, fTemp, phi0), errnum);
860
861
      if (ret == XLAL_FAILURE)
      {
862
        XLALPrintError(" ERROR in XLALSimInspiralTransformPrecessingNewInitialConditions(): error converting angles. errnum=%d: %s\n",errnum, XLALErrorString(errnum) );
863
864
        return;
      }
865
  }
866

867

NV KRISHNENDU's avatar
NV KRISHNENDU committed
868
  
869
/* ==== Spin induced quadrupole moment PARAMETERS ==== */
870

871
 if(LALInferenceCheckVariable(model->params, "dQuadMonS")&&LALInferenceCheckVariable(model->params, "dQuadMonA")){
872
873
    REAL8 dQuadMon1=0.;
    REAL8 dQuadMon2=0.;
874
875
    REAL8 dQuadMonS = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonS");
    REAL8 dQuadMonA = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonA");
876
877
878
    LALInferencedQuadMonSdQuadMonA(dQuadMonS,dQuadMonA,&dQuadMon1,&dQuadMon2);
    XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
    XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
879
880
        fprintf(stdout,"Both dQuadMonS and dQaudMonA are  sampled");
    fprintf(stdout,"dQM1: %e, dQM2: %e, dQMS: %e, dQMA: %e \n",dQuadMon1,dQuadMon2,dQuadMonS,dQuadMonA);
881
  }
NV KRISHNENDU's avatar
NV KRISHNENDU committed
882
 else  if(LALInferenceCheckVariable(model->params, "dQuadMon1")&&LALInferenceCheckVariable(model->params, "dQuadMon2")){
883
884
885
886
    REAL8 dQuadMon1 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon1");
    REAL8 dQuadMon2 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon2");
    XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
    XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
NV KRISHNENDU's avatar
NV KRISHNENDU committed
887
888
        fprintf(stdout,"Both dQuadMon1 and dQaudMon2 are  sampled");
    fprintf(stdout,"dQM1: %e, dQM2: %e \n",dQuadMon1,dQuadMon2);
889
  }
890
891


NV KRISHNENDU's avatar
NV KRISHNENDU committed
892
893
894
895
896
  



/* ==== TIDAL PARAMETERS ==== */
897
898
899
900
  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"));
901
902
903
  REAL8 lambdaT = 0.;
  REAL8 dLambdaT = 0.;
  REAL8 sym_mass_ratio_eta = 0.;
904
  if(LALInferenceCheckVariable(model->params, "lambdaT")&&LALInferenceCheckVariable(model->params, "dLambdaT")){
905
906
    REAL8 lambda1=0.;
    REAL8 lambda2=0.;
907
908
    lambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "lambdaT");
    dLambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "dLambdaT");
909
910
    sym_mass_ratio_eta = m1*m2/((m1+m2)*(m1+m2));
    LALInferenceLambdaTsEta2Lambdas(lambdaT,dLambdaT,sym_mass_ratio_eta,&lambda1,&lambda2);
911
912
    XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
    XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
913
  }
914
  if(model->eos_fam)
John Douglas Veitch's avatar
John Douglas Veitch committed
915
  {
916
      LALSimNeutronStarFamily *eos_fam = model->eos_fam;
John Douglas Veitch's avatar
John Douglas Veitch committed
917
      REAL8 r1=0, r2=0, k2_1=0, k2_2=0, lambda1=0, lambda2=0;
John Douglas Veitch's avatar
John Douglas Veitch committed
918
919
      REAL8 mass_max = XLALSimNeutronStarMaximumMass(eos_fam) / LAL_MSUN_SI;
      REAL8 mass_min = XLALSimNeutronStarFamMinimumMass(eos_fam) / LAL_MSUN_SI;
920

John Douglas Veitch's avatar
John Douglas Veitch committed
921
922
923
      if(m1<mass_max && m1>mass_min)
      {
        /* Compute l1, l2 from mass and EOS */
John Douglas Veitch's avatar
John Douglas Veitch committed
924
925
926
        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);
John Douglas Veitch's avatar
John Douglas Veitch committed
927
928
929
      }
      if(m2<mass_max && m2>mass_min)
      {
John Douglas Veitch's avatar
John Douglas Veitch committed
930
931
932
         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);          
John Douglas Veitch's avatar
John Douglas Veitch committed
933
      }
John Douglas Veitch's avatar
John Douglas Veitch committed
934
935
936
      /* Set waveform params */
      XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
      XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
937
      REAL8 dQuadMon1 = dquadmon_from_lambda(lambda1);
John Douglas Veitch's avatar
John Douglas Veitch committed
938
      REAL8 dQuadMon2 = dquadmon_from_lambda(lambda2);
939
940
      XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars, dQuadMon1);
      XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars, dQuadMon2);
Ian Harry's avatar
Ian Harry committed
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960

      /* 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;
        }
      }
John Douglas Veitch's avatar
John Douglas Veitch committed