LALInference  4.1.6.1-89842e6
LALInferenceTemplate.c
Go to the documentation of this file.
1 /*
2  * LALInferenceTemplate.c: Bayesian Followup, template calls to LAL
3  * template functions. Temporary GeneratePPN
4  *
5  * Copyright (C) 2009 Ilya Mandel, Vivien Raymond, Christian Roever,
6  * Marc van der Sluys, John Veitch, Will M. Farr and Salvatore Vitale
7  *
8  *
9  * This program is free software; you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation; either version 2 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with with program; see the file COPYING. If not, write to the
21  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
22  * MA 02110-1301 USA
23  */
24 
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <lal/LALInspiral.h>
28 #include <lal/SeqFactories.h>
29 #include <lal/TimeSeries.h>
30 #include <lal/FrequencySeries.h>
31 #include <lal/Date.h>
32 #include <lal/VectorOps.h>
33 #include <lal/TimeFreqFFT.h>
34 #include <lal/GenerateInspiral.h>
35 #include <lal/GenerateInspRing.h>
36 #include <lal/LALStatusMacros.h>
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>
43 #include <lal/LALSimNeutronStar.h>
44 #include <lal/LALSimInspiralTestingGRCorrections.h>
45 
46 #include <lal/LALSimSphHarmMode.h>
47 
48 /* LIB imports*/
49 #include <lal/LALInferenceBurstRoutines.h>
50 
51 
52 #define PROGRAM_NAME "LALInferenceTemplate.c"
53 #define CVS_ID_STRING "$Id$"
54 #define CVS_REVISION "$Revision$"
55 #define CVS_SOURCE "$Source$"
56 #define CVS_DATE "$Date$"
57 #define CVS_NAME_STRING "$Name$"
58 
59 #ifdef __GNUC__
60 #define UNUSED __attribute__ ((unused))
61 #else
62 #define UNUSED
63 #endif
64 
65 /* Max amplitude orders found in LALSimulation (not accessible from outside of LALSim) */
66 #define MAX_NONPRECESSING_AMP_PN_ORDER 6
67 #define MAX_PRECESSING_AMP_PN_ORDER 3
68 
69 #define Pi_p2 9.8696044010893586188344909998761511
70 #define Pi_p2by3 2.1450293971110256000774441009412356
71 #define log4 1.3862943611198906188344642429163531
72 
73 static void q2masses(double mc, double q, double *m1, double *m2);
74 
75 /* list of testing GR parameters to be passed to the waveform */
76 /* 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 */
77 
78 
79 const char list_extra_parameters[76][16] = {"dchiMinus2","dchiMinus1","dchi0","dchi1","dchi2","dchi3","dchi3S","dchi3NS","dchi4","dchi4S","dchi4NS","dchi5","dchi5S","dchi5NS","dchi5l","dchi5lS","dchi5lNS","dchi6","dchi6S","dchi6NS","dchi6l","dchi7","dchi7S","dchi7NS","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","dchikappaS","dchikappaA","domega220","dtau220","domega210","dtau210","domega330","dtau330","domega440","dtau440","domega550","dtau550","db1","db2","db3","db4","dc1","dc2","dc4","dcl"};
80 
81 const UINT4 N_extra_params = 76;
82 
83 const char list_FTA_parameters[26][16] = {"dchiMinus2","dchiMinus1","dchi0","dchi1","dchi2","dchi3","dchi3S","dchi3NS","dchi4","dchi4S","dchi4NS","dchi5","dchi5S","dchi5NS","dchi5l","dchi5lS","dchi5lNS","dchi6","dchi6S","dchi6NS","dchi6l","dchi7","dchi7S","dchi7NS","dchikappaS","dchikappaA"};
84 
85 const UINT4 N_FTA_params = 26;
86 
87 
88 /* Return the quadrupole moment of a neutron star given its lambda
89  * We use the relations defined here. https://arxiv.org/pdf/1302.4499.pdf.
90  * Note that the convention we use is that:
91  * \f$\mathrm{dquadmon} = \bar{Q} - 1.\f$
92  * Where \f$\bar{Q}\f$ (dimensionless) is the reduced quadrupole moment.
93  */
94 static REAL8 dquadmon_from_lambda(REAL8 lambdav);
95 static REAL8 dquadmon_from_lambda(REAL8 lambdav)
96 {
97 
98  double ll = log(lambdav);
99  double ai = .194, bi = .0936, ci = 0.0474, di = -4.21e-3, ei = 1.23e-4;
100  double ln_quad_moment = ai + bi*ll + ci*ll*ll + di*pow(ll,3.0) + ei*pow(ll,4.0);
101  return(exp(ln_quad_moment) - 1.0);
102 }
103 
104 static int InterpolateWaveform(REAL8Vector *freqs, COMPLEX16FrequencySeries *src, COMPLEX16FrequencySeries *dest);
105 static int InterpolateWaveform(REAL8Vector *freqs, COMPLEX16FrequencySeries *src, COMPLEX16FrequencySeries *dest)
106 {
107  REAL8 deltaF = dest->deltaF;
108  UINT4 j=ceil(freqs->data[0] / deltaF);
109  COMPLEX16 *d=dest->data->data;
110  memset(d, 0, sizeof(*(d))*j);
111 
112  /* Loop over reduced frequency set */
113  for(UINT4 i=0;i<freqs->length-1;i++)
114  {
115  double startpsi = carg(src->data->data[i]);
116  double startamp = cabs(src->data->data[i]);
117  double endpsi = carg(src->data->data[i+1]);
118  double endamp = cabs(src->data->data[i+1]);
119 
120  double startf=freqs->data[i];
121  double endf=freqs->data[i+1];
122 
123  double df = endf - startf; /* Big freq step */
124 
125  /* linear interpolation setup */
126  double dpsi = (endpsi - startpsi);
127 
128  /* Catch case where phase wraps around */
129  /* NOTE: If changing this check that waveforms are not corrupted
130  * at high frequencies when dpsi/df can go slightly -ve without
131  * the phase wrapping around (e.g. TF2 1.4-1.4 srate=4096)
132  */
133  if (dpsi/df<-LAL_PI ) {dpsi+=LAL_TWOPI;}
134 
135  double dpsidf = dpsi/df;
136  double dampdf = (endamp - startamp)/df;
137 
138  double damp = dampdf *deltaF;
139 
140  const double dim = sin(dpsidf*deltaF);
141  const double dre = 2.0*sin(dpsidf*deltaF*0.5)*sin(dpsidf*deltaF*0.5);
142 
143  /* Loop variables */
144  double newRe,newIm,f,re,im,a;
145  for(f=j*deltaF,
146  re = cos(startpsi), im = sin(startpsi),
147  a = startamp;
148 
149  f<endf;
150 
151  j++, f+=deltaF,
152  newRe = re - dre*re-dim*im,
153  newIm = im + re*dim-dre*im,
154  re=newRe, im = newIm,
155  a += damp )
156  {
157  d[j] = a * (re + I*im);
158  }
159  }
160  memset(&(d[j]), 0, sizeof(d[j])*(dest->data->length - j));
161  return 0;
162 }
163 
164 
165 void LALInferenceTemplateNullFreqdomain(LALInferenceModel *model)
166 /**********************************************/
167 /* returns a frequency-domain 'null' template */
168 /* (all zeroes, implying no signal present). */
169 /**********************************************/
170 {
171  UINT4 i;
172  if ((model->freqhPlus==NULL) || (model->freqhCross==NULL)) {
173  XLALPrintError(" ERROR in templateNullFreqdomain(): encountered unallocated 'freqModelhPlus/-Cross'.\n");
174  XLAL_ERROR_VOID(XLAL_EFAULT);
175  }
176  for (i=0; i<model->freqhPlus->data->length; ++i){
177  model->freqhPlus->data->data[i] = 0.0;
178  model->freqhCross->data->data[i] = 0.0;
179  }
180  model->domain = LAL_SIM_DOMAIN_FREQUENCY;
181  return;
182 }
183 
184 
185 
186 void LALInferenceTemplateNullTimedomain(LALInferenceModel *model)
187 /*********************************************/
188 /* returns a time-domain 'null' template */
189 /* (all zeroes, implying no signal present). */
190 /*********************************************/
191 {
192  UINT4 i;
193  if ((model->timehPlus==NULL) || (model->timehCross==NULL)) {
194  XLALPrintError(" ERROR in templateNullTimedomain(): encountered unallocated 'timeModelhPlus/-Cross'.\n");
195  XLAL_ERROR_VOID(XLAL_EFAULT);
196  }
197  for (i=0; i<model->timehPlus->data->length; ++i){
198  model->timehPlus->data->data[i] = 0.0;
199  model->timehCross->data->data[i] = 0.0;
200  }
201  model->domain = LAL_SIM_DOMAIN_TIME;
202  return;
203 }
204 
205 
206 
207 /* ============ LAL template wrapper function: ========== */
208 
209 
210 static void mc2masses(double mc, double eta, double *m1, double *m2);
211 
212 static void mc2masses(double mc, double eta, double *m1, double *m2)
213 /* Compute individual companion masses (m1, m2) */
214 /* for given chirp mass (m_c) & mass ratio (eta) */
215 /* (note: m1 >= m2). */
216 {
217  double root = sqrt(0.25-eta);
218  double fraction = (0.5+root) / (0.5-root);
219  *m2 = mc * (pow(1+fraction,0.2) / pow(fraction,0.6));
220  *m1 = mc * (pow(1+1.0/fraction,0.2) / pow(1.0/fraction,0.6));
221  return;
222 }
223 
224 static void q2masses(double mc, double q, double *m1, double *m2)
225 /* Compute individual companion masses (m1, m2) */
226 /* for given chirp mass (m_c) & asymmetric mass */
227 /* ratio (q). note: q = m2/m1, where m1 >= m2 */
228 {
229  *m1 = mc * pow(q, -3.0/5.0) * pow(q+1, 1.0/5.0);
230  *m2 = (*m1) * q;
231  return;
232 }
233 
234 void LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence(LALInferenceModel *model){
235 /*************************************************************************************************************************/
236  Approximant approximant = (Approximant) 0;
237 
238  int ret=0;
239  INT4 errnum=0;
240 
241  model->roq->hptildeLinear=NULL, model->roq->hctildeLinear=NULL;
242  model->roq->hptildeQuadratic=NULL, model->roq->hctildeQuadratic=NULL;
243  REAL8 mc;
244  REAL8 phi0, m1, m2, distance, inclination;
245 
246  REAL8 *m1_p,*m2_p;
247 
248  if (LALInferenceCheckVariable(model->params, "LAL_APPROXIMANT"))
249  approximant = *(Approximant*) LALInferenceGetVariable(model->params, "LAL_APPROXIMANT");
250  else {
251  XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_APPROXIMANT\" parameter not provided!\n");
252  XLAL_ERROR_VOID(XLAL_EDATA);
253  }
254 
255  if (LALInferenceCheckVariable(model->params, "LAL_PNORDER"))
256  XLALSimInspiralWaveformParamsInsertPNPhaseOrder(model->LALpars, *(INT4 *) LALInferenceGetVariable(model->params, "LAL_PNORDER"));
257  else {
258  XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_PNORDER\" parameter not provided!\n");
259  XLAL_ERROR_VOID(XLAL_EDATA);
260  }
261 
262  /* Explicitly set the default amplitude order if one is not specified.
263  * This serves two purposes:
264  * 1) The default behavior of the code won't change unexpectedly due to changes in LALSimulation.
265  * 2) We need to know the amplitude order in order to set the starting frequency of the waveform properly. */
266  if (LALInferenceCheckVariable(model->params, "LAL_AMPORDER"))
267  XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars, *(INT4 *) LALInferenceGetVariable(model->params, "LAL_AMPORDER"));
268  else
269  if (!XLALSimInspiralWaveformParamsPNAmplitudeOrderIsDefault(model->LALpars))
270  XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,-1);
271 
272  REAL8 f_ref = 100.0;
273  if (LALInferenceCheckVariable(model->params, "f_ref")) f_ref = *(REAL8 *)LALInferenceGetVariable(model->params, "f_ref");
274 
275  REAL8 fTemp = f_ref;
276  if(LALInferenceCheckVariable(model->params,"chirpmass"))
277  {
278  mc = *(REAL8*) LALInferenceGetVariable(model->params, "chirpmass");
279  if (LALInferenceCheckVariable(model->params,"q")) {
280  REAL8 q = *(REAL8 *)LALInferenceGetVariable(model->params,"q");
281  q2masses(mc, q, &m1, &m2);
282  } else {
283  REAL8 eta = *(REAL8*) LALInferenceGetVariable(model->params, "eta");
284  mc2masses(mc, eta, &m1, &m2);
285  }
286  }
287  else if((m1_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass1")) && (m2_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass2")))
288  {
289  m1=*m1_p;
290  m2=*m2_p;
291  }
292  else
293  {
294  fprintf(stderr,"No mass parameters found!");
295  exit(0);
296  }
297  if(LALInferenceCheckVariable(model->params,"logdistance"))
298  {
299  distance = exp(LALInferenceGetREAL8Variable(model->params, "logdistance"))* LAL_PC_SI * 1.0e6; /* distance (1 Mpc) in units of metres */
300  }
301  else
302  {
303  distance = LAL_PC_SI * 1.0e6; /* 1Mpc default */
304  }
305 
306  phi0 = LALInferenceGetREAL8Variable(model->params, "phase"); /* START phase as per lalsimulation convention, radians*/
307  /* Zenith angle between J and N in radians. Also known as inclination angle when spins are aligned */
308  REAL8 thetaJN = acos(LALInferenceGetREAL8Variable(model->params, "costheta_jn")); /* zenith angle between J and N in radians */
309 
310  /* ==== SPINS ==== */
311  /* We will default to spinless signal and then add in the spin components if required */
312  /* If there are non-aligned spins, we must convert between the System Frame coordinates
313  * and the cartestian coordinates */
314 
315  /* The cartesian spin coordinates (default 0), as passed to LALSimulation */
316  REAL8 spin1x = 0.0;
317  REAL8 spin1y = 0.0;
318  REAL8 spin1z = 0.0;
319  REAL8 spin2x = 0.0;
320  REAL8 spin2y = 0.0;
321  REAL8 spin2z = 0.0;
322 
323  /* System frame coordinates as used for jump proposals */
324  REAL8 a_spin1 = 0.0; /* Magnitude of spin1 */
325  REAL8 a_spin2 = 0.0; /* Magnitude of spin2 */
326  REAL8 phiJL = 0.0; /* azimuthal angle of L_N on its cone about J radians */
327  REAL8 tilt1 = 0.0; /* zenith angle between S1 and LNhat in radians */
328  REAL8 tilt2 = 0.0; /* zenith angle between S2 and LNhat in radians */
329  REAL8 phi12 = 0.0; /* difference in azimuthal angle btwn S1, S2 in radians */
330 
331  /* Now check if we have spin amplitudes */
332  if(LALInferenceCheckVariable(model->params, "a_spin1")) a_spin1 = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin1");
333  if(LALInferenceCheckVariable(model->params, "a_spin2")) a_spin2 = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin2");
334 
335  /* Check if we have spin angles too */
336  if(LALInferenceCheckVariable(model->params, "phi_jl"))
337  phiJL = LALInferenceGetREAL8Variable(model->params, "phi_jl");
338  if(LALInferenceCheckVariable(model->params, "tilt_spin1"))
339  tilt1 = LALInferenceGetREAL8Variable(model->params, "tilt_spin1");
340  if(LALInferenceCheckVariable(model->params, "tilt_spin2"))
341  tilt2 = LALInferenceGetREAL8Variable(model->params, "tilt_spin2");
342  if(LALInferenceCheckVariable(model->params, "phi12"))
343  phi12 = LALInferenceGetREAL8Variable(model->params, "phi12");
344 
345  /* If we have tilt angles zero, then the spins are aligned and we just set the z component */
346  /* However, if the waveform supports precession then we still need to get the right coordinate components */
347  if(tilt1==0.0 && tilt2==0.0)
348  {
349  spin1z=a_spin1;
350  spin2z=a_spin2;
351  inclination = thetaJN; /* Inclination angle is just thetaJN */
352  }
353  else
354  { /* Template is not aligned-spin only. */
355  /* Set all the other spin components according to the angles we received above */
356  /* The transformation function doesn't know fLow, so f_ref==0 isn't interpretted as a request to use the starting frequency for reference. */
357  XLAL_TRY(ret=XLALSimInspiralTransformPrecessingNewInitialConditions(
358  &inclination, &spin1x, &spin1y, &spin1z, &spin2x, &spin2y, &spin2z,
359  thetaJN, phiJL, tilt1, tilt2, phi12, a_spin1, a_spin2, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, fTemp, phi0), errnum);
360  if (ret == XLAL_FAILURE)
361  {
362  XLALPrintError(" ERROR in XLALSimInspiralTransformPrecessingNewInitialConditions(): error converting angles. errnum=%d\n",errnum );
363  return;
364  }
365  }
366 /* ==== Spin induced quadrupole moment PARAMETERS ==== */
367  if(LALInferenceCheckVariable(model->params, "dQuadMonS")&&LALInferenceCheckVariable(model->params, "dQuadMonA")){
368  REAL8 dQuadMon1=0.;
369  REAL8 dQuadMon2=0.;
370  REAL8 dQuadMonS = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonS");
371  REAL8 dQuadMonA = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonA");
372  LALInferencedQuadMonSdQuadMonA(dQuadMonS,dQuadMonA,&dQuadMon1,&dQuadMon2);
373  XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
374  XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
375  }
376  else if(LALInferenceCheckVariable(model->params, "dQuadMon1")&&LALInferenceCheckVariable(model->params, "dQuadMon2")){
377  REAL8 dQuadMon1 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon1");
378  REAL8 dQuadMon2 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon2");
379  XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
380  XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
381  }
382  /* ==== TIDAL PARAMETERS ==== */
383  if(LALInferenceCheckVariable(model->params, "lambda1"))
384  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda1"));
385  if(LALInferenceCheckVariable(model->params, "lambda2"))
386  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda2"));
387  REAL8 lambdaT = 0.;
388  REAL8 dLambdaT = 0.;
389  REAL8 sym_mass_ratio_eta = 0.;
390 
391  if(LALInferenceCheckVariable(model->params, "lambdaT")&&LALInferenceCheckVariable(model->params, "dLambdaT")){
392  REAL8 lambda1=0.;
393  REAL8 lambda2=0.;
394  lambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "lambdaT");
395  dLambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "dLambdaT");
396  sym_mass_ratio_eta = m1*m2/((m1+m2)*(m1+m2));
397  LALInferenceLambdaTsEta2Lambdas(lambdaT,dLambdaT,sym_mass_ratio_eta,&lambda1,&lambda2);
398  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars,lambda1);
399  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars,lambda2);
400  }
401 
402  /* ==== 4-PIECE POLYTROPE EOS PARAMETERS ==== */
403  REAL8 logp1 = 0.;
404  REAL8 gamma1 = 0.;
405  REAL8 gamma2 = 0.;
406  REAL8 gamma3 = 0.;
407  if(LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "gamma1")&&LALInferenceCheckVariable(model->params, "gamma2")&&LALInferenceCheckVariable(model->params, "gamma3")){
408  REAL8 lambda1 = 0.;
409  REAL8 lambda2 = 0.;
410  logp1 = *(REAL8*) LALInferenceGetVariable(model->params, "logp1");
411  gamma1 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma1");
412  gamma2 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma2");
413  gamma3 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma3");
414  // Find lambda1,2(m1,2|eos)
415  LALInferenceLogp1GammasMasses2Lambdas(logp1,gamma1,gamma2,gamma3,m1,m2,&lambda1,&lambda2);
416  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
417  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
418  }
419 
420  /* ==== SPECTRAL DECOMPOSITION PARAMETERS ==== */
421  REAL8 SDgamma0 = 0.;
422  REAL8 SDgamma1 = 0.;
423  REAL8 SDgamma2 = 0.;
424  REAL8 SDgamma3 = 0.;
425  /* Checks for 4 spectral parameters */
426  if(!LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "SDgamma0")&&LALInferenceCheckVariable(model->params, "SDgamma1")&&LALInferenceCheckVariable(model->params, "SDgamma2")&&LALInferenceCheckVariable(model->params,"SDgamma3")){
427  REAL8 lambda1 = 0.;
428  REAL8 lambda2 = 0.;
429  SDgamma0 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma0");
430  SDgamma1 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma1");
431  SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
432  SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
433  REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
434  LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
435  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
436  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
437  }
438 
439  /* ==== BINARY_LOVE PARAMETERS ==== */
440  if(LALInferenceCheckVariable(model->params, "lambdaS")){
441  REAL8 lambda1=0.;
442  REAL8 lambda2=0.;
443  LALInferenceBinaryLove(model->params, &lambda1, &lambda2);
444  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
445  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
446  LALInferenceAddVariable(model->params, "lambda1", &lambda1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
447  LALInferenceAddVariable(model->params, "lambda2", &lambda2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
448  }
449 
450  /* Only use GR templates */
451  /* Fill in the extra parameters for testing GR, if necessary */
452  for (UINT4 k=0; k<N_extra_params; k++)
453  {
454  if(LALInferenceCheckVariable(model->params,list_extra_parameters[k]))
455  {
456  XLALDictInsert(model->LALpars, list_extra_parameters[k], (void *)LALInferenceGetVariable(model->params,list_extra_parameters[k]), sizeof(double), LAL_D_TYPE_CODE);
457 
458  //XLALSimInspiralAddTestGRParam(&nonGRparams,list_extra_parameters[k],*(REAL8 *)LALInferenceGetVariable(model->params,list_extra_parameters[k]));
459  }
460  }
461  /* Fill in PPE params if they are available */
462  char PPEparam[64]="";
463  const char *PPEnames[]={"aPPE","alphaPPE","bPPE","betaPPE",NULL};
464  for(UINT4 idx=0;PPEnames[idx];idx++)
465  {
466  for(UINT4 ppeidx=0;;ppeidx++)
467  {
468  sprintf(PPEparam, "%s%d",PPEnames[idx],ppeidx);
469  if(LALInferenceCheckVariable(model->params,PPEparam))
470  XLALDictInsert(model->LALpars, PPEparam, (void *)LALInferenceGetVariable(model->params,PPEparam), sizeof(double), LAL_D_TYPE_CODE);
471  else
472  break;
473  }
474  }
475 
476  /* ==== Call the waveform generator ==== */
477  /* Correct distance to account for renormalisation of data due to window RMS */
478  double corrected_distance = distance * sqrt(model->window->sumofsquares/model->window->data->length);
479  XLAL_TRY(ret=XLALSimInspiralChooseFDWaveformSequence (&(model->roq->hptildeLinear), &(model->roq->hctildeLinear), phi0, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI,
480  spin1x, spin1y, spin1z, spin2x, spin2y, spin2z, f_ref, corrected_distance, inclination, model->LALpars, approximant, (model->roq->frequencyNodesLinear)), errnum);
481 
482  XLAL_TRY(ret=XLALSimInspiralChooseFDWaveformSequence (&(model->roq->hptildeQuadratic), &(model->roq->hctildeQuadratic), phi0, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI,
483  spin1x, spin1y, spin1z, spin2x, spin2y, spin2z, f_ref, corrected_distance, inclination, model->LALpars, approximant, (model->roq->frequencyNodesQuadratic)), errnum);
484 
485  REAL8 instant = model->freqhPlus->epoch.gpsSeconds + 1e-9*model->freqhPlus->epoch.gpsNanoSeconds;
486  LALInferenceSetVariable(model->params, "time", &instant);
487 
488  return;
489 }
490 
491 void LALInferenceTemplateSineGaussian(LALInferenceModel *model)
492 /*****************************************************/
493 /* Sine-Gaussian (burst) template. */
494 /* Signal is (by now?) linearly polarised, */
495 /* i.e., the cross-waveform remains zero. */
496 /* * * * * * * * * * * * * * * * * * * * * * * * * * */
497 /* The (plus-) waveform is: */
498 /* a * exp(-((t-mu)/sigma)^2) * sin(2*pi*f*t-phi) */
499 /* Note that by setting f=0, phi=pi/2 you also get */
500 /* a `pure' Gaussian template. */
501 /* */
502 /* * * * * * * * * * * * * * * * * * * * * * * * * * ************************************/
503 /* Required (`model->params') parameters are: */
504 /* - "time" (the "mu" parameter of the Gaussian part; REAL8, GPS sec.) */
505 /* - "sigma" (width, the "sigma" parameter of the Gaussian part; REAL8, seconds) */
506 /* - "frequency" (frequency of the sine part; REAL8, Hertz) */
507 /* - "phase" (phase (at above "mu"); REAL8, radians) */
508 /* - "amplitude" (amplitude, REAL8) */
509 /****************************************************************************************/
510 {
511  double endtime = *(REAL8*) LALInferenceGetVariable(model->params, "time"); /* time parameter ("mu"), GPS sec. */
512  double sigma = *(REAL8*) LALInferenceGetVariable(model->params, "sigma"); /* width parameter, seconds */
513  double f = *(REAL8*) LALInferenceGetVariable(model->params, "frequency"); /* frequency, Hz */
514  double phi = *(REAL8*) LALInferenceGetVariable(model->params, "phase"); /* phase, rad */
515  double a = *(REAL8*) LALInferenceGetVariable(model->params, "amplitude"); /* amplitude */
516  double t, tsigma, twopif = LAL_TWOPI*f;
517  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
518  unsigned long i;
519  if (sigma <= 0.0) {
520  fprintf(stderr, " ERROR in templateSineGaussian(): zero or negative \"sigma\" parameter (sigma=%e).\n", sigma);
521  exit(1);
522  }
523  if (f < 0.0)
524  fprintf(stderr, " WARNING in templateSineGaussian(): negative \"frequency\" parameter (f=%e).\n", f);
525  if (a < 0.0)
526  fprintf(stderr, " WARNING in templateSineGaussian(): negative \"amplitude\" parameter (a=%e).\n", a);
527  for (i=0; i<model->timehPlus->data->length; ++i){
528  t = ((double)i)*model->deltaT + (epochGPS-endtime); /* t-mu */
529  tsigma = t/sigma; /* (t-mu)/sigma */
530  if (fabs(tsigma) < 5.0) /* (only do computations within a 10 sigma range) */
531  model->timehPlus->data->data[i] = a * exp(-0.5*tsigma*tsigma) * sin(twopif*t+phi);
532  else
533  model->timehPlus->data->data[i] = 0.0;
534  model->timehCross->data->data[i] = 0.0;
535  }
536  model->domain = LAL_SIM_DOMAIN_TIME;
537  return;
538 }
539 
540 void LALInferenceTemplateDampedSinusoid(LALInferenceModel *model)
541 /*****************************************************/
542 /* Damped Sinusoid (burst) template. */
543 /* Signal is linearly polarized, */
544 /* i.e., cross term is zero. */
545 /* * * * * * * * * * * * * * * * * * * * * * * * * * */
546 /* The (plus-) waveform is an exponentially decaying */
547 /* sine wave: */
548 /* a * exp((t-time)/tau) * sin(2*pi*f*(t-time)) */
549 /* where "time" is the time parameter denoting the */
550 /* instant where the signal starts. */
551 /* * * * * * * * * * * * * * * * * * * * * * * * * * **************************/
552 /* Required (`model->params') parameters are: */
553 /* - "time" (the instant at which the signal starts; REAL8, GPS sec.) */
554 /* - "tau" (width parameter; REAL8, seconds) */
555 /* - "frequency" (frequency of the sine part; REAL8, Hertz) */
556 /* - "amplitude" (amplitude, REAL8) */
557 /******************************************************************************/
558 {
559  double endtime = *(REAL8*) LALInferenceGetVariable(model->params, "time"); /* time parameter ("mu"), GPS sec. */
560  double tau = *(REAL8*) LALInferenceGetVariable(model->params, "tau"); /* width parameter, seconds */
561  double f = *(REAL8*) LALInferenceGetVariable(model->params, "frequency"); /* frequency, Hz */
562  double a = *(REAL8*) LALInferenceGetVariable(model->params, "amplitude"); /* amplitude */
563  double t, ttau, twopif = LAL_TWOPI*f;
564  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
565  unsigned long i;
566  if (tau <= 0.0) {
567  fprintf(stderr, " ERROR in templateDampedSinusoid(): zero or negative \"tau\" parameter (tau=%e).\n", tau);
568  exit(1);
569  }
570  if (f < 0.0)
571  fprintf(stderr, " WARNING in templateDampedSinusoid(): negative \"frequency\" parameter (f=%e).\n", f);
572  for (i=0; i<model->timehPlus->data->length; ++i){
573  t = ((double)i)*model->deltaT + (epochGPS-endtime); /* t-mu */
574  if ((t>0.0) && ((ttau=t/tau) < 10.0)) /* (only do computations within a 10 tau range) */
575  model->timehPlus->data->data[i] = a * exp(-ttau) * sin(twopif*t);
576  else
577  model->timehPlus->data->data[i] = 0.0;
578  model->timehCross->data->data[i] = 0.0;
579  }
580  model->domain = LAL_SIM_DOMAIN_TIME;
581  return;
582 }
583 
584 
585 
586 void LALInferenceTemplateSinc(LALInferenceModel *model)
587 /*****************************************************/
588 /* Sinc function (burst) template. */
589 /* Signal is linearly polarized, */
590 /* i.e., cross term is zero. */
591 /* * * * * * * * * * * * * * * * * * * * * * * * * * */
592 /* The (plus-) waveform is a sinc function of given */
593 /* frequency: */
594 /* a * sinc(2*pi*f*(t-time)) */
595 /* = a * sin(2*pi*f*(t-time)) / (2*pi*f*(t-time)) */
596 /* where "time" is the time parameter denoting the */
597 /* signal's central peak location. */
598 /* * * * * * * * * * * * * * * * * * * * * * * * * * *************************/
599 /* Required (`model->params') parameters are: */
600 /* - "time" (the instant at which the signal peaks; REAL8, GPS sec.) */
601 /* - "frequency" (frequency of the sine part; REAL8, Hertz) */
602 /* - "amplitude" (amplitude, REAL8) */
603 /*****************************************************************************/
604 {
605  double endtime = *(REAL8*) LALInferenceGetVariable(model->params, "time"); /* time parameter ("mu"), GPS sec. */
606  double f = *(REAL8*) LALInferenceGetVariable(model->params, "frequency"); /* frequency, Hz */
607  double a = *(REAL8*) LALInferenceGetVariable(model->params, "amplitude"); /* amplitude */
608  double t, sinArg, sinc, twopif = LAL_TWOPI*f;
609  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
610  unsigned long i;
611  if (f < 0.0)
612  fprintf(stderr, " WARNING in templateSinc(): negative \"frequency\" parameter (f=%e).\n", f);
613  for (i=0; i<model->timehPlus->data->length; ++i){
614  t = ((double)i)*model->deltaT + (epochGPS-endtime); /* t-mu */
615  sinArg = twopif*t;
616  sinc = (sinArg==0.0) ? 1.0 : sin(sinArg)/sinArg;
617  model->timehPlus->data->data[i] = a * sinc;
618  model->timehCross->data->data[i] = 0.0;
619  }
620  model->domain = LAL_SIM_DOMAIN_TIME;
621  return;
622 }
623 
624 
625 void LALInferenceTemplateASinOmegaT(LALInferenceModel *model)
626 /************************************************************/
627 /* Trivial h(t)=A*sin(Omega*t) template */
628 /* Required (`model->params') parameters are: */
629 /* - "A" (dimensionless amplitude, REAL8) */
630 /* - "Omega" (frequency; REAL8, radians/sec) */
631 /************************************************************/
632 {
633  double A = *(REAL8*) LALInferenceGetVariable(model->params, "A"); /* dim-less */
634  double Omega = *(REAL8*) LALInferenceGetVariable(model->params, "Omega"); /* rad/sec */
635  double t;
636  double epochGPS = XLALGPSGetREAL8(&(model->timehPlus->epoch));
637 
638  unsigned long i;
639  for (i=0; i<model->timehPlus->data->length; ++i){
640  t = ((double)i)*model->deltaT + (epochGPS); /* t-mu */
641  model->timehPlus->data->data[i] = A * sin(Omega*t);
642  model->timehCross->data->data[i] = 0.0;
643  }
644  model->domain = LAL_SIM_DOMAIN_TIME;
645  return;
646 }
647 
648 void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
649 /*************************************************************************************************************************/
650 /* Wrapper for LALSimulation waveforms: */
651 /* XLALSimInspiralChooseFDWaveform() and XLALSimInspiralChooseTDWaveform(). */
652 /* */
653 /* model->params parameters are: */
654 /* - "name" description; type OPTIONAL (default value) */
655 /* */
656 /* MODEL PARAMETERS */
657 /* - "LAL_APPROXIMANT Approximant; Approximant */
658 /* - "LAL_PNORDER" Phase PN order; INT4 */
659 /* - "LAL_AMPORDER" Amplitude PN order; INT4 OPTIONAL (-1) */
660 /* - "spinO" Spin order; LALSimInspiralSpinOrder OPTIONAL (LAL_SIM_INSPIRAL_SPIN_ORDER_DEFAULT) */
661 /* - "tideO" Tidal order; LALSimInspiralTidalOrder OPTIONAL (LAL_SIM_INSPIRAL_TIDAL_ORDER_DEFAULT) */
662 /* - "f_ref" frequency at which the (frequency dependent) parameters are defined; REAL8 OPTIONAL (0.0) */
663 /* - "fLow" lower frequency bound; REAL8 OPTIONAL (model->fLow) */
664 /* */
665 /* MASS PARAMETERS; either: */
666 /* - "mass1" mass of object 1 in solar mass; REAL8 */
667 /* - "mass2" mass of object 1 in solar mass; REAL8 */
668 /* OR */
669 /* - "chirpmass" chirpmass in solar mass; REAL8 */
670 /* - "q" asymmetric mass ratio m2/m1, 0<q<1; REAL8 */
671 /* OR */
672 /* - "chirpmass" chirpmass in solar mass; REAL8 */
673 /* - "eta" symmetric mass ratio (m1*m2)/(m1+m2)^2; REAL8 */
674 /* */
675 /* ORIENTATION AND SPIN PARAMETERS */
676 /* - "phi0" reference phase as per LALSimulation convention; REAL8 */
677 /* - "logdistance" distance in Mpc */
678 /* - "costheta_jn"); cos of zenith angle between J and N in radians; REAL8 */
679 /* - "phi_jl"); azimuthal angle of L_N on its cone about J radians; REAL8 */
680 /* - "tilt_spin1"); zenith angle between S1 and LNhat in radians; REAL8 */
681 /* - "tilt_spin2"); zenith angle between S2 and LNhat in radians; REAL8 */
682 /* - "phi12"); difference in azimuthal angle between S1, S2 in radians; REAL8 */
683 /* - "a_spin1" magnitude of spin 1 in general configuration, -1<a_spin1<1; REAL8 OPTIONAL (0.0) */
684 /* - "a_spin2" magnitude of spin 2 in general configuration, -1<a_spin1<1; REAL8 OPTIONAL (0.0) */
685 /* */
686 /* OTHER PARAMETERS */
687 /* - "lambda1" tidal parameter of object 1; REAL8 OPTIONAL (0.0) */
688 /* - "lambda2" tidal parameter of object 1; REAL8 OPTIONAL (0.0) */
689 /* */
690 /* - "time" used as an OUTPUT only; REAL8 */
691 /* */
692 /* */
693 /* model needs to also contain: */
694 /* - model->fLow Unless - "fLow" OPTIONAL */
695 /* - model->deltaT */
696 /* - if model->domain == LAL_SIM_DOMAIN_FREQUENCY */
697 /* - model->deltaF */
698 /* - model->freqhCross */
699 /* - model->freqhPlus */
700 /* - else */
701 /* - model->timehPlus */
702 /* - model->timehCross */
703 /*************************************************************************************************************************/
704 {
705 
706  Approximant approximant = (Approximant) 0;
707 
708  INT4 generic_fd_correction;
709 
710  static int sizeWarning = 0;
711  int ret=0;
712  INT4 errnum=0;
713 
714  REAL8TimeSeries *hplus=NULL; /**< +-polarization waveform [returned] */
715  REAL8TimeSeries *hcross=NULL; /**< x-polarization waveform [returned] */
716  COMPLEX16FrequencySeries *hptilde=NULL, *hctilde=NULL;
717  REAL8 mc;
718  REAL8 phi0, deltaT, m1, m2, f_low, f_start, distance, inclination;
719 
720  REAL8 *m1_p,*m2_p;
721  REAL8 deltaF, f_max;
722 
723  /* Sampling rate for time domain models */
724  deltaT = model->deltaT;
725 
726  if (LALInferenceCheckVariable(model->params, "LAL_APPROXIMANT"))
727  approximant = *(Approximant*) LALInferenceGetVariable(model->params, "LAL_APPROXIMANT");
728  else {
729  XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_APPROXIMANT\" parameter not provided!\n");
730  XLAL_ERROR_VOID(XLAL_EDATA);
731  }
732 
733  if (LALInferenceCheckVariable(model->params, "LAL_PNORDER"))
734  XLALSimInspiralWaveformParamsInsertPNPhaseOrder(model->LALpars, *(INT4*) LALInferenceGetVariable(model->params, "LAL_PNORDER"));
735  else {
736  XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_PNORDER\" parameter not provided!\n");
737  XLAL_ERROR_VOID(XLAL_EDATA);
738  }
739 
740  /* Explicitly set the default amplitude order if one is not specified.
741  * This serves two purposes:
742  * 1) The default behavior of the code won't change unexpectedly due to changes in LALSimulation.
743  * 2) We need to know the amplitude order in order to set the starting frequency of the waveform properly. */
744  if (LALInferenceCheckVariable(model->params, "LAL_AMPORDER"))
745  XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,*(INT4*) LALInferenceGetVariable(model->params, "LAL_AMPORDER"));
746  else
747  if (!XLALSimInspiralWaveformParamsPNAmplitudeOrderIsDefault(model->LALpars))
748  XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,-1);
749 
750  /* Check if the user requested modify waveforms with FTA */
751  if (LALInferenceCheckVariable(model->params, "generic_fd_correction"))
752  generic_fd_correction = *(INT4*) LALInferenceGetVariable(model->params, "generic_fd_correction");
753  else
754  generic_fd_correction = 0;
755 
756  REAL8 f_ref = 100.0;
757  if (LALInferenceCheckVariable(model->params, "f_ref")) f_ref = *(REAL8 *)LALInferenceGetVariable(model->params, "f_ref");
758 
759  REAL8 fTemp = f_ref;
760 
761  if(LALInferenceCheckVariable(model->params,"chirpmass"))
762  {
763  mc = *(REAL8*) LALInferenceGetVariable(model->params, "chirpmass");
764  if (LALInferenceCheckVariable(model->params,"q")) {
765  REAL8 q = *(REAL8 *)LALInferenceGetVariable(model->params,"q");
766  q2masses(mc, q, &m1, &m2);
767  } else {
768  REAL8 eta = *(REAL8*) LALInferenceGetVariable(model->params, "eta");
769  mc2masses(mc, eta, &m1, &m2);
770  }
771  }
772  else if((m1_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass1")) && (m2_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass2")))
773  {
774  m1=*m1_p;
775  m2=*m2_p;
776  }
777  else
778  {
779  fprintf(stderr,"No mass parameters found!");
780  exit(0);
781  }
782 
783  if(!LALInferenceCheckVariable(model->params,"logdistance")) distance=LAL_PC_SI * 1e6; /* If distance not given, 1Mpc used */
784  else
785  {
786  distance = exp(LALInferenceGetREAL8Variable(model->params, "logdistance"))* LAL_PC_SI * 1.0e6; /* distance (1 Mpc) in units of metres */
787  }
788  phi0 = LALInferenceGetREAL8Variable(model->params, "phase"); /* START phase as per lalsimulation convention, radians*/
789 
790  /* Zenith angle between J and N in radians. Also known as inclination angle when spins are aligned */
791  REAL8 thetaJN = acos(LALInferenceGetREAL8Variable(model->params, "costheta_jn")); /* zenith angle between J and N in radians */
792 
793  /* Check if fLow is a model parameter, otherwise use data structure definition */
794  if(LALInferenceCheckVariable(model->params, "flow"))
795  f_low = *(REAL8*) LALInferenceGetVariable(model->params, "flow");
796  else
797  f_low = model->fLow;
798 
799  f_start = XLALSimInspiralfLow2fStart(f_low, XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(model->LALpars), approximant);
800 
801  /* Don't let TaylorF2 generate unphysical inspiral up to Nyquist */
802  if (approximant == TaylorF2)
803  f_max = 0.0; /* this will stop at ISCO */
804  else
805  f_max = model->fHigh; /* this will be the highest frequency used across the network */
806 
807  /* ==== SPINS ==== */
808  /* We will default to spinless signal and then add in the spin components if required */
809  /* If there are non-aligned spins, we must convert between the System Frame coordinates
810  * and the cartestian coordinates */
811 
812  /* The cartesian spin coordinates (default 0), as passed to LALSimulation */
813  REAL8 spin1x = 0.0;
814  REAL8 spin1y = 0.0;
815  REAL8 spin1z = 0.0;
816  REAL8 spin2x = 0.0;
817  REAL8 spin2y = 0.0;
818  REAL8 spin2z = 0.0;
819 
820  /* System frame coordinates as used for jump proposals */
821  REAL8 a_spin1 = 0.0; /* Magnitude of spin1 */
822  REAL8 a_spin2 = 0.0; /* Magnitude of spin2 */
823  REAL8 phiJL = 0.0; /* azimuthal angle of L_N on its cone about J radians */
824  REAL8 tilt1 = 0.0; /* zenith angle between S1 and LNhat in radians */
825  REAL8 tilt2 = 0.0; /* zenith angle between S2 and LNhat in radians */
826  REAL8 phi12 = 0.0; /* difference in azimuthal angle btwn S1, S2 in radians */
827 
828  /* Now check if we have spin amplitudes */
829  if(LALInferenceCheckVariable(model->params, "a_spin1")) a_spin1 = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin1");
830  if(LALInferenceCheckVariable(model->params, "a_spin2")) a_spin2 = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin2");
831 
832  /* Check if we have spin angles too */
833  if(LALInferenceCheckVariable(model->params, "phi_jl"))
834  phiJL = LALInferenceGetREAL8Variable(model->params, "phi_jl");
835  if(LALInferenceCheckVariable(model->params, "tilt_spin1"))
836  tilt1 = LALInferenceGetREAL8Variable(model->params, "tilt_spin1");
837  if(LALInferenceCheckVariable(model->params, "tilt_spin2"))
838  tilt2 = LALInferenceGetREAL8Variable(model->params, "tilt_spin2");
839  if(LALInferenceCheckVariable(model->params, "phi12"))
840  phi12 = LALInferenceGetREAL8Variable(model->params, "phi12");
841 
842  /* If we have tilt angles zero, then the spins are aligned and we just set the z component */
843  /* However, if the waveform supports precession then we still need to get the right coordinate components */
844  if(tilt1==0.0 && tilt2==0.0)
845  {
846  spin1z=a_spin1;
847  spin2z=a_spin2;
848  inclination = thetaJN; /* Inclination angle is just thetaJN */
849  }
850  else
851  { /* Template is not aligned-spin only. */
852  /* Set all the other spin components according to the angles we received above */
853  /* The transformation function doesn't know fLow, so f_ref==0 isn't interpretted as a request to use the starting frequency for reference. */
854  if(fTemp==0.0)
855  fTemp = f_start;
856  // If we have a *time-domain* AND *precessing* EOB template, set fTemp to f_start.
857  // This is done since EOB only uses f_start and ignores f_ref.
858  if(model->domain == LAL_SIM_DOMAIN_TIME && (strstr(XLALSimInspiralGetStringFromApproximant(approximant),"EOB") != NULL)){
859  fTemp = f_start;
860  }
861  XLAL_TRY(ret=XLALSimInspiralTransformPrecessingNewInitialConditions(
862  &inclination, &spin1x, &spin1y, &spin1z, &spin2x, &spin2y, &spin2z,
863  thetaJN, phiJL, tilt1, tilt2, phi12, a_spin1, a_spin2, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, fTemp, phi0), errnum);
864  if (ret == XLAL_FAILURE)
865  {
866  XLALPrintError(" ERROR in XLALSimInspiralTransformPrecessingNewInitialConditions(): error converting angles. errnum=%d: %s\n",errnum, XLALErrorString(errnum) );
867  return;
868  }
869  }
870 /* ==== Spin induced quadrupole moment PARAMETERS ==== */
871  if(LALInferenceCheckVariable(model->params, "dQuadMonS")&&LALInferenceCheckVariable(model->params, "dQuadMonA")){
872  REAL8 dQuadMon1=0.;
873  REAL8 dQuadMon2=0.;
874  REAL8 dQuadMonS = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonS");
875  REAL8 dQuadMonA = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonA");
876  LALInferencedQuadMonSdQuadMonA(dQuadMonS,dQuadMonA,&dQuadMon1,&dQuadMon2);
877  XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
878  XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
879  }
880  else if(LALInferenceCheckVariable(model->params, "dQuadMon1")&&LALInferenceCheckVariable(model->params, "dQuadMon2")){
881  REAL8 dQuadMon1 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon1");
882  REAL8 dQuadMon2 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon2");
883  XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
884  XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
885  }
886 /* ==== TIDAL PARAMETERS ==== */
887  if(LALInferenceCheckVariable(model->params, "lambda1"))
888  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda1"));
889  if(LALInferenceCheckVariable(model->params, "lambda2"))
890  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda2"));
891  REAL8 lambdaT = 0.;
892  REAL8 dLambdaT = 0.;
893  REAL8 sym_mass_ratio_eta = 0.;
894  if(LALInferenceCheckVariable(model->params, "lambdaT")&&LALInferenceCheckVariable(model->params, "dLambdaT")){
895  REAL8 lambda1=0.;
896  REAL8 lambda2=0.;
897  lambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "lambdaT");
898  dLambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "dLambdaT");
899  sym_mass_ratio_eta = m1*m2/((m1+m2)*(m1+m2));
900  LALInferenceLambdaTsEta2Lambdas(lambdaT,dLambdaT,sym_mass_ratio_eta,&lambda1,&lambda2);
901  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
902  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
903  }
904  if(model->eos_fam)
905  {
906  LALSimNeutronStarFamily *eos_fam = model->eos_fam;
907  REAL8 r1=0, r2=0, k2_1=0, k2_2=0, lambda1=0, lambda2=0;
908  REAL8 mass_max = XLALSimNeutronStarMaximumMass(eos_fam) / LAL_MSUN_SI;
909  REAL8 mass_min = XLALSimNeutronStarFamMinimumMass(eos_fam) / LAL_MSUN_SI;
910 
911  if(m1<mass_max && m1>mass_min)
912  {
913  /* Compute l1, l2 from mass and EOS */
914  r1 = XLALSimNeutronStarRadius(m1*LAL_MSUN_SI, eos_fam);
915  k2_1 = XLALSimNeutronStarLoveNumberK2(m1*LAL_MSUN_SI, eos_fam);
916  lambda1 = (2./3.)*k2_1 * pow(r1/(m1*LAL_MRSUN_SI), 5.0);
917  }
918  if(m2<mass_max && m2>mass_min)
919  {
920  r2 = XLALSimNeutronStarRadius(m2*LAL_MSUN_SI, eos_fam);
921  k2_2 = XLALSimNeutronStarLoveNumberK2(m2*LAL_MSUN_SI, eos_fam);
922  lambda2 = (2./3.)*k2_2 * pow(r2/(m2*LAL_MRSUN_SI), 5.0);
923  }
924  /* Set waveform params */
925  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
926  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
927  REAL8 dQuadMon1 = dquadmon_from_lambda(lambda1);
928  REAL8 dQuadMon2 = dquadmon_from_lambda(lambda2);
929  XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars, dQuadMon1);
930  XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars, dQuadMon2);
931 
932  /* Calculate maximum frequency */
933  /* If both lambdas are non-zero compute EOS-dependent f_max */
934  if((lambda1 > 0) && (lambda2 > 0) && (approximant == TaylorF2))
935  {
936  /* Start with ISCO */
937  f_max = 1. / (pow(6,1.5) * LAL_PI * (m1*LAL_MTSUN_SI + m2*LAL_MTSUN_SI));
938  REAL8 log_lambda1 = log(lambda1);
939  REAL8 log_lambda2 = log(lambda2);
940  REAL8 compactness1 = 0.371 - 0.0391 * log_lambda1 + 0.001056 * log_lambda1 * log_lambda1;
941  REAL8 compactness2 = 0.371 - 0.0391 * log_lambda2 + 0.001056 * log_lambda2 * log_lambda2;
942  REAL8 rad1 = m1*LAL_MTSUN_SI / compactness1;
943  REAL8 rad2 = m2*LAL_MTSUN_SI / compactness2;
944  /* Use the smaller of ISCO and the EOS-dependent termination */
945  REAL8 fmax_eos = 1. / LAL_PI * pow((m1*LAL_MTSUN_SI + m2*LAL_MTSUN_SI) / pow((rad1 + rad2),3.0),0.5);
946  if (fmax_eos < f_max)
947  {
948  f_max = fmax_eos;
949  }
950  }
951 
952  /* Add derived quantities for output */
953  LALInferenceAddVariable(model->params, "radius1", &r1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
954  LALInferenceAddVariable(model->params, "radius2", &r2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
955  LALInferenceAddVariable(model->params, "lambda1", &lambda1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
956  LALInferenceAddVariable(model->params, "lambda2", &lambda2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
957  }
958 
959  /* ==== 4-PIECE POLYTROPE EOS PARAMETERS ==== */
960  REAL8 logp1 = 0.;
961  REAL8 gamma1 = 0.;
962  REAL8 gamma2 = 0.;
963  REAL8 gamma3 = 0.;
964  if(LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "gamma1")&&LALInferenceCheckVariable(model->params, "gamma2")&&LALInferenceCheckVariable(model->params, "gamma3")){
965  REAL8 lambda1 = 0.;
966  REAL8 lambda2 = 0.;
967  logp1 = *(REAL8*) LALInferenceGetVariable(model->params, "logp1");
968  gamma1 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma1");
969  gamma2 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma2");
970  gamma3 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma3");
971  // Find lambda1,2(m1,2|eos)
972  LALInferenceLogp1GammasMasses2Lambdas(logp1,gamma1,gamma2,gamma3,m1,m2,&lambda1,&lambda2);
973  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
974  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
975  }
976 
977  /* ==== SPECTRAL DECOMPOSITION PARAMETERS ==== */
978  REAL8 SDgamma0 = 0.;
979  REAL8 SDgamma1 = 0.;
980  REAL8 SDgamma2 = 0.;
981  REAL8 SDgamma3 = 0.;
982  /* Checks for 4 spectral parameters */
983  if(!LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "SDgamma0")&&LALInferenceCheckVariable(model->params, "SDgamma1")&&LALInferenceCheckVariable(model->params, "SDgamma2")&&LALInferenceCheckVariable(model->params,"SDgamma3")){
984  REAL8 lambda1 = 0.;
985  REAL8 lambda2 = 0.;
986  SDgamma0 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma0");
987  SDgamma1 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma1");
988  SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
989  SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
990  REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
991  LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
992  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
993  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
994  }
995 
996  /* ==== BINARY_LOVE PARAMETERS ==== */
997  if(LALInferenceCheckVariable(model->params, "lambdaS")){
998  REAL8 lambda1=0.;
999  REAL8 lambda2=0.;
1000  LALInferenceBinaryLove(model->params, &lambda1, &lambda2);
1001  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
1002  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
1003  LALInferenceAddVariable(model->params, "lambda1", &lambda1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
1004  LALInferenceAddVariable(model->params, "lambda2", &lambda2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
1005  }
1006 
1007  /* Only use GR templates */
1008  /* Fill in the extra parameters for testing GR, if necessary */
1009  for (UINT4 k=0; k<N_extra_params; k++)
1010  {
1011  if(LALInferenceCheckVariable(model->params,list_extra_parameters[k]))
1012  {
1013  XLALDictInsert(model->LALpars, list_extra_parameters[k], (void *)LALInferenceGetVariable(model->params,list_extra_parameters[k]), sizeof(double), LAL_D_TYPE_CODE);
1014  }
1015  }
1016  /* Fill in PPE params if they are available */
1017  char PPEparam[64]="";
1018  const char *PPEnames[]={"aPPE","alphaPPE","bPPE","betaPPE",NULL};
1019  for(UINT4 idx=0;PPEnames[idx];idx++)
1020  {
1021  for(UINT4 ppeidx=0;;ppeidx++)
1022  {
1023  sprintf(PPEparam, "%s%d",PPEnames[idx],ppeidx);
1024  if(LALInferenceCheckVariable(model->params,PPEparam))
1025  XLALDictInsert(model->LALpars, PPEparam, (void *)LALInferenceGetVariable(model->params,PPEparam), sizeof(double), LAL_D_TYPE_CODE);
1026  else
1027  break;
1028  }
1029  }
1030 
1031 
1032 
1033  /* ==== Call the waveform generator ==== */
1034  if(model->domain == LAL_SIM_DOMAIN_FREQUENCY) {
1035  deltaF = model->deltaF;
1036  /* Correct distance to account for renormalisation of data due to window RMS */
1037  double corrected_distance = distance * sqrt(model->window->sumofsquares/model->window->data->length);
1038 
1039  /* check if the user requested a generic phase correction */
1040  if (generic_fd_correction != 1)
1041  {
1042  XLAL_TRY(ret=XLALSimInspiralChooseFDWaveformFromCache(&hptilde, &hctilde, phi0,
1043  deltaF, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, spin1x, spin1y, spin1z,
1044  spin2x, spin2y, spin2z, f_start, f_max, f_ref, corrected_distance, inclination, model->LALpars,
1045  approximant,model->waveformCache, NULL), errnum);
1046  }
1047  else
1048  {
1049 
1050  /* create a duplicate of LALpars */
1051  LALDict *grParams = XLALCreateDict();
1052  LALDictEntry *param;
1053  LALDictIter iter;
1054  XLALDictIterInit(&iter, model->LALpars);
1055  while ((param = XLALDictIterNext(&iter)) != NULL) {
1056  XLALDictInsertValue( grParams , XLALDictEntryGetKey(param), XLALDictEntryGetValue(param) );
1057  }
1058 
1059  /* set all FTA parameters in grParams to their default GR value*/
1060  REAL8 default_GR_value = 0.0;
1061  for (UINT4 k=0; k<N_FTA_params; k++)
1062  {
1063  XLALDictInsert(grParams, list_FTA_parameters[k], (void *) (&default_GR_value), sizeof(double), LAL_D_TYPE_CODE);
1064  }
1065 
1066  REAL8 correction_window = 1.;
1067  if(LALInferenceCheckVariable(model->params, "correction_window"))
1068  correction_window = *(REAL8*) LALInferenceGetVariable(model->params, "correction_window");
1069 
1070  REAL8 correction_ncycles_taper = 1.;
1071  if(LALInferenceCheckVariable(model->params, "correction_ncycles_taper"))
1072  correction_ncycles_taper = *(REAL8*) LALInferenceGetVariable(model->params, "correction_ncycles_taper");
1073 
1074  /* check if HM version of non-GR correction should be used */
1075  INT4 generic_fd_correction_hm = 0;
1076  if (LALInferenceCheckVariable(model->params, "generic_fd_correction_hm"))
1077  generic_fd_correction_hm = *(INT4*) LALInferenceGetVariable(model->params, "generic_fd_correction_hm");
1078 
1079  if(generic_fd_correction_hm == 0){
1080  /* generate waveform with non-GR parameters set to default (zero) */
1081  XLAL_TRY(ret=XLALSimInspiralChooseFDWaveformFromCache(&hptilde, &hctilde, phi0,
1082  deltaF, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, spin1x, spin1y, spin1z,
1083  spin2x, spin2y, spin2z, f_start, f_max, f_ref, corrected_distance, inclination, grParams,
1084  approximant,model->waveformCache, NULL), errnum);
1085 
1086  /* apply FTA corrections */
1087  XLAL_TRY(ret = XLALSimInspiralTestingGRCorrections(hptilde,2,2,m1*LAL_MSUN_SI,m2*LAL_MSUN_SI, spin1z, spin2z, f_start, f_ref, correction_window, correction_ncycles_taper, model->LALpars), errnum);
1088 
1089  XLAL_TRY(ret = XLALSimInspiralTestingGRCorrections(hctilde,2,2,m1*LAL_MSUN_SI,m2*LAL_MSUN_SI, spin1z, spin2z, f_start, f_ref, correction_window, correction_ncycles_taper, model->LALpars), errnum);
1090  }
1091  else{
1092  /* generate modes with non-GR parameters set to default (zero) */
1093  SphHarmFrequencySeries *hlms = NULL;
1094  XLAL_TRY(hlms=XLALSimInspiralChooseFDModes(m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z, deltaF, f_start, f_max, f_ref, phi0, corrected_distance, 0., grParams, approximant), errnum);
1095 
1096  if(hlms==NULL)
1097  ret=XLAL_FAILURE;
1098  else{
1099  ret=XLAL_SUCCESS;
1100 
1101  /* apply FTA corrections */
1102  COMPLEX16FrequencySeries *hlm_mode = NULL;
1103  SphHarmFrequencySeries *hlms_temp = hlms;
1104  INT4 length;
1105  while(hlms_temp){
1106  /* Extract the mode from the SphericalHarmFrequencySeries */
1107  hlm_mode = XLALSphHarmFrequencySeriesGetMode(hlms_temp, hlms_temp->l,hlms_temp->m);
1108 
1109  /* Resize the modes to keep only positive frequencies. ChooseFDModes returns both positive and negative frequencies */
1110  length = hlm_mode->data->length -1;
1111  hlm_mode = XLALResizeCOMPLEX16FrequencySeries(hlm_mode, (INT4) length/2 + 1, length);
1112 
1113  /* Apply FTA correction */
1114  XLAL_TRY(ret = XLALSimInspiralTestingGRCorrections(hlm_mode,hlms_temp->l,abs(hlms_temp->m),m1*LAL_MSUN_SI,m2*LAL_MSUN_SI, spin1z, spin2z, f_start, f_ref, correction_window, correction_ncycles_taper, model->LALpars), errnum);
1115  hlms_temp = hlms_temp->next;
1116  }
1117 
1118  /* Construct hp and hc from hlms */
1119  size_t npts_wave = hlm_mode->data->length;
1120 
1121  hptilde = XLALCreateCOMPLEX16FrequencySeries("FD hplus",
1122  &(hlm_mode->epoch), 0.0, deltaF,
1123  &(hlm_mode->sampleUnits), npts_wave);
1124  hctilde = XLALCreateCOMPLEX16FrequencySeries("FD hcross",
1125  &(hlm_mode->epoch), 0.0, deltaF,
1126  &(hlm_mode->sampleUnits), npts_wave);
1127  memset(hptilde->data->data, 0, npts_wave * sizeof(COMPLEX16));
1128  memset(hctilde->data->data, 0, npts_wave * sizeof(COMPLEX16));
1129 
1130  hlms_temp = hlms;
1131  while(hlms_temp){
1132  /* Add the modes to hptilde and hctilde */
1133  XLAL_TRY(ret=XLALSimAddModeFD(hptilde, hctilde, hlms_temp->mode, inclination, LAL_PI/2. - phi0, hlms_temp->l, hlms_temp->m, 1), errnum);
1134  hlms_temp = hlms_temp->next;
1135  }
1136 
1137  XLALDestroySphHarmFrequencySeries(hlms);
1138  }
1139  }
1140 
1141 
1142  }
1143  /* if the waveform failed to generate, fill the buffer with zeros
1144  * so that the previous waveform is not left there
1145  */
1146  if(ret!=XLAL_SUCCESS){
1147  memset(model->freqhPlus->data->data,0,sizeof(model->freqhPlus->data->data[0])*model->freqhPlus->data->length);
1148  memset(model->freqhCross->data->data,0,sizeof(model->freqhCross->data->data[0])*model->freqhCross->data->length);
1149  if ( hptilde ) XLALDestroyCOMPLEX16FrequencySeries(hptilde);
1150  if ( hctilde ) XLALDestroyCOMPLEX16FrequencySeries(hctilde);
1151  errnum&=~XLAL_EFUNC; /* Mask out the internal function failure bit */
1152  switch(errnum)
1153  {
1154  case XLAL_EDOM:
1155  /* The waveform was called outside its domain. Return an empty vector but not an error */
1156  XLAL_ERROR_VOID(XLAL_EUSR0);
1157  default:
1158  /* Another error occurred that we can't handle. Propogate upward */
1159  XLALSetErrno(errnum);
1160  XLAL_ERROR_VOID(errnum,"%s: Template generation failed in XLALSimInspiralChooseFDWaveformFromCache:\n\
1161 XLALSimInspiralChooseFDWaveformFromCache(&hptilde, &hctilde, \
1162 %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, \
1163 ,model->LALpars,model->waveformCache)\n",__func__,
1164  phi0, deltaF, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
1165  f_start, f_max, f_ref, distance, inclination);
1166  }
1167  }
1168 
1169  if (hptilde==NULL || hptilde->data==NULL || hptilde->data->data==NULL ) {
1170  XLALPrintError(" ERROR in %s: encountered unallocated 'hptilde'.\n",__func__);
1171  XLAL_ERROR_VOID(XLAL_EFAULT);
1172  }
1173  if (hctilde==NULL || hctilde->data==NULL || hctilde->data->data==NULL ) {
1174  XLALPrintError(" ERROR in %s: encountered unallocated 'hctilde'.\n",__func__);
1175  XLAL_ERROR_VOID(XLAL_EFAULT);
1176  }
1177 
1178  INT4 rem=0;
1179  UINT4 size=hptilde->data->length;
1180  if(size>model->freqhPlus->data->length) size=model->freqhPlus->data->length;
1181  memcpy(model->freqhPlus->data->data,hptilde->data->data,sizeof(hptilde->data->data[0])*size);
1182  if( (rem=(model->freqhPlus->data->length - size)) > 0)
1183  memset(&(model->freqhPlus->data->data[size]),0, rem*sizeof(hptilde->data->data[0]) );
1184 
1185  size=hctilde->data->length;
1186  if(size>model->freqhCross->data->length) size=model->freqhCross->data->length;
1187  memcpy(model->freqhCross->data->data,hctilde->data->data,sizeof(hctilde->data->data[0])*size);
1188  if( (rem=(model->freqhCross->data->length - size)) > 0)
1189  memset(&(model->freqhCross->data->data[size]),0, rem*sizeof(hctilde->data->data[0]) );
1190 
1191  REAL8 instant = model->freqhPlus->epoch.gpsSeconds + 1e-9*model->freqhPlus->epoch.gpsNanoSeconds;
1192  LALInferenceSetVariable(model->params, "time", &instant);
1193 
1194  } else {
1195 
1196  XLAL_TRY(ret=XLALSimInspiralChooseTDWaveformFromCache(&hplus, &hcross, phi0, deltaT,
1197  m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, spin1x, spin1y, spin1z,
1198  spin2x, spin2y, spin2z, f_start, f_ref, distance,
1199  inclination, model->LALpars, approximant,model->waveformCache), errnum);
1200  /* if the waveform failed to generate, fill the buffer with zeros
1201  * so that the previous waveform is not left there
1202  */
1203  if(ret!=XLAL_SUCCESS){
1204  memset(model->freqhPlus->data->data,0,sizeof(model->freqhPlus->data->data[0])*model->freqhPlus->data->length);
1205  memset(model->freqhCross->data->data,0,sizeof(model->freqhCross->data->data[0])*model->freqhCross->data->length);
1206  if ( hptilde ) XLALDestroyCOMPLEX16FrequencySeries(hptilde);
1207  if ( hctilde ) XLALDestroyCOMPLEX16FrequencySeries(hctilde);
1208  if ( hplus) XLALDestroyREAL8TimeSeries(hplus);
1209  if ( hcross ) XLALDestroyREAL8TimeSeries(hcross);
1210  errnum&=~XLAL_EFUNC; /* Mask out the internal function failure bit */
1211  switch(errnum)
1212  {
1213  case XLAL_EDOM:
1214  /* The waveform was called outside its domain. Return an empty vector but not an error */
1215  XLAL_ERROR_VOID(XLAL_EUSR0);
1216  default:
1217  /* Another error occurred that we can't handle. Propogate upward */
1218  XLALSetErrno(errnum);
1219  XLAL_ERROR_VOID(errnum,"%s: Template generation failed in XLALSimInspiralChooseTDWaveformFromCache\n\
1220 XLALSimInspiralChooseTDWaveformFromCache(&hplus, &hcross, \
1221 %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, \
1222 model->LALpars,model->waveformCache)\n",__func__,
1223  phi0, deltaT, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
1224  f_start, f_ref, distance, inclination);
1225  }
1226  }
1227 
1228 
1229  /* The following complicated mess is a result of the following considerations:
1230 
1231  1) The discrete time samples of the template and the timeModel
1232  buffers will not, in general line up.
1233 
1234  2) The likelihood function will timeshift the template in the
1235  frequency domain to align it properly with the desired tc in
1236  each detector (these are different because the detectors
1237  receive the signal at different times). Because this
1238  timeshifting is done in the frequency domain, the effective
1239  time-domain template is periodic. We want to avoid the
1240  possibility of non-zero template samples wrapping around from
1241  the start/end of the buffer, since real templates are not
1242  periodic!
1243 
1244  3) If the template apporaches the ends of the timeModel buffer,
1245  then it should be tapered in the same way as the timeData
1246  (currently 0.4 seconds, hard-coded! Tukey window; see
1247  LALInferenceReadData.c, near line 233) so that template and
1248  signal in the data match. However, as an optimization, we
1249  perform only one tapering and FFT-ing in the likelihood
1250  function; subsequent timeshifts for the different detectors
1251  will cause the tapered regions of the template and data to
1252  become mis-aligned.
1253 
1254  The algorthim we use is the following:
1255 
1256  1) Inject the template to align with the nearest sample in the
1257  timeModel buffer to the desired geocent_end time.
1258 
1259  2) Check whether either the start or the end of the template
1260  overlaps the tapered region, plus a safety buffer corresponding
1261  to a conservative estimate of the largest geocenter <-->
1262  detector timeshift.
1263 
1264  a) If there is no overlap at the start or end of the buffer,
1265  we're done.
1266 
1267  b) If there is an overlap, issue one warning per process
1268  (which can be disabled by setting the LAL debug level) about
1269  a too-short segment length, and return.
1270 */
1271 
1272  size_t waveLength = hplus->data->length;
1273  size_t bufLength = model->timehPlus->data->length;
1274 
1275  /* 2*Rearth/(c*deltaT)---2 is safety factor---is the maximum time
1276  shift for any earth-based detector. */
1277  size_t maxShift = (size_t)lround(4.255e-2/deltaT);
1278 
1279  /* Taper 0.4 seconds at start and end (hard-coded! in
1280  LALInferenceReadData.c, around line 233). */
1281  size_t taperLength = (size_t)lround(0.4/deltaT);
1282 
1283  /* Within unsafeLength of ends of buffer, possible danger of
1284  wrapping and/or tapering interactions. */
1285  size_t unsafeLength = taperLength + maxShift;
1286 
1287  REAL8 desiredTc = *(REAL8 *)LALInferenceGetVariable(model->params, "time");
1288  REAL8 tStart = XLALGPSGetREAL8(&(model->timehPlus->epoch));
1289  REAL8 tEnd = tStart + deltaT * model->timehPlus->data->length;
1290 
1291  if (desiredTc < tStart || desiredTc > tEnd) {
1292  XLALDestroyREAL8TimeSeries(hplus);
1293  XLALDestroyREAL8TimeSeries(hcross);
1294 
1295  XLAL_PRINT_ERROR("desired tc (%.4f) outside data buffer\n", desiredTc);
1296  XLAL_ERROR_VOID(XLAL_EDOM);
1297  }
1298 
1299  /* The nearest sample in model buffer to the desired tc. */
1300  size_t tcSample = (size_t)lround((desiredTc - XLALGPSGetREAL8(&(model->timehPlus->epoch)))/deltaT);
1301 
1302  /* The acutal coalescence time that corresponds to the buffer
1303  sample on which the waveform's tC lands. */
1304  REAL8 injTc = XLALGPSGetREAL8(&(model->timehPlus->epoch)) + tcSample*deltaT;
1305 
1306  /* The sample at which the waveform reaches tc. */
1307  size_t waveTcSample = (size_t)lround(-XLALGPSGetREAL8(&(hplus->epoch))/deltaT);
1308 
1309  /* 1 + (number of samples post-tc in waveform) */
1310  size_t wavePostTc = waveLength - waveTcSample;
1311 
1312  size_t bufStartIndex = (tcSample >= waveTcSample ? tcSample - waveTcSample : 0);
1313  size_t bufEndIndex = (wavePostTc + tcSample <= bufLength ? wavePostTc + tcSample : bufLength);
1314  size_t bufWaveLength = bufEndIndex - bufStartIndex;
1315  size_t waveStartIndex = (tcSample >= waveTcSample ? 0 : waveTcSample - tcSample);
1316 
1317  if (bufStartIndex < unsafeLength || (bufLength - bufEndIndex) <= unsafeLength) {
1318  /* The waveform could be timeshifted into a region where it will
1319  be tapered improperly, or even wrap around from the periodic
1320  timeshift. Issue warning. */
1321  if (!sizeWarning) {
1322  fprintf(stderr, "WARNING: Generated template is too long to guarantee that it will not\n");
1323  fprintf(stderr, "WARNING: (a) lie in a tapered region of the time-domain buffer\n");
1324  fprintf(stderr, "WARNING: (b) wrap periodically when timeshifted in likelihood computation\n");
1325  fprintf(stderr, "WARNING: Either of these may cause differences between the template and the\n");
1326  fprintf(stderr, "WARNING: correct GW waveform in each detector.\n");
1327  fprintf(stderr, "WARNING: Parameter estimation will continue, but you should consider\n");
1328  fprintf(stderr, "WARNING: increasing the data segment length (using the --seglen) option.\n");
1329  sizeWarning = 1;
1330  }
1331  }
1332 
1333  /* Clear model buffers */
1334  memset(model->timehPlus->data->data, 0, sizeof(REAL8)*model->timehPlus->data->length);
1335  memset(model->timehCross->data->data, 0, sizeof(REAL8)*model->timehCross->data->length);
1336 
1337  /* Inject */
1338  memcpy(model->timehPlus->data->data + bufStartIndex,
1339  hplus->data->data + waveStartIndex,
1340  bufWaveLength*sizeof(REAL8));
1341  memcpy(model->timehCross->data->data + bufStartIndex,
1342  hcross->data->data + waveStartIndex,
1343  bufWaveLength*sizeof(REAL8));
1344 
1345  LALInferenceSetVariable(model->params, "time", &injTc);
1346  }
1347 
1348  if ( hplus ) XLALDestroyREAL8TimeSeries(hplus);
1349  if ( hcross ) XLALDestroyREAL8TimeSeries(hcross);
1350  if ( hptilde ) XLALDestroyCOMPLEX16FrequencySeries(hptilde);
1351  if ( hctilde ) XLALDestroyCOMPLEX16FrequencySeries(hctilde);
1352 
1353  return;
1354 }
1355 
1356 
1357 void LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated(LALInferenceModel *model)
1358 /*************************************************************************************************************************/
1359 /* Wrapper for LALSimulation waveforms: */
1360 /* XLALSimInspiralChooseFDWaveform() and XLALSimInspiralChooseTDWaveform(). */
1361 /* */
1362 /* model->params parameters are: */
1363 /* - "name" description; type OPTIONAL (default value) */
1364 /* */
1365 /* MODEL PARAMETERS */
1366 /* - "LAL_APPROXIMANT Approximant; Approximant */
1367 /* - "LAL_PNORDER" Phase PN order; INT4 */
1368 /* - "LAL_AMPORDER" Amplitude PN order; INT4 OPTIONAL (-1) */
1369 /* - "spinO" Spin order; LALSimInspiralSpinOrder OPTIONAL (LAL_SIM_INSPIRAL_SPIN_ORDER_DEFAULT) */
1370 /* - "tideO" Tidal order; LALSimInspiralTidalOrder OPTIONAL (LAL_SIM_INSPIRAL_TIDAL_ORDER_DEFAULT) */
1371 /* - "f_ref" frequency at which the (frequency dependent) parameters are defined; REAL8 OPTIONAL (0.0) */
1372 /* - "fLow" lower frequency bound; REAL8 OPTIONAL (model->fLow) */
1373 /* */
1374 /* MASS PARAMETERS; either: */
1375 /* - "mass1" mass of object 1 in solar mass; REAL8 */
1376 /* - "mass2" mass of object 1 in solar mass; REAL8 */
1377 /* OR */
1378 /* - "chirpmass" chirpmass in solar mass; REAL8 */
1379 /* - "q" asymmetric mass ratio m2/m1, 0<q<1; REAL8 */
1380 /* OR */
1381 /* - "chirpmass" chirpmass in solar mass; REAL8 */
1382 /* - "eta" symmetric mass ratio (m1*m2)/(m1+m2)^2; REAL8 */
1383 /* */
1384 /* ORIENTATION AND SPIN PARAMETERS */
1385 /* - "phi0" reference phase as per LALSimulation convention; REAL8 */
1386 /* - "logdistance" distance in Mpc */
1387 /* - "costheta_jn"); cos of zenith angle between J and N in radians; REAL8 */
1388 /* - "phi_jl"); azimuthal angle of L_N on its cone about J radians; REAL8 */
1389 /* - "tilt_spin1"); zenith angle between S1 and LNhat in radians; REAL8 */
1390 /* - "tilt_spin2"); zenith angle between S2 and LNhat in radians; REAL8 */
1391 /* - "phi12"); difference in azimuthal angle between S1, S2 in radians; REAL8 */
1392 /* - "a_spin1" magnitude of spin 1 in general configuration, -1<a_spin1<1; REAL8 OPTIONAL (0.0) */
1393 /* - "a_spin2" magnitude of spin 2 in general configuration, -1<a_spin1<1; REAL8 OPTIONAL (0.0) */
1394 /* */
1395 /* OTHER PARAMETERS */
1396 /* - "lambda1" tidal parameter of object 1; REAL8 OPTIONAL (0.0) */
1397 /* - "lambda2" tidal parameter of object 1; REAL8 OPTIONAL (0.0) */
1398 /* */
1399 /* - "time" used as an OUTPUT only; REAL8 */
1400 /* */
1401 /* */
1402 /* model needs to also contain: */
1403 /* - model->fLow Unless - "fLow" OPTIONAL */
1404 /* - model->deltaT */
1405 /* - if model->domain == LAL_SIM_DOMAIN_FREQUENCY */
1406 /* - model->deltaF */
1407 /* - model->freqhCross */
1408 /* - model->freqhPlus */
1409 /* - else */
1410 /* - model->timehPlus */
1411 /* - model->timehCross */
1412 /*************************************************************************************************************************/
1413 {
1414  Approximant approximant = (Approximant) 0;
1415 
1416  static int sizeWarning = 0;
1417  int ret=0;
1418  INT4 errnum=0;
1419  REAL8TimeSeries *hplus=NULL; /**< +-polarization waveform [returned] */
1420  REAL8TimeSeries *hcross=NULL; /**< x-polarization waveform [returned] */
1421  COMPLEX16FrequencySeries *hptilde=NULL, *hctilde=NULL;
1422  REAL8 mc;
1423  REAL8 phi0, deltaT, m1, m2, f_low, f_start, distance, inclination;
1424 
1425  REAL8 *m1_p,*m2_p;
1426  REAL8 f_max;
1427 
1428  /* Sampling rate for time domain models */
1429  deltaT = model->deltaT;
1430 
1431  if (LALInferenceCheckVariable(model->params, "LAL_APPROXIMANT"))
1432  approximant = *(Approximant*) LALInferenceGetVariable(model->params, "LAL_APPROXIMANT");
1433  else {
1434  XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_APPROXIMANT\" parameter not provided!\n");
1435  XLAL_ERROR_VOID(XLAL_EDATA);
1436  }
1437 
1438  if (LALInferenceCheckVariable(model->params, "LAL_PNORDER")) {
1439  if (XLALSimInspiralWaveformParamsLookupPNPhaseOrder(model->LALpars)!=(*(INT4*) LALInferenceGetVariable(model->params, "LAL_PNORDER")))
1440  XLALSimInspiralWaveformParamsInsertPNPhaseOrder(model->LALpars, *(INT4*) LALInferenceGetVariable(model->params, "LAL_PNORDER"));
1441  }
1442  else {
1443  XLALPrintError(" ERROR in templateLALGenerateInspiral(): (INT4) \"LAL_PNORDER\" parameter not provided!\n");
1444  XLAL_ERROR_VOID(XLAL_EDATA);
1445  }
1446 
1447  /* Explicitly set the default amplitude order if one is not specified.
1448  * This serves two purposes:
1449  * 1) The default behavior of the code won't change unexpectedly due to changes in LALSimulation.
1450  * 2) We need to know the amplitude order in order to set the starting frequency of the waveform properly. */
1451 
1452  if (LALInferenceCheckVariable(model->params, "LAL_AMPORDER")) {
1453  if ( XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(model->LALpars)!= *(INT4*) LALInferenceGetVariable(model->params, "LAL_AMPORDER"))
1454  XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars, *(INT4*) LALInferenceGetVariable(model->params, "LAL_AMPORDER"));
1455  }
1456  else
1457  if (!XLALSimInspiralWaveformParamsPNAmplitudeOrderIsDefault(model->LALpars))
1458  XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(model->LALpars,-1);
1459 
1460  REAL8 f_ref = 100.0;
1461  if (LALInferenceCheckVariable(model->params, "f_ref")) f_ref = *(REAL8 *)LALInferenceGetVariable(model->params, "f_ref");
1462 
1463  REAL8 fTemp = f_ref;
1464 
1465  if(LALInferenceCheckVariable(model->params,"chirpmass"))
1466  {
1467  mc = *(REAL8*) LALInferenceGetVariable(model->params, "chirpmass");
1468  if (LALInferenceCheckVariable(model->params,"q")) {
1469  REAL8 q = *(REAL8 *)LALInferenceGetVariable(model->params,"q");
1470  q2masses(mc, q, &m1, &m2);
1471  } else {
1472  REAL8 eta = *(REAL8*) LALInferenceGetVariable(model->params, "eta");
1473  mc2masses(mc, eta, &m1, &m2);
1474  }
1475  }
1476  else if((m1_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass1")) && (m2_p=(REAL8 *)LALInferenceGetVariable(model->params, "mass2")))
1477  {
1478  m1=*m1_p;
1479  m2=*m2_p;
1480  }
1481  else
1482  {
1483  fprintf(stderr,"No mass parameters found!");
1484  exit(0);
1485  }
1486 
1487  distance = exp(LALInferenceGetREAL8Variable(model->params, "logdistance"))* LAL_PC_SI * 1.0e6; /* distance (1 Mpc) in units of metres */
1488 
1489  phi0 = LALInferenceGetREAL8Variable(model->params, "phase"); /* START phase as per lalsimulation convention, radians*/
1490 
1491  /* Zenith angle between J and N in radians. Also known as inclination angle when spins are aligned */
1492  REAL8 thetaJN = acos(LALInferenceGetREAL8Variable(model->params, "costheta_jn")); /* zenith angle between J and N in radians */
1493 
1494  /* Check if fLow is a model parameter, otherwise use data structure definition */
1495  if(LALInferenceCheckVariable(model->params, "flow"))
1496  f_low = *(REAL8*) LALInferenceGetVariable(model->params, "flow");
1497  else
1498  f_low = model->fLow;
1499 
1500  f_start = XLALSimInspiralfLow2fStart(f_low, XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(model->LALpars), approximant);
1501 
1502  /* Don't let TaylorF2 generate unphysical inspiral up to Nyquist */
1503  if (approximant == TaylorF2)
1504  f_max = 0.0; /* this will stop at ISCO */
1505  else
1506  f_max = model->fHigh; /* this will be the highest frequency used across the network */
1507 
1508  /* ==== SPINS ==== */
1509  /* We will default to spinless signal and then add in the spin components if required */
1510  /* If there are non-aligned spins, we must convert between the System Frame coordinates
1511  * and the cartestian coordinates */
1512 
1513  /* The cartesian spin coordinates (default 0), as passed to LALSimulation */
1514  REAL8 spin1x = 0.0;
1515  REAL8 spin1y = 0.0;
1516  REAL8 spin1z = 0.0;
1517  REAL8 spin2x = 0.0;
1518  REAL8 spin2y = 0.0;
1519  REAL8 spin2z = 0.0;
1520 
1521  /* System frame coordinates as used for jump proposals */
1522  REAL8 a_spin1 = 0.0; /* Magnitude of spin1 */
1523  REAL8 a_spin2 = 0.0; /* Magnitude of spin2 */
1524  REAL8 phiJL = 0.0; /* azimuthal angle of L_N on its cone about J radians */
1525  REAL8 tilt1 = 0.0; /* zenith angle between S1 and LNhat in radians */
1526  REAL8 tilt2 = 0.0; /* zenith angle between S2 and LNhat in radians */
1527  REAL8 phi12 = 0.0; /* difference in azimuthal angle btwn S1, S2 in radians */
1528 
1529  /* Now check if we have spin amplitudes */
1530  if(LALInferenceCheckVariable(model->params, "a_spin1")) a_spin1 = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin1");
1531  if(LALInferenceCheckVariable(model->params, "a_spin2")) a_spin2 = *(REAL8*) LALInferenceGetVariable(model->params, "a_spin2");
1532 
1533  /* Check if we have spin angles too */
1534  if(LALInferenceCheckVariable(model->params, "phi_jl"))
1535  phiJL = LALInferenceGetREAL8Variable(model->params, "phi_jl");
1536  if(LALInferenceCheckVariable(model->params, "tilt_spin1"))
1537  tilt1 = LALInferenceGetREAL8Variable(model->params, "tilt_spin1");
1538  if(LALInferenceCheckVariable(model->params, "tilt_spin2"))
1539  tilt2 = LALInferenceGetREAL8Variable(model->params, "tilt_spin2");
1540  if(LALInferenceCheckVariable(model->params, "phi12"))
1541  phi12 = LALInferenceGetREAL8Variable(model->params, "phi12");
1542 
1543  /* If we have tilt angles zero, then the spins are aligned and we just set the z component */
1544  /* However, if the waveform supports precession then we still need to get the right coordinate components */
1545  if(tilt1==0.0 && tilt2==0.0)
1546  {
1547  spin1z=a_spin1;
1548  spin2z=a_spin2;
1549  inclination = thetaJN; /* Inclination angle is just thetaJN */
1550  }
1551  else
1552  { /* Template is not aligned-spin only. */
1553  /* Set all the other spin components according to the angles we received above */
1554  /* The transformation function doesn't know fLow, so f_ref==0 isn't interpretted as a request to use the starting frequency for reference. */
1555  if(fTemp==0.0)
1556  fTemp = f_start;
1557  // If we have a *time-domain* AND *precessing* EOB template, set fTemp to f_start.
1558  // This is done since EOB only uses f_start and ignores f_ref.
1559  if(model->domain == LAL_SIM_DOMAIN_TIME && (strstr(XLALSimInspiralGetStringFromApproximant(approximant),"EOB") != NULL)){
1560  fTemp = f_start;
1561  }
1562  XLAL_TRY(ret=XLALSimInspiralTransformPrecessingNewInitialConditions(
1563  &inclination, &spin1x, &spin1y, &spin1z, &spin2x, &spin2y, &spin2z,
1564  thetaJN, phiJL, tilt1, tilt2, phi12, a_spin1, a_spin2, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, fTemp, phi0), errnum);
1565  if (ret == XLAL_FAILURE)
1566  {
1567  XLALPrintError(" ERROR in XLALSimInspiralTransformPrecessingNewInitialConditions(): error converting angles. errnum=%d: %s\n",errnum, XLALErrorString(errnum) );
1568  return;
1569  }
1570  }
1571 /* ==== Spin induced quadrupole moment PARAMETERS ==== */
1572  if(LALInferenceCheckVariable(model->params, "dQuadMonS")&&LALInferenceCheckVariable(model->params, "dQuadMonA")){
1573  REAL8 dQuadMon1=0.;
1574  REAL8 dQuadMon2=0.;
1575  REAL8 dQuadMonS = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonS");
1576  REAL8 dQuadMonA = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMonA");
1577  LALInferencedQuadMonSdQuadMonA(dQuadMonS,dQuadMonA,&dQuadMon1,&dQuadMon2);
1578  XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
1579  XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
1580  }
1581  else if(LALInferenceCheckVariable(model->params, "dQuadMon1")&&LALInferenceCheckVariable(model->params, "dQuadMon2")){
1582  REAL8 dQuadMon1 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon1");
1583  REAL8 dQuadMon2 = *(REAL8*) LALInferenceGetVariable(model->params, "dQuadMon2");
1584  XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars,dQuadMon1);
1585  XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars,dQuadMon2);
1586  }
1587  /* ==== TIDAL PARAMETERS ==== */
1588  if(LALInferenceCheckVariable(model->params, "lambda1"))
1589  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda1"));
1590  if(LALInferenceCheckVariable(model->params, "lambda2"))
1591  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, *(REAL8*) LALInferenceGetVariable(model->params, "lambda2"));
1592  REAL8 lambdaT = 0.;
1593  REAL8 dLambdaT = 0.;
1594  REAL8 sym_mass_ratio_eta = 0.;
1595  if(LALInferenceCheckVariable(model->params, "lambdaT")&&LALInferenceCheckVariable(model->params, "dLambdaT")){
1596  REAL8 lambda1=0.;
1597  REAL8 lambda2=0.;
1598  lambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "lambdaT");
1599  dLambdaT = *(REAL8*) LALInferenceGetVariable(model->params, "dLambdaT");
1600  sym_mass_ratio_eta = m1*m2/((m1+m2)*(m1+m2));
1601  LALInferenceLambdaTsEta2Lambdas(lambdaT,dLambdaT,sym_mass_ratio_eta,&lambda1,&lambda2);
1602  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
1603  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
1604  }
1605  /* ==== BINARY_LOVE PARAMETERS ==== */
1606  if(LALInferenceCheckVariable(model->params, "lambdaS")){
1607  REAL8 lambda1=0.;
1608  REAL8 lambda2=0.;
1609  LALInferenceBinaryLove(model->params, &lambda1, &lambda2);
1610  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
1611  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
1612  LALInferenceAddVariable(model->params, "lambda1", &lambda1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
1613  LALInferenceAddVariable(model->params, "lambda2", &lambda2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
1614  }
1615 
1616  /* ==== 4-PIECE POLYTROPE EOS PARAMETERS ==== */
1617  REAL8 logp1 = 0.;
1618  REAL8 gamma1 = 0.;
1619  REAL8 gamma2 = 0.;
1620  REAL8 gamma3 = 0.;
1621  if(LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "gamma1")&&LALInferenceCheckVariable(model->params, "gamma2")&&LALInferenceCheckVariable(model->params, "gamma3")){
1622  REAL8 lambda1 = 0.;
1623  REAL8 lambda2 = 0.;
1624  logp1 = *(REAL8*) LALInferenceGetVariable(model->params, "logp1");
1625  gamma1 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma1");
1626  gamma2 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma2");
1627  gamma3 = *(REAL8*) LALInferenceGetVariable(model->params, "gamma3");
1628  // Find lambda1,2(m1,2|eos)
1629  LALInferenceLogp1GammasMasses2Lambdas(logp1,gamma1,gamma2,gamma3,m1,m2,&lambda1,&lambda2);
1630  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
1631  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
1632  }
1633 
1634  /* ==== SPECTRAL DECOMPOSITION PARAMETERS ==== */
1635  REAL8 SDgamma0 = 0.;
1636  REAL8 SDgamma1 = 0.;
1637  REAL8 SDgamma2 = 0.;
1638  REAL8 SDgamma3 = 0.;
1639  if(!LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "SDgamma0")&&LALInferenceCheckVariable(model->params, "SDgamma1")&&LALInferenceCheckVariable(model->params, "SDgamma2")&&LALInferenceCheckVariable(model->params,"SDgamma3")){
1640  REAL8 lambda1 = 0.;
1641  REAL8 lambda2 = 0.;
1642  SDgamma0 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma0");
1643  SDgamma1 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma1");
1644  SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
1645  SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
1646  REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
1647  LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
1648  XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
1649  XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
1650  }
1651 
1652 
1653  /* Only use GR templates */
1654  /* Fill in the extra parameters for testing GR, if necessary */
1655  for (UINT4 k=0; k<N_extra_params; k++)
1656  {
1657  if(LALInferenceCheckVariable(model->params,list_extra_parameters[k]))
1658  {
1659  XLALDictInsert(model->LALpars, list_extra_parameters[k], (void *)LALInferenceGetVariable(model->params,list_extra_parameters[k]), sizeof(double), LAL_D_TYPE_CODE);
1660  }
1661  }
1662  /* Fill in PPE params if they are available */
1663  char PPEparam[64]="";
1664  const char *PPEnames[]= {"aPPE","alphaPPE","bPPE","betaPPE",NULL};
1665  for(UINT4 idx=0; PPEnames[idx]; idx++)
1666  {
1667  for(UINT4 ppeidx=0;; ppeidx++)
1668  {
1669  sprintf(PPEparam, "%s%d",PPEnames[idx],ppeidx);
1670  if(LALInferenceCheckVariable(model->params,PPEparam))
1671  XLALDictInsert(model->LALpars, PPEparam, (void *)LALInferenceGetVariable(model->params,PPEparam), sizeof(double), LAL_D_TYPE_CODE);
1672  else
1673  break;
1674  }
1675  }
1676  INT4 Nbands=-1; /* Use optimum number of bands */
1677  double mc_min=1.0/pow(2,0.2); /* For min 1.0-1.0 waveform */
1678 
1679  /* Vector of frequencies at which to compute FD template */
1680  static REAL8Sequence *frequencies = NULL;
1681 
1682  /* ==== Call the waveform generator ==== */
1683  if(model->domain == LAL_SIM_DOMAIN_FREQUENCY) {
1684  if(!frequencies) frequencies = LALInferenceMultibandFrequencies(Nbands,f_start,0.5/deltaT, model->deltaF, mc_min);
1685  double corrected_distance = distance * sqrt(model->window->sumofsquares/model->window->data->length);
1686 
1687 
1688  XLAL_TRY(ret=XLALSimInspiralChooseFDWaveformFromCache(&hptilde, &hctilde, phi0,
1689  0.0, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, spin1x, spin1y, spin1z,
1690  spin2x, spin2y, spin2z, f_start, f_max, f_ref, corrected_distance, inclination, model->LALpars,
1691  approximant,model->waveformCache, frequencies), errnum);
1692 
1693  /* if the waveform failed to generate, fill the buffer with zeros
1694  * so that the previous waveform is not left there
1695  */
1696  if(ret!=XLAL_SUCCESS){
1697  memset(model->freqhPlus->data->data,0,sizeof(model->freqhPlus->data->data[0])*model->freqhPlus->data->length);
1698  memset(model->freqhCross->data->data,0,sizeof(model->freqhCross->data->data[0])*model->freqhCross->data->length);
1699  if ( hptilde ) XLALDestroyCOMPLEX16FrequencySeries(hptilde);
1700  if ( hctilde ) XLALDestroyCOMPLEX16FrequencySeries(hctilde);
1701  errnum&=~XLAL_EFUNC; /* Mask out the internal function failure bit */
1702  switch(errnum)
1703  {
1704  case XLAL_EDOM:
1705  /* The waveform was called outside its domain. Return an empty vector but not an error */
1706  XLAL_ERROR_VOID(XLAL_EUSR0);
1707  default:
1708  /* Another error occurred that we can't handle. Propogate upward */
1709  XLALSetErrno(errnum);
1710  XLAL_ERROR_VOID(errnum,"%s: Template generation failed in XLALSimInspiralChooseFDWaveformFromCache:\n\
1711 XLALSimInspiralChooseFDWaveformFromCache(&hptilde, &hctilde, \
1712 %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, \
1713 model->LALpars,%d,model->waveformCache)\n",__func__,
1714  phi0, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
1715  f_start, f_max, f_ref, distance, inclination, approximant);
1716  }
1717  }
1718 
1719  if (hptilde==NULL || hptilde->data==NULL || hptilde->data->data==NULL ) {
1720  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimInspiralChooseWaveform(): encountered unallocated 'hptilde'.\n");
1721  XLAL_ERROR_VOID(XLAL_EFAULT);
1722  }
1723  if (hctilde==NULL || hctilde->data==NULL || hctilde->data->data==NULL ) {
1724  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimInspiralChooseWaveform(): encountered unallocated 'hctilde'.\n");
1725  XLAL_ERROR_VOID(XLAL_EFAULT);
1726  }
1727 
1728 
1729  InterpolateWaveform(frequencies, hptilde, model->freqhPlus);
1730  InterpolateWaveform(frequencies, hctilde, model->freqhCross);
1731 
1732  REAL8 instant = model->freqhPlus->epoch.gpsSeconds + 1e-9*model->freqhPlus->epoch.gpsNanoSeconds;
1733  LALInferenceSetVariable(model->params, "time", &instant);
1734 
1735 
1736  } else {
1737 
1738  XLAL_TRY(ret=XLALSimInspiralChooseTDWaveformFromCache(&hplus, &hcross, phi0, deltaT,
1739  m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, spin1x, spin1y, spin1z,
1740  spin2x, spin2y, spin2z, f_start, f_ref, distance,
1741  inclination, model->LALpars,
1742  approximant,model->waveformCache), errnum);
1743  /* if the waveform failed to generate, fill the buffer with zeros
1744  * so that the previous waveform is not left there
1745  */
1746  if(ret!=XLAL_SUCCESS){
1747  memset(model->freqhPlus->data->data,0,sizeof(model->freqhPlus->data->data[0])*model->freqhPlus->data->length);
1748  memset(model->freqhCross->data->data,0,sizeof(model->freqhCross->data->data[0])*model->freqhCross->data->length);
1749  if ( hptilde ) XLALDestroyCOMPLEX16FrequencySeries(hptilde);
1750  if ( hctilde ) XLALDestroyCOMPLEX16FrequencySeries(hctilde);
1751  errnum&=~XLAL_EFUNC; /* Mask out the internal function failure bit */
1752  switch(errnum)
1753  {
1754  case XLAL_EDOM:
1755  /* The waveform was called outside its domain. Return an empty vector but not an error */
1756  XLAL_ERROR_VOID(XLAL_EUSR0);
1757  default:
1758  /* Another error occurred that we can't handle. Propogate upward */
1759  XLALSetErrno(errnum);
1760  XLAL_ERROR_VOID(errnum,"%s: Template generation failed in XLALSimInspiralChooseTDWaveformFromCache",__func__);
1761  }
1762  }
1763 
1764  /* The following complicated mess is a result of the following considerations:
1765 
1766  1) The discrete time samples of the template and the timeModel
1767  buffers will not, in general line up.
1768 
1769  2) The likelihood function will timeshift the template in the
1770  frequency domain to align it properly with the desired tc in
1771  each detector (these are different because the detectors
1772  receive the signal at different times). Because this
1773  timeshifting is done in the frequency domain, the effective
1774  time-domain template is periodic. We want to avoid the
1775  possibility of non-zero template samples wrapping around from
1776  the start/end of the buffer, since real templates are not
1777  periodic!
1778 
1779  3) If the template apporaches the ends of the timeModel buffer,
1780  then it should be tapered in the same way as the timeData
1781  (currently 0.4 seconds, hard-coded! Tukey window; see
1782  LALInferenceReadData.c, near line 233) so that template and
1783  signal in the data match. However, as an optimization, we
1784  perform only one tapering and FFT-ing in the likelihood
1785  function; subsequent timeshifts for the different detectors
1786  will cause the tapered regions of the template and data to
1787  become mis-aligned.
1788 
1789  The algorthim we use is the following:
1790 
1791  1) Inject the template to align with the nearest sample in the
1792  timeModel buffer to the desired geocent_end time.
1793 
1794  2) Check whether either the start or the end of the template
1795  overlaps the tapered region, plus a safety buffer corresponding
1796  to a conservative estimate of the largest geocenter <-->
1797  detector timeshift.
1798 
1799  a) If there is no overlap at the start or end of the buffer,
1800  we're done.
1801 
1802  b) If there is an overlap, issue one warning per process
1803  (which can be disabled by setting the LAL debug level) about
1804  a too-short segment length, and return.
1805  */
1806 
1807  size_t waveLength = hplus->data->length;
1808  size_t bufLength = model->timehPlus->data->length;
1809 
1810  /* 2*Rearth/(c*deltaT)---2 is safety factor---is the maximum time
1811  shift for any earth-based detector. */
1812  size_t maxShift = (size_t)lround(4.255e-2/deltaT);
1813 
1814  /* Taper 0.4 seconds at start and end (hard-coded! in
1815  LALInferenceReadData.c, around line 233). */
1816  size_t taperLength = (size_t)lround(0.4/deltaT);
1817 
1818  /* Within unsafeLength of ends of buffer, possible danger of
1819  wrapping and/or tapering interactions. */
1820  size_t unsafeLength = taperLength + maxShift;
1821 
1822  REAL8 desiredTc = *(REAL8 *)LALInferenceGetVariable(model->params, "time");
1823  REAL8 tStart = XLALGPSGetREAL8(&(model->timehPlus->epoch));
1824  REAL8 tEnd = tStart + deltaT * model->timehPlus->data->length;
1825 
1826  if (desiredTc < tStart || desiredTc > tEnd) {
1827  XLALDestroyREAL8TimeSeries(hplus);
1828  XLALDestroyREAL8TimeSeries(hcross);
1829 
1830  XLAL_PRINT_ERROR("desired tc (%.4f) outside data buffer\n", desiredTc);
1831  XLAL_ERROR_VOID(XLAL_EDOM);
1832  }
1833 
1834  /* The nearest sample in model buffer to the desired tc. */
1835  size_t tcSample = (size_t)lround((desiredTc - XLALGPSGetREAL8(&(model->timehPlus->epoch)))/deltaT);
1836 
1837  /* The acutal coalescence time that corresponds to the buffer
1838  sample on which the waveform's tC lands. */
1839  REAL8 injTc = XLALGPSGetREAL8(&(model->timehPlus->epoch)) + tcSample*deltaT;
1840 
1841  /* The sample at which the waveform reaches tc. */
1842  size_t waveTcSample = (size_t)lround(-XLALGPSGetREAL8(&(hplus->epoch))/deltaT);
1843 
1844  /* 1 + (number of samples post-tc in waveform) */
1845  size_t wavePostTc = waveLength - waveTcSample;
1846 
1847  size_t bufStartIndex = (tcSample >= waveTcSample ? tcSample - waveTcSample : 0);
1848  size_t bufEndIndex = (wavePostTc + tcSample <= bufLength ? wavePostTc + tcSample : bufLength);
1849  size_t bufWaveLength = bufEndIndex - bufStartIndex;
1850  size_t waveStartIndex = (tcSample >= waveTcSample ? 0 : waveTcSample - tcSample);
1851 
1852  if (bufStartIndex < unsafeLength || (bufLength - bufEndIndex) <= unsafeLength) {
1853  /* The waveform could be timeshifted into a region where it will
1854  be tapered improperly, or even wrap around from the periodic
1855  timeshift. Issue warning. */
1856  if (!sizeWarning) {
1857  fprintf(stderr, "WARNING: Generated template is too long to guarantee that it will not\n");
1858  fprintf(stderr, "WARNING: (a) lie in a tapered region of the time-domain buffer\n");
1859  fprintf(stderr, "WARNING: (b) wrap periodically when timeshifted in likelihood computation\n");
1860  fprintf(stderr, "WARNING: Either of these may cause differences between the template and the\n");
1861  fprintf(stderr, "WARNING: correct GW waveform in each detector.\n");
1862  fprintf(stderr, "WARNING: Parameter estimation will continue, but you should consider\n");
1863  fprintf(stderr, "WARNING: increasing the data segment length (using the --seglen) option.\n");
1864  sizeWarning = 1;
1865  }
1866  }
1867 
1868  /* Clear model buffers */
1869  memset(model->timehPlus->data->data, 0, sizeof(REAL8)*model->timehPlus->data->length);
1870  memset(model->timehCross->data->data, 0, sizeof(REAL8)*model->timehCross->data->length);
1871 
1872  /* Inject */
1873  memcpy(model->timehPlus->data->data + bufStartIndex,
1874  hplus->data->data + waveStartIndex,
1875  bufWaveLength*sizeof(REAL8));
1876  memcpy(model->timehCross->data->data + bufStartIndex,
1877  hcross->data->data + waveStartIndex,
1878  bufWaveLength*sizeof(REAL8));
1879 
1880  LALInferenceSetVariable(model->params, "time", &injTc);
1881  }
1882 
1883  if ( hplus ) XLALDestroyREAL8TimeSeries(hplus);
1884  if ( hcross ) XLALDestroyREAL8TimeSeries(hcross);
1885  if ( hptilde ) XLALDestroyCOMPLEX16FrequencySeries(hptilde);
1886  if ( hctilde ) XLALDestroyCOMPLEX16FrequencySeries(hctilde);
1887 
1888  return;
1889 }
1890 
1891 void LALInferenceTemplateXLALSimBurstChooseWaveform(LALInferenceModel *model)
1892 /*************************************************************************************************************************/
1893 /* Wrapper for LALSimulation waveforms: */
1894 /* XLALSimBurstChooseFDWaveform() and XLALSimBurstChooseTDWaveform(). */
1895 /* */
1896 /* model->params parameters are: */
1897 /* - "name" description; type OPTIONAL (default value) */
1898 /* "LAL_APPROXIMANT" burst approximant, BurstApproximant */
1899 /* "frequency" central frequency, REAL8 */
1900 /* "Q" quality, REAL8 (optional, depending on the WF) */
1901 /* "duration" duration, REAL8 (optional, depending on the WF) */
1902 /* "polar_angle" ellipticity polar angle, REAL8 */
1903 /* "polar_eccentricity" ellipticity ellipse eccentricity, REAL8 (optional) */
1904 /* */
1905 /*************************************************************************************************************************/
1906 {
1907 
1908  BurstApproximant approximant = (BurstApproximant) 0;
1909  unsigned long i;
1910  static int sizeWarning = 0;
1911  int ret=0;
1912  INT4 errnum=0;
1913  REAL8 instant;
1914 
1915 
1916  REAL8TimeSeries *hplus=NULL; /**< +-polarization waveform [returned] */
1917  REAL8TimeSeries *hcross=NULL; /**< x-polarization waveform [returned] */
1918  COMPLEX16FrequencySeries *hptilde=NULL, *hctilde=NULL;
1919  REAL8 deltaT,deltaF,
1920  freq=0.0,
1921  quality=0.0,
1922  duration=0.0,
1923  f_low, f_max,
1924  hrss=1.0,
1925  polar_ecc=1.0,polar_angle=LAL_PI/2.;
1926  LALSimBurstExtraParam *extraParams = NULL;
1927 
1928  if (LALInferenceCheckVariable(model->params, "LAL_APPROXIMANT"))
1929  approximant = *(BurstApproximant*) LALInferenceGetVariable(model->params, "LAL_APPROXIMANT");
1930  else {
1931  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimBurstChooseWaveform(): (INT4) \"LAL_APPROXIMANT\" parameter not provided!\n");
1932  XLAL_ERROR_VOID(XLAL_EDATA);
1933  }
1934 
1935  if(LALInferenceCheckVariable(model->params,"frequency"))
1936  {
1937  freq=*(REAL8*) LALInferenceGetVariable(model->params, "frequency");
1938  }
1939 
1940  if(LALInferenceCheckVariable(model->params,"quality"))
1941  {
1942  quality=*(REAL8*) LALInferenceGetVariable(model->params, "quality");
1943  }
1944  if(LALInferenceCheckVariable(model->params,"duration"))
1945  {
1946  duration=*(REAL8*) LALInferenceGetVariable(model->params, "duration");
1947  }
1948  /*
1949  * not needed since template is calculated at hrss=1 and scaled back in the likelihood
1950  if(LALInferenceCheckVariable(model->params,"hrss"))
1951  {
1952  hrss=*(REAL8*) LALInferenceGetVariable(model->params, "hrss");
1953  } else if (LALInferenceCheckVariable(model->params,"loghrss"))
1954  hrss=exp(*(REAL8*) LALInferenceGetVariable(model->params, "hrss"));
1955 
1956  if(LALInferenceCheckVariable(model->params,"alpha"))
1957  {
1958  alpha=*(REAL8*) LALInferenceGetVariable(model->params, "alpha");
1959  if (!extraParams) extraParams=XLALSimBurstCreateExtraParam("alpha",alpha);
1960  else XLALSimBurstAddExtraParam(&extraParams,"alpha",alpha);
1961  polar_angle=alpha;
1962  }
1963  */
1964  /* If someone wants to use old parametrization, allow for */
1965  if(LALInferenceCheckVariable(model->params,"polar_angle"))
1966  polar_angle=*(REAL8*) LALInferenceGetVariable(model->params, "polar_angle");
1967  if(LALInferenceCheckVariable(model->params,"polar_eccentricity"))
1968  polar_ecc=*(REAL8*) LALInferenceGetVariable(model->params, "polar_eccentricity");
1969 
1970  f_low=0.0; // These are not really used for burst signals.
1971  f_max = model->fHigh; /* this will be the highest frequency used across the network */
1972 
1973 
1974  if (model->timehCross==NULL) {
1975  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimInspiralChooseWaveform(): encountered unallocated 'timeData'.\n");
1976  XLAL_ERROR_VOID(XLAL_EFAULT);
1977  }
1978  deltaT = model->timehCross->deltaT;
1979 
1980  if(model->domain == LAL_SIM_DOMAIN_FREQUENCY) {
1981  if (model->freqhCross==NULL) {
1982  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimInspiralChooseWaveform(): encountered unallocated 'freqhCross'.\n");
1983  XLAL_ERROR_VOID(XLAL_EFAULT);
1984  }
1985 
1986  deltaF = model->deltaF;
1987 
1988  /*Create BurstExtra params here and set depending on approx or let chooseFD do that*/
1989 
1990  /* Correct hrss to account for bias introduced by window normalisation */
1991  double corrected_hrss = hrss / sqrt(model->window->sumofsquares/model->window->data->length);
1992 
1993  XLAL_TRY(ret=XLALSimBurstChooseFDWaveformFromCache(&hptilde, &hctilde, deltaF,deltaT,freq,quality,duration,f_low,f_max,corrected_hrss,polar_angle,polar_ecc,extraParams,approximant,model->burstWaveformCache), errnum);
1994  //XLAL_TRY(ret=XLALSimBurstChooseFDWaveform(&hptilde, &hctilde, deltaF,deltaT,freq,quality,duration,f_low,f_max,hrss,polar_angle,polar_ecc,extraParams,approximant), errnum);
1995  if (ret == XLAL_FAILURE)
1996  {
1997  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimBurstChooseWaveform(). errnum=%d: %s\n",errnum ,XLALErrorString(errnum));
1998  return;
1999  }
2000  if (hptilde==NULL || hptilde->data==NULL || hptilde->data->data==NULL ) {
2001  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimBurstChooseWaveform: encountered unallocated 'hptilde'.\n");
2002  XLAL_ERROR_VOID(XLAL_EFAULT);
2003  }
2004  if (hctilde==NULL || hctilde->data==NULL || hctilde->data->data==NULL ) {
2005  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimBurstChooseWaveform: encountered unallocated 'hctilde'.\n");
2006  XLAL_ERROR_VOID(XLAL_EFAULT);
2007  }
2008 
2009  COMPLEX16 *dataPtr = hptilde->data->data;
2010 
2011  for (i=0; i<model->freqhCross->data->length; ++i) {
2012  dataPtr = hptilde->data->data;
2013  if(i < hptilde->data->length){
2014  model->freqhPlus->data->data[i] = dataPtr[i];
2015  }else{
2016  model->freqhPlus->data->data[i] = 0.0;
2017  }
2018  }
2019  for (i=0; i<model->freqhCross->data->length; ++i) {
2020  dataPtr = hctilde->data->data;
2021  if(i < hctilde->data->length){
2022  model->freqhCross->data->data[i] = dataPtr[i];
2023  }else{
2024  model->freqhCross->data->data[i] = 0.0;
2025  }
2026  }
2027 
2028  /* Destroy the extra params */
2029  XLALSimBurstDestroyExtraParam(extraParams);
2030 
2031  instant= (model->timehCross->epoch.gpsSeconds + 1e-9*model->timehCross->epoch.gpsNanoSeconds);
2032  LALInferenceSetVariable(model->params, "time", &instant);
2033  }
2034  else {
2035  /*Time domain WF*/
2036 
2037  // XLAL_TRY(ret=XLALSimBurstChooseTDWaveform(&hplus, &hcross,deltaT,freq,quality,duration,f_low,f_max,hrss,polar_angle,polar_ecc,extraParams,approximant), errnum);
2038  XLAL_TRY(ret=XLALSimBurstChooseTDWaveformFromCache(&hplus, &hcross,deltaT,freq,quality,duration,f_low,f_max,hrss,polar_angle,polar_ecc,extraParams,approximant,model->burstWaveformCache), errnum);
2039  XLALSimBurstDestroyExtraParam(extraParams);
2040  if (ret == XLAL_FAILURE || hplus == NULL || hcross == NULL)
2041  {
2042  XLALPrintError(" ERROR in XLALSimBurstChooseTDWaveform(): error generating waveform. errnum=%d: %s\n",errnum ,XLALErrorString(errnum));
2043  for (i=0; i<model->timehCross->data->length; i++){
2044  model->timehPlus->data->data[i] = 0.0;
2045  model->timehCross->data->data[i] = 0.0;
2046  }
2047  return;
2048  }
2049 
2050  /* The following complicated mess is a result of the following considerations:
2051 
2052  1) The discrete time samples of the template and the timeModel
2053  buffers will not, in general line up.
2054 
2055  2) The likelihood function will timeshift the template in the
2056  frequency domain to align it properly with the desired tc in
2057  each detector (these are different because the detectors
2058  receive the signal at different times). Because this
2059  timeshifting is done in the frequency domain, the effective
2060  time-domain template is periodic. We want to avoid the
2061  possibility of non-zero template samples wrapping around from
2062  the start/end of the buffer, since real templates are not
2063  periodic!
2064 
2065  3) If the template apporaches the ends of the timeModel buffer,
2066  then it should be tapered in the same way as the timeData
2067  (currently 0.4 seconds, hard-coded! Tukey window; see
2068  LALInferenceReadData.c, near line 233) so that template and
2069  signal in the data match. However, as an optimization, we
2070  perform only one tapering and FFT-ing in the likelihood
2071  function; subsequent timeshifts for the different detectors
2072  will cause the tapered regions of the template and data to
2073  become mis-aligned.
2074 
2075  The algorthim we use is the following:
2076 
2077  1) Inject the template to align with the nearest sample in the
2078  timeModel buffer to the desired geocent_end time.
2079 
2080  2) Check whether either the start or the end of the template
2081  overlaps the tapered region, plus a safety buffer corresponding
2082  to a conservative estimate of the largest geocenter <-->
2083  detector timeshift.
2084 
2085  a) If there is no overlap at the start or end of the buffer,
2086  we're done.
2087 
2088  b) If there is an overlap, issue one warning per process
2089  (which can be disabled by setting the LAL debug level) about
2090  a too-short segment length, and return.
2091 */
2092 
2093  size_t waveLength = hplus->data->length;
2094  size_t bufLength = model->timehCross->data->length;
2095 
2096  /* 2*Rearth/(c*deltaT)---2 is safety factor---is the maximum time
2097  shift for any earth-based detector. */
2098  size_t maxShift = (size_t)lround(4.255e-2/hplus->deltaT);
2099 
2100  /* Taper 0.4 seconds at start and end (hard-coded! in
2101  LALInferenceReadData.c, around line 233). */
2102  REAL8 pad=model->padding;
2103  size_t taperLength = (size_t)lround(pad/hplus->deltaT);
2104 
2105  /* Within unsafeLength of ends of buffer, possible danger of
2106  wrapping and/or tapering interactions. */
2107  size_t unsafeLength = taperLength + maxShift;
2108 
2109  REAL8 desiredTc = *(REAL8 *)LALInferenceGetVariable(model->params, "time");
2110  REAL8 tStart = XLALGPSGetREAL8(&(model->timehCross->epoch));
2111  REAL8 tEnd = tStart + model->deltaT * model->timehCross->data->length;
2112 
2113  if (desiredTc < tStart || desiredTc > tEnd) {
2114  XLALDestroyREAL8TimeSeries(hplus);
2115  XLALDestroyREAL8TimeSeries(hcross);
2116 
2117  XLAL_PRINT_ERROR("desired tc (%.4f) outside data buffer\n", desiredTc);
2118  XLAL_ERROR_VOID(XLAL_EDOM);
2119  }
2120 
2121  /* The nearest sample in model buffer to the desired tc. */
2122  size_t tcSample = (size_t)lround((desiredTc - XLALGPSGetREAL8(&(model->timehPlus->epoch)))/model->deltaT);
2123 
2124  /* The acutal coalescence time that corresponds to the buffer
2125  sample on which the waveform's tC lands. */
2126  REAL8 injTc = XLALGPSGetREAL8(&(model->timehPlus->epoch)) + tcSample*model->deltaT;
2127 
2128  /* The sample at which the waveform reaches tc. */
2129  size_t waveTcSample = (size_t)lround(-XLALGPSGetREAL8(&(hplus->epoch))/hplus->deltaT);
2130 
2131  /* 1 + (number of samples post-tc in waveform) */
2132  size_t wavePostTc = waveLength - waveTcSample;
2133 
2134  size_t bufStartIndex = (tcSample >= waveTcSample ? tcSample - waveTcSample : 0);
2135  size_t bufEndIndex = (wavePostTc + tcSample <= bufLength ? wavePostTc + tcSample : bufLength);
2136  size_t bufWaveLength = bufEndIndex - bufStartIndex;
2137  size_t waveStartIndex = (tcSample >= waveTcSample ? 0 : waveTcSample - tcSample);
2138 
2139  if (bufStartIndex < unsafeLength || (bufLength - bufEndIndex) <= unsafeLength) {
2140  /* The waveform could be timeshifted into a region where it will
2141  be tapered improperly, or even wrap around from the periodic
2142  timeshift. Issue warning. */
2143  if (!sizeWarning) {
2144  fprintf(stderr, "WARNING: Generated template is too long to guarantee that it will not\n");
2145  fprintf(stderr, "WARNING: (a) lie in a tapered region of the time-domain buffer\n");
2146  fprintf(stderr, "WARNING: (b) wrap periodically when timeshifted in likelihood computation\n");
2147  fprintf(stderr, "WARNING: Either of these may cause differences between the template and the\n");
2148  fprintf(stderr, "WARNING: correct GW waveform in each detector.\n");
2149  fprintf(stderr, "WARNING: Parameter estimation will continue, but you should consider\n");
2150  fprintf(stderr, "WARNING: increasing the data segment length (using the --seglen) option.\n");
2151  sizeWarning = 1;
2152 
2153  }
2154  }
2155 
2156  /* Clear IFOdata buffers */
2157  memset(model->timehPlus->data->data, 0, sizeof(REAL8)*model->timehPlus->data->length);
2158  memset(model->timehCross->data->data, 0, sizeof(REAL8)*model->timehCross->data->length);
2159 
2160  /* Inject */
2161  memcpy(model->timehPlus->data->data + bufStartIndex,
2162  hplus->data->data + waveStartIndex,
2163  bufWaveLength*sizeof(REAL8));
2164  memcpy(model->timehCross->data->data + bufStartIndex,
2165  hcross->data->data + waveStartIndex,
2166  bufWaveLength*sizeof(REAL8));
2167 
2168  LALInferenceSetVariable(model->params, "time", &injTc);
2169  }
2170 
2171 
2172  if ( hplus ) XLALDestroyREAL8TimeSeries(hplus);
2173  if ( hcross ) XLALDestroyREAL8TimeSeries(hcross);
2174  if ( hptilde ) XLALDestroyCOMPLEX16FrequencySeries(hptilde);
2175  if ( hctilde ) XLALDestroyCOMPLEX16FrequencySeries(hctilde);
2176 
2177  return;
2178 }
2179 
2180 
2181 void LALInferenceTemplateXLALSimBurstSineGaussianF(LALInferenceModel *model)
2182 /*************************************************************************************************************************/
2183 /* Wrapper for LALSimulation waveforms: */
2184 /* XLALInferenceBurstChooseFDWaveform() and XLALInferenceBurstChooseTDWaveform(). */
2185 /* */
2186 /* model->params parameters are: */
2187 /* - "name" description; type OPTIONAL (default value) */
2188 /* "LAL_APPROXIMANT" burst approximant, BurstApproximant */
2189 /* "frequency" central frequency, REAL8 */
2190 /* "Q" quality, REAL8 (optional, depending on the WF) */
2191 /* "duration" duration, REAL8 (optional, depending on the WF) */
2192 /* "polar_angle" ellipticity polar angle, REAL8 */
2193 /* "polar_eccentricity" ellipticity ellipse eccentricity, REAL8 (optional) */
2194 /* */
2195 /*************************************************************************************************************************/
2196 {
2197 
2198  unsigned long i;
2199  int ret=0;
2200  INT4 errnum=0;
2201  REAL8 instant;
2202 
2203 
2204  REAL8TimeSeries *hplus=NULL; /**< +-polarization waveform [returned] */
2205  REAL8TimeSeries *hcross=NULL; /**< x-polarization waveform [returned] */
2206  COMPLEX16FrequencySeries *hptilde=NULL, *hctilde=NULL;
2207  REAL8 deltaT,deltaF,
2208  freq=0.0,
2209  quality=0.0,tau=0.0,
2210  hrss=1.0;
2211 
2212  freq=*(REAL8*) LALInferenceGetVariable(model->params, "frequency");
2213  quality=*(REAL8*) LALInferenceGetVariable(model->params, "quality");
2214  if(LALInferenceCheckVariable(model->params,"duration"))
2215  {
2216  tau=*(REAL8*) LALInferenceGetVariable(model->params, "duration");
2217  quality=tau*freq*LAL_SQRT2*LAL_PI;
2218  }
2219  REAL8 polar_angle=*(REAL8*) LALInferenceGetVariable(model->params, "polar_angle");
2220  /* If someone wants to use old parametrization, allow for */
2221  REAL8 polar_ecc=*(REAL8*) LALInferenceGetVariable(model->params, "polar_eccentricity");
2222 
2223  if (model->timehCross==NULL) {
2224  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimInspiralChooseWaveform(): encountered unallocated 'timeData'.\n");
2225  XLAL_ERROR_VOID(XLAL_EFAULT);
2226  }
2227  deltaT = model->timehCross->deltaT;
2228 
2229  if (model->freqhCross==NULL) {
2230  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimInspiralChooseWaveform(): encountered unallocated 'freqhCross'.\n");
2231  XLAL_ERROR_VOID(XLAL_EFAULT);
2232  }
2233 
2234  deltaF = model->deltaF;
2235  /* Compute corrected hrss to account for bias introduced by window normalisation */
2236  double corrected_hrss = hrss / sqrt(model->window->sumofsquares/model->window->data->length);
2237 
2238  XLAL_TRY(ret=XLALInferenceBurstSineGaussianFFast(&hptilde, &hctilde, quality,freq, corrected_hrss, polar_ecc, polar_angle,deltaF,deltaT), errnum);
2239  if (ret == XLAL_FAILURE)
2240  {
2241  XLALPrintError(" ERROR in LALInferenceTemplateXLALSimBurstChooseWaveform(). errnum=%d: %s\n",errnum,XLALErrorString(errnum) );
2242  return;
2243  }
2244  if (hptilde==NULL || hptilde->data==NULL || hptilde->data->data==NULL ) {
2245  XLALPrintError(" ERROR in LALInferenceTemplateXLALInferenceBurstChooseWaveform: encountered unallocated 'hptilde'.\n");
2246  XLAL_ERROR_VOID(XLAL_EFAULT);
2247  }
2248  if (hctilde==NULL || hctilde->data==NULL || hctilde->data->data==NULL ) {
2249  XLALPrintError(" ERROR in LALInferenceTemplateXLALInferenceBurstChooseWaveform: encountered unallocated 'hctilde'.\n");
2250  XLAL_ERROR_VOID(XLAL_EFAULT);
2251  }
2252  size_t lower =(size_t) ( hctilde->f0/hctilde->deltaF);
2253  size_t upper= (size_t) ( hctilde->data->length + lower);
2254  COMPLEX16 *dataPtrP = hptilde->data->data;
2255  COMPLEX16 *dataPtrC = hctilde->data->data;
2256  for (i=0; i<lower; ++i) {
2257  model->freqhPlus->data->data[i] = 0.0;
2258  model->freqhCross->data->data[i] = 0.0;
2259  }
2260 
2261  if (upper>model->freqhPlus->data->length)
2262  upper=model->freqhPlus->data->length;
2263  else{
2264  for (i=upper; i<model->freqhPlus->data->length; ++i) {
2265  model->freqhPlus->data->data[i] = 0.0;
2266  model->freqhCross->data->data[i] = 0.0;
2267  }
2268  }
2269  for (i=lower; i<upper; ++i) {
2270  model->freqhPlus->data->data[i] = dataPtrP[i-lower];
2271  model->freqhCross->data->data[i] = dataPtrC[i-lower];
2272  }
2273 
2274  instant= (model->timehCross->epoch.gpsSeconds + 1e-9*model->timehCross->epoch.gpsNanoSeconds);
2275  LALInferenceSetVariable(model->params, "time", &instant);
2276 
2277 
2278  if ( hplus ) XLALDestroyREAL8TimeSeries(hplus);
2279  if ( hcross ) XLALDestroyREAL8TimeSeries(hcross);
2280  if ( hptilde ) XLALDestroyCOMPLEX16FrequencySeries(hptilde);
2281  if ( hctilde ) XLALDestroyCOMPLEX16FrequencySeries(hctilde);
2282 
2283  return;
2284 }
2285 
2286 
2287 void LALInferenceDumptemplateFreqDomain(LALInferenceVariables *currentParams,
2288  LALInferenceModel *model,
2289  const char *filename)
2290 /* de-bugging function writing (frequency-domain) template to a CSV file */
2291 /* File contains real & imaginary parts of plus & cross components. */
2292 /* Template amplitude is scaled to 1Mpc distance. */
2293 {
2294  FILE *outfile=NULL;
2295  REAL8 f;
2296  UINT4 i;
2297 
2298  REAL8 deltaF = model->deltaF;
2299 
2300  LALInferenceCopyVariables(currentParams, model->params);
2301  model->templt(model);
2302  if (model->domain == LAL_SIM_DOMAIN_TIME)
2303  LALInferenceExecuteFT(model);
2304 
2305  outfile = fopen(filename, "w");
2306  /*fprintf(outfile, "f PSD dataRe dataIm signalPlusRe signalPlusIm signalCrossRe signalCrossIm\n");*/
2307  fprintf(outfile, "\"f\",\"signalPlusRe\",\"signalPlusIm\",\"signalCrossRe\",\"signalCrossIm\"\n");
2308  for (i=0; i<model->freqhPlus->data->length; ++i){
2309  f = ((double) i) * deltaF;
2310  fprintf(outfile, "%f,%e,%e,%e,%e\n",
2311  f,
2312  creal(model->freqhPlus->data->data[i]),
2313  cimag(model->freqhPlus->data->data[i]),
2314  creal(model->freqhCross->data->data[i]),
2315  cimag(model->freqhCross->data->data[i]));
2316  }
2317  fclose(outfile);
2318  fprintf(stdout, " wrote (frequency-domain) template to CSV file \"%s\".\n", filename);
2319 }
2320 
2321 
2322 void LALInferenceDumptemplateTimeDomain(LALInferenceVariables *currentParams,
2323  LALInferenceModel *model,
2324  const char *filename)
2325 /* de-bugging function writing (frequency-domain) template to a CSV file */
2326 /* File contains real & imaginary parts of plus & cross components. */
2327 /* Template amplitude is scaled to 1Mpc distance. */
2328 {
2329  FILE *outfile=NULL;
2330  REAL8 deltaT, t, epoch; // deltaF - set but not used
2331  UINT4 i;
2332 
2333  LALInferenceCopyVariables(currentParams, model->params);
2334  model->templt(model);
2335  if (model->domain == LAL_SIM_DOMAIN_FREQUENCY)
2336  LALInferenceExecuteInvFT(model);
2337 
2338  outfile = fopen(filename, "w");
2339  fprintf(outfile, "\"t\",\"signalPlus\",\"signalCross\"\n");
2340  deltaT = model->deltaT;
2341  epoch = XLALGPSGetREAL8(&model->timehPlus->epoch);
2342  for (i=0; i<model->timehPlus->data->length; ++i){
2343  t = epoch + ((double) i) * deltaT;
2344  fprintf(outfile, "%f,%e,%e\n",
2345  t,
2346  model->timehPlus->data->data[i],
2347  model->timehCross->data->data[i]);
2348  }
2349  fclose(outfile);
2350  fprintf(stdout, " wrote (time-domain) template to CSV file \"%s\".\n", filename);
2351 }
XLALDictEntryGetKey
const char * XLALDictEntryGetKey(const LALDictEntry *entry)
XLALDictInsert
int XLALDictInsert(LALDict *dict, const char *key, const void *data, size_t size, LALTYPECODE type)
XLALDictIterNext
LALDictEntry * XLALDictIterNext(LALDictIter *iter)
XLALDictIterInit
void XLALDictIterInit(LALDictIter *iter, LALDict *dict)
XLALDictEntryGetValue
const LALValue * XLALDictEntryGetValue(const LALDictEntry *entry)
XLALDictInsertValue
int XLALDictInsertValue(LALDict *dict, const char *key, const LALValue *value)
XLALCreateDict
LALDict * XLALCreateDict(void)
XLALSimBurstChooseFDWaveformFromCache
int XLALSimBurstChooseFDWaveformFromCache(COMPLEX16FrequencySeries **hptilde, COMPLEX16FrequencySeries **hctilde, REAL8 deltaF, REAL8 deltaT, REAL8 f0, REAL8 q, REAL8 tau, REAL8 f_min, REAL8 f_max, REAL8 hrss, REAL8 polar_angle, REAL8 polar_ecc, LALSimBurstExtraParam *extraParams, BurstApproximant approximant, LALSimBurstWaveformCache *cache)
Chooses between different approximants when requesting a waveform to be generated Returns the wavefor...
Definition: LALInferenceBurstRoutines.c:1161
XLALSimBurstDestroyExtraParam
void XLALSimBurstDestroyExtraParam(LALSimBurstExtraParam *parameter)
Function that destroys the whole burst extra params linked list.
Definition: LALInferenceBurstRoutines.c:1634
XLALInferenceBurstSineGaussianFFast
int XLALInferenceBurstSineGaussianFFast(COMPLEX16FrequencySeries **hplus, COMPLEX16FrequencySeries **hcross, REAL8 Q, REAL8 centre_frequency, REAL8 hrss, REAL8 eccentricity, REAL8 phase, REAL8 deltaF, REAL8 deltaT)
Definition: LALInferenceBurstRoutines.c:558
XLALSimBurstChooseTDWaveformFromCache
int XLALSimBurstChooseTDWaveformFromCache(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 deltaT, REAL8 f0, REAL8 q, REAL8 tau, REAL8 f_min, REAL8 f_max, REAL8 hrss, REAL8 polar_angle, REAL8 polar_ecc, LALSimBurstExtraParam *extraParams, BurstApproximant approximant, LALSimBurstWaveformCache *cache)
Chooses between different approximants when requesting a waveform to be generated Returns the wavefor...
Definition: LALInferenceBurstRoutines.c:1056
BurstApproximant
BurstApproximant
Enum that specifies the PN approximant to be used in computing the waveform.
Definition: LALInferenceBurstRoutines.h:53
LALInferenceMultibandFrequencies
REAL8Sequence * LALInferenceMultibandFrequencies(int NBands, double f_min, double f_max, double deltaF0, double mc)
Create a list of frequencies to use in multiband template generation, between f_min and f_max mc is m...
Definition: LALInferenceMultibanding.c:41
InterpolateWaveform
static int InterpolateWaveform(REAL8Vector *freqs, COMPLEX16FrequencySeries *src, COMPLEX16FrequencySeries *dest)
Definition: LALInferenceTemplate.c:105
list_FTA_parameters
const char list_FTA_parameters[26][16]
Definition: LALInferenceTemplate.c:83
dquadmon_from_lambda
static REAL8 dquadmon_from_lambda(REAL8 lambdav)
Definition: LALInferenceTemplate.c:95
list_extra_parameters
const char list_extra_parameters[76][16]
Definition: LALInferenceTemplate.c:79
mc2masses
static void mc2masses(double mc, double eta, double *m1, double *m2)
Definition: LALInferenceTemplate.c:212
N_FTA_params
const UINT4 N_FTA_params
Definition: LALInferenceTemplate.c:85
N_extra_params
const UINT4 N_extra_params
Definition: LALInferenceTemplate.c:81
q2masses
static void q2masses(double mc, double q, double *m1, double *m2)
Definition: LALInferenceTemplate.c:224
currentParams
LALInferenceVariables currentParams
Definition: LALInferenceTest.c:368
j
int j
Definition: LALInferenceTest.c:381
k
int k
Definition: LALInferenceTest.c:381
tEnd
#define tEnd
XLALSimInspiralGetStringFromApproximant
const char * XLALSimInspiralGetStringFromApproximant(Approximant approximant)
XLALSimInspiralfLow2fStart
REAL8 XLALSimInspiralfLow2fStart(REAL8 fLow, INT4 ampOrder, INT4 approximant)
gamma
#define gamma
XLALSimInspiralTestingGRCorrections
int XLALSimInspiralTestingGRCorrections(COMPLEX16FrequencySeries *htilde, const UINT4 l, const UINT4 m, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi1z, const REAL8 chi2z, const REAL8 f_low, const REAL8 f_ref, const REAL8 f_window_div_f_Peak, const REAL8 NCyclesStep, LALDict *LALpars)
XLALSimInspiralChooseFDWaveformSequence
int XLALSimInspiralChooseFDWaveformSequence(COMPLEX16FrequencySeries **hptilde, COMPLEX16FrequencySeries **hctilde, REAL8 phiRef, REAL8 m1, REAL8 m2, REAL8 S1x, REAL8 S1y, REAL8 S1z, REAL8 S2x, REAL8 S2y, REAL8 S2z, REAL8 f_ref, REAL8 distance, REAL8 inclination, LALDict *LALpars, Approximant approximant, REAL8Sequence *frequencies)
XLALSimInspiralWaveformParamsInsertPNPhaseOrder
int XLALSimInspiralWaveformParamsInsertPNPhaseOrder(LALDict *params, INT4 value)
XLALSimInspiralWaveformParamsPNAmplitudeOrderIsDefault
int XLALSimInspiralWaveformParamsPNAmplitudeOrderIsDefault(LALDict *params)
XLALSimInspiralWaveformParamsInsertTidalLambda1
int XLALSimInspiralWaveformParamsInsertTidalLambda1(LALDict *params, REAL8 value)
XLALSimInspiralWaveformParamsLookupPNPhaseOrder
INT4 XLALSimInspiralWaveformParamsLookupPNPhaseOrder(LALDict *params)
XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder
INT4 XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(LALDict *params)
XLALSimInspiralWaveformParamsInsertTidalLambda2
int XLALSimInspiralWaveformParamsInsertTidalLambda2(LALDict *params, REAL8 value)
XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder
int XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(LALDict *params, INT4 value)
XLALSimInspiralWaveformParamsInsertdQuadMon2
int XLALSimInspiralWaveformParamsInsertdQuadMon2(LALDict *params, REAL8 value)
XLALSimInspiralWaveformParamsInsertdQuadMon1
int XLALSimInspiralWaveformParamsInsertdQuadMon1(LALDict *params, REAL8 value)
LALSimNeutronStarFamily
struct tagLALSimNeutronStarFamily LALSimNeutronStarFamily
fprintf
#define fprintf
e
double e
duration
double duration
r2
const double r2
data
sigmaKerr data[0]
XLALResizeCOMPLEX16FrequencySeries
COMPLEX16FrequencySeries * XLALResizeCOMPLEX16FrequencySeries(COMPLEX16FrequencySeries *series, int first, size_t length)
XLALCreateCOMPLEX16FrequencySeries
COMPLEX16FrequencySeries * XLALCreateCOMPLEX16FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
XLALDestroyCOMPLEX16FrequencySeries
void XLALDestroyCOMPLEX16FrequencySeries(COMPLEX16FrequencySeries *series)
LAL_MSUN_SI
#define LAL_MSUN_SI
LAL_PI
#define LAL_PI
LAL_TWOPI
#define LAL_TWOPI
LAL_MTSUN_SI
#define LAL_MTSUN_SI
LAL_PC_SI
#define LAL_PC_SI
LAL_SQRT2
#define LAL_SQRT2
LAL_MRSUN_SI
#define LAL_MRSUN_SI
COMPLEX16
double complex COMPLEX16
REAL8
double REAL8
UINT4
uint32_t UINT4
INT4
int32_t INT4
LAL_D_TYPE_CODE
LAL_D_TYPE_CODE
LALInferencedQuadMonSdQuadMonA
void LALInferencedQuadMonSdQuadMonA(REAL8 dQuadMonS, REAL8 dQuadMonA, REAL8 *dQuadMon1, REAL8 *dQuadMon2)
Convert from dQuadMonS, dQuadMonA to dQuadMon1, dQuadMon2.
Definition: LALInference.c:2322
LALInferenceAddVariable
void LALInferenceAddVariable(LALInferenceVariables *vars, const char *name, const void *value, LALInferenceVariableType type, LALInferenceParamVaryType vary)
Add a variable named name to vars with initial value referenced by value.
Definition: LALInference.c:395
LALInferenceLogp1GammasMasses2Lambdas
void LALInferenceLogp1GammasMasses2Lambdas(REAL8 logp1, REAL8 gamma1, REAL8 gamma2, REAL8 gamma3, REAL8 mass1, REAL8 mass2, REAL8 *lambda1, REAL8 *lambda2)
Calculate lambda1,2(m1,2|eos(logp1,gamma1,gamma2,gamma3))
Definition: LALInference.c:2340
LALInferenceGetREAL8Variable
REAL8 LALInferenceGetREAL8Variable(LALInferenceVariables *vars, const char *name)
Definition: LALInference.c:4039
LALInferenceExecuteFT
void LALInferenceExecuteFT(LALInferenceModel *model)
Execute FFT for data in IFOdata.
Definition: LALInference.c:1945
LALInferenceLambdaTsEta2Lambdas
void LALInferenceLambdaTsEta2Lambdas(REAL8 lambdaT, REAL8 dLambdaT, REAL8 eta, REAL8 *lambda1, REAL8 *lambda2)
Convert from lambdaT, dLambdaT, and eta to lambda1 and lambda2.
Definition: LALInference.c:2329
LALInferenceCopyVariables
void LALInferenceCopyVariables(LALInferenceVariables *origin, LALInferenceVariables *target)
Deep copy the variables from one to another LALInferenceVariables structure.
Definition: LALInference.c:558
LALInferenceSDGammasMasses2Lambdas
void LALInferenceSDGammasMasses2Lambdas(REAL8 gamma[], REAL8 mass1, REAL8 mass2, REAL8 *lambda1, REAL8 *lambda2, int size)
Convert from spectral parameters to lambda1, lambda2.
Definition: LALInference.c:2370
LALInferenceBinaryLove
void LALInferenceBinaryLove(LALInferenceVariables *vars, REAL8 *lambda1, REAL8 *lambda2)
Compute Tidal deformabilities following BinaryLove Universal relations.
Definition: LALInference.c:4488
LALInferenceGetVariable
void * LALInferenceGetVariable(const LALInferenceVariables *vars, const char *name)
Return a pointer to the memory the variable vars is stored in specified by name User must cast this p...
Definition: LALInference.c:238
LALInferenceExecuteInvFT
void LALInferenceExecuteInvFT(LALInferenceModel *model)
Execute Inverse FFT for data in IFOdata.
Definition: LALInference.c:2038
LALInferenceCheckVariable
int LALInferenceCheckVariable(LALInferenceVariables *vars, const char *name)
Checks for name being present in vars returns 1(==true) or 0.
Definition: LALInference.c:525
LALInferenceSetVariable
void LALInferenceSetVariable(LALInferenceVariables *vars, const char *name, const void *value)
Set a variable named name in vars with a value.
Definition: LALInference.c:352
LALINFERENCE_PARAM_OUTPUT
@ LALINFERENCE_PARAM_OUTPUT
A parameter that never changes, functions should respect this.
Definition: LALInference.h:131
LALINFERENCE_REAL8_t
@ LALINFERENCE_REAL8_t
Definition: LALInference.h:109
LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence
void LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence(LALInferenceModel *model)
Definition: LALInferenceTemplate.c:234
LALInferenceTemplateXLALSimBurstSineGaussianF
void LALInferenceTemplateXLALSimBurstSineGaussianF(LALInferenceModel *model)
Definition: LALInferenceTemplate.c:2181
LALInferenceTemplateASinOmegaT
void LALInferenceTemplateASinOmegaT(LALInferenceModel *model)
Trivial h(t) = A*sin(Omega*t) template.
Definition: LALInferenceTemplate.c:625
LALInferenceTemplateNullFreqdomain
void LALInferenceTemplateNullFreqdomain(LALInferenceModel *model)
Returns a frequency-domain 'null' template (all zeroes, implying no signal present).
Definition: LALInferenceTemplate.c:165
LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated
void LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated(LALInferenceModel *model)
Definition: LALInferenceTemplate.c:1357
LALInferenceDumptemplateFreqDomain
void LALInferenceDumptemplateFreqDomain(LALInferenceVariables *currentParams, LALInferenceModel *model, const char *filename)
De-bugging function writing a (frequency-domain) signal template to a CSV file.
Definition: LALInferenceTemplate.c:2287
LALInferenceDumptemplateTimeDomain
void LALInferenceDumptemplateTimeDomain(LALInferenceVariables *currentParams, LALInferenceModel *model, const char *filename)
De-bugging function writing a (time-domain) signal template to a CSV file.
Definition: LALInferenceTemplate.c:2322
LALInferenceTemplateSineGaussian
void LALInferenceTemplateSineGaussian(LALInferenceModel *model)
Sine-Gaussian (burst) template.
Definition: LALInferenceTemplate.c:491
LALInferenceTemplateXLALSimBurstChooseWaveform
void LALInferenceTemplateXLALSimBurstChooseWaveform(LALInferenceModel *model)
Definition: LALInferenceTemplate.c:1891
LALInferenceTemplateXLALSimInspiralChooseWaveform
void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
"XLALSimInspiralChooseWaveform{TD,FD}" wrapper.
Definition: LALInferenceTemplate.c:648
LALInferenceTemplateSinc
void LALInferenceTemplateSinc(LALInferenceModel *model)
Sinc function (burst) template.
Definition: LALInferenceTemplate.c:586
LALInferenceTemplateNullTimedomain
void LALInferenceTemplateNullTimedomain(LALInferenceModel *model)
Returns a time-domain 'null' template (all zeroes, implying no signal present).
Definition: LALInferenceTemplate.c:186
LALInferenceTemplateDampedSinusoid
void LALInferenceTemplateDampedSinusoid(LALInferenceModel *model)
Damped Sinusoid template.
Definition: LALInferenceTemplate.c:540
XLALSimInspiralChooseFDModes
SphHarmFrequencySeries * XLALSimInspiralChooseFDModes(REAL8 m1, REAL8 m2, REAL8 S1x, REAL8 S1y, REAL8 S1z, REAL8 S2x, REAL8 S2y, REAL8 S2z, REAL8 deltaF, REAL8 f_min, REAL8 f_max, REAL8 f_ref, REAL8 phiRef, REAL8 distance, REAL8 inclination, LALDict *params, Approximant approximant)
Approximant
Approximant
LAL_SIM_DOMAIN_TIME
LAL_SIM_DOMAIN_TIME
LAL_SIM_DOMAIN_FREQUENCY
LAL_SIM_DOMAIN_FREQUENCY
TaylorF2
TaylorF2
XLALSimInspiralChooseTDWaveformFromCache
int XLALSimInspiralChooseTDWaveformFromCache(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 phiRef, REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 S1x, REAL8 S1y, REAL8 S1z, REAL8 S2x, REAL8 S2y, REAL8 S2z, REAL8 f_min, REAL8 f_ref, REAL8 r, REAL8 i, LALDict *LALpars, Approximant approximant, LALSimInspiralWaveformCache *cache)
XLALSimInspiralChooseFDWaveformFromCache
int XLALSimInspiralChooseFDWaveformFromCache(COMPLEX16FrequencySeries **hptilde, COMPLEX16FrequencySeries **hctilde, REAL8 phiRef, REAL8 deltaF, REAL8 m1, REAL8 m2, REAL8 S1x, REAL8 S1y, REAL8 S1z, REAL8 S2x, REAL8 S2y, REAL8 S2z, REAL8 f_min, REAL8 f_max, REAL8 f_ref, REAL8 r, REAL8 i, LALDict *LALpars, Approximant approximant, LALSimInspiralWaveformCache *cache, REAL8Sequence *frequencies)
XLALSimNeutronStarLoveNumberK2
double XLALSimNeutronStarLoveNumberK2(double m, LALSimNeutronStarFamily *fam)
XLALSimNeutronStarRadius
double XLALSimNeutronStarRadius(double m, LALSimNeutronStarFamily *fam)
XLALSimNeutronStarFamMinimumMass
double XLALSimNeutronStarFamMinimumMass(LALSimNeutronStarFamily *fam)
XLALSimNeutronStarMaximumMass
double XLALSimNeutronStarMaximumMass(LALSimNeutronStarFamily *fam)
XLALSimAddModeFD
int XLALSimAddModeFD(COMPLEX16FrequencySeries *hptilde, COMPLEX16FrequencySeries *hctilde, COMPLEX16FrequencySeries *hlmtilde, REAL8 theta, REAL8 phi, INT4 l, INT4 m, INT4 sym)
XLALSphHarmFrequencySeriesGetMode
COMPLEX16FrequencySeries * XLALSphHarmFrequencySeriesGetMode(SphHarmFrequencySeries *ts, UINT4 l, INT4 m)
XLALDestroySphHarmFrequencySeries
void XLALDestroySphHarmFrequencySeries(SphHarmFrequencySeries *ts)
q
static const INT4 q
a
static const INT4 a
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
XLAL_ERROR_VOID
#define XLAL_ERROR_VOID(...)
XLALErrorString
const char * XLALErrorString(int errnum)
XLALSetErrno
int XLALSetErrno(int errnum)
XLAL_TRY
#define XLAL_TRY(statement, errnum)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_PRINT_ERROR
#define XLAL_PRINT_ERROR(...)
XLAL_EDATA
XLAL_EDATA
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFAULT
XLAL_EUSR0
XLAL_EUSR0
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_FAILURE
XLAL_FAILURE
XLALGPSGetREAL8
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
XLALSimInspiralTransformPrecessingNewInitialConditions
int XLALSimInspiralTransformPrecessingNewInitialConditions(REAL8 *incl, REAL8 *S1x, REAL8 *S1y, REAL8 *S1z, REAL8 *S2x, REAL8 *S2y, REAL8 *S2z, const REAL8 thetaJN, const REAL8 phiJL, const REAL8 theta1, const REAL8 theta2, const REAL8 phi12, const REAL8 chi1, const REAL8 chi2, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 fRef, const REAL8 phiRef)
bayestar_cosmology.filename
filename
Definition: bayestar_cosmology.py:84
cbcBayesBurstPostProc.dest
dest
Definition: cbcBayesBurstPostProc.py:1198
src
src
lalinference.imrtgr.pneqns.Omega
def Omega(v, m1, m2, S1, S2, Ln)
Preccesion frequency spins Eqs.
Definition: pneqns.py:17
lalinference_compute_roq_weights.i
int i
Definition: lalinference_compute_roq_weights.py:269
lalinference_compute_roq_weights.deltaF
int deltaF
Definition: lalinference_compute_roq_weights.py:278
lalinference_evolve_spins_and_append_samples.f_start
f_start
Definition: lalinference_evolve_spins_and_append_samples.py:127
lalinference_evolve_spins_and_append_samples.tilt1
tilt1
Definition: lalinference_evolve_spins_and_append_samples.py:109
lalinference_evolve_spins_and_append_samples.phi12
phi12
Definition: lalinference_evolve_spins_and_append_samples.py:113
lalinference_evolve_spins_and_append_samples.lambda1
lambda1
Definition: lalinference_evolve_spins_and_append_samples.py:140
lalinference_evolve_spins_and_append_samples.lambda2
lambda2
Definition: lalinference_evolve_spins_and_append_samples.py:140
lalinference_evolve_spins_and_append_samples.deltaT
deltaT
Definition: lalinference_evolve_spins_and_append_samples.py:140
lalinference_evolve_spins_and_append_samples.tilt2
tilt2
Definition: lalinference_evolve_spins_and_append_samples.py:109
COMPLEX16FrequencySeries::data
COMPLEX16Sequence * data
COMPLEX16Vector::data
COMPLEX16 * data
LALInferenceModel
Structure to constain a model and its parameters.
Definition: LALInference.h:436
LALInferenceModel::freqhPlus
COMPLEX16FrequencySeries * freqhPlus
Time series model buffers.
Definition: LALInference.h:454
LALInferenceModel::eos_fam
LALSimNeutronStarFamily * eos_fam
Is ROQ enabled.
Definition: LALInference.h:465
LALInferenceModel::fHigh
REAL8 fHigh
Start frequency for waveform generation.
Definition: LALInference.h:449
LALInferenceModel::timehCross
REAL8TimeSeries * timehCross
Definition: LALInference.h:453
LALInferenceModel::burstWaveformCache
LALSimBurstWaveformCache * burstWaveformCache
Waveform cache.
Definition: LALInference.h:459
LALInferenceModel::waveformCache
LALSimInspiralWaveformCache * waveformCache
Definition: LALInference.h:458
LALInferenceModel::timehPlus
REAL8TimeSeries * timehPlus
Definition: LALInference.h:453
LALInferenceModel::padding
REAL8 padding
A window.
Definition: LALInference.h:462
LALInferenceModel::params
LALInferenceVariables * params
Definition: LALInference.h:437
LALInferenceModel::deltaT
REAL8 deltaT
End frequency for waveform generation.
Definition: LALInference.h:450
LALInferenceModel::freqhCross
COMPLEX16FrequencySeries * freqhCross
Definition: LALInference.h:454
LALInferenceModel::LALpars
LALDict * LALpars
Projected freq series model buffers.
Definition: LALInference.h:457
LALInferenceModel::window
REAL8Window * window
Pre-calculated FFT plans for forward and reverse FFTs.
Definition: LALInference.h:461
LALInferenceModel::fLow
REAL8 fLow
Array of single-IFO SNRs at params
Definition: LALInference.h:448
LALInferenceModel::roq
struct tagLALInferenceROQModel * roq
The padding of the above window.
Definition: LALInference.h:463
LALInferenceModel::templt
LALInferenceTemplateFunction templt
Domain of model.
Definition: LALInference.h:440
LALInferenceModel::domain
LALSimulationDomain domain
IFO-dependent parameters and buffers.
Definition: LALInference.h:439
LALInferenceVariables
The LALInferenceVariables structure to contain a set of parameters Implemented as a linked list of LA...
Definition: LALInference.h:170
LALSimBurstExtraParam
Linked list of any number of parameters for testing GR.
Definition: LALInferenceBurstRoutines.h:169
LIGOTimeGPS::gpsNanoSeconds
INT4 gpsNanoSeconds
REAL8TimeSeries::data
REAL8Sequence * data
REAL8TimeSeries::epoch
LIGOTimeGPS epoch
REAL8Vector::data
REAL8 * data
REAL8Window::data
REAL8Sequence * data
REAL8Window::sumofsquares
REAL8 sumofsquares
SphHarmFrequencySeries::next
struct tagSphHarmFrequencySeries * next
SphHarmFrequencySeries::mode
COMPLEX16FrequencySeries * mode
epoch
enum @1 epoch
hrss
double hrss
f_max
double f_max
df
double df
outfile
char * outfile
pad
double pad