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