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