LALSimInspiral.c 317 KB
Newer Older
1
/*
2
 * Copyright (C) 2008 J. Creighton, S. Fairhurst, B. Krishnan, L. Santamaria, D. Keppel, Evan Ochsner, C. Pankow, 2014 A. Klein
3
4
5
6
7
8
9
10
11
12
13
14
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
Adam Mercer's avatar
Adam Mercer committed
15
16
 *  along with with program; see the file COPYING. If not, write to the
 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
17
18
19
 *  MA  02111-1307  USA
 */

20
#include <complex.h>
21
22
23
24
25
26
27
#include <math.h>

#include <gsl/gsl_const.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_odeiv.h>

28
#include <lal/SphericalHarmonics.h>
29
#include <lal/LALSimInspiral.h>
30
#include <lal/LALSimIMR.h>
31
#include <lal/LALSimSphHarmMode.h>
32
33
#include <lal/LALConstants.h>
#include <lal/LALStdlib.h>
34
#include <lal/LALString.h>
35
#include <lal/Sequence.h>
36
#include <lal/TimeSeries.h>
37
#include <lal/FrequencySeries.h>
38
#include <lal/TimeFreqFFT.h>
39
#include <lal/BandPassTimeSeries.h>
40
#include <lal/Units.h>
41
#include <lal/LALSimBlackHoleRingdown.h>
42

43
#include "LALSimInspiralPNCoefficients.c"
44
#include "check_series_macros.h"
45
#include "check_waveform_macros.h"
46

47
48
49
50
51
52
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

53
54
55
56
/**
 * (Twice) the highest known PN order of amplitude correction for
 * non-precessing binaries.
 */
57
#define MAX_NONPRECESSING_AMP_PN_ORDER 6
58
59
60
61
62

/**
 * (Twice) the highest known PN order of amplitude correction for
 * precessing binaries.
 */
63
#define MAX_PRECESSING_AMP_PN_ORDER 3
64

65
/* Macro functions to rotate the components of a vector about an axis */
66
67
68
69
#define ROTATEZ(angle, vx, vy, vz)\
	tmp1 = vx*cos(angle) - vy*sin(angle);\
	tmp2 = vx*sin(angle) + vy*cos(angle);\
	vx = tmp1;\
70
	vy = tmp2
71
72

#define ROTATEY(angle, vx, vy, vz)\
73
74
	tmp1 = vx*cos(angle) + vz*sin(angle);\
	tmp2 = - vx*sin(angle) + vz*cos(angle);\
75
	vx = tmp1;\
76
	vz = tmp2
77

78

79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#define INITIALIZE_NAME(a) [a] = #a
/* TODO: UPDATE ME WHENEVER A NEW APPROXIMANT IS ADDED */
static const char *lalSimulationApproximantNames[] = {
    INITIALIZE_NAME(TaylorT1),
    INITIALIZE_NAME(TaylorT2),
    INITIALIZE_NAME(TaylorT3),
    INITIALIZE_NAME(TaylorF1),
    INITIALIZE_NAME(TaylorF2),
    INITIALIZE_NAME(TaylorR2F4),
    INITIALIZE_NAME(TaylorF2RedSpin),
    INITIALIZE_NAME(TaylorF2RedSpinTidal),
    INITIALIZE_NAME(PadeT1),
    INITIALIZE_NAME(PadeF1),
    INITIALIZE_NAME(EOB),
    INITIALIZE_NAME(BCV),
    INITIALIZE_NAME(BCVSpin),
    INITIALIZE_NAME(SpinTaylorT1),
    INITIALIZE_NAME(SpinTaylorT2),
    INITIALIZE_NAME(SpinTaylorT3),
    INITIALIZE_NAME(SpinTaylorT4),
    INITIALIZE_NAME(SpinTaylorT5),
    INITIALIZE_NAME(SpinTaylorF2),
    INITIALIZE_NAME(SpinTaylorFrameless),
    INITIALIZE_NAME(SpinTaylor),
    INITIALIZE_NAME(PhenSpinTaylor),
    INITIALIZE_NAME(PhenSpinTaylorRD),
    INITIALIZE_NAME(SpinQuadTaylor),
    INITIALIZE_NAME(FindChirpSP),
    INITIALIZE_NAME(FindChirpPTF),
    INITIALIZE_NAME(GeneratePPN),
    INITIALIZE_NAME(BCVC),
    INITIALIZE_NAME(FrameFile),
    INITIALIZE_NAME(AmpCorPPN),
    INITIALIZE_NAME(NumRel),
    INITIALIZE_NAME(NumRelNinja2),
    INITIALIZE_NAME(EccentricFD),
    INITIALIZE_NAME(Eccentricity),
    INITIALIZE_NAME(EOBNR),
    INITIALIZE_NAME(EOBNRv2),
    INITIALIZE_NAME(EOBNRv2HM),
119
120
    INITIALIZE_NAME(EOBNRv2_ROM),
    INITIALIZE_NAME(EOBNRv2HM_ROM),
121
122
    INITIALIZE_NAME(SEOBNRv1),
    INITIALIZE_NAME(SEOBNRv2),
123
    INITIALIZE_NAME(SEOBNRv2_opt),
124
    INITIALIZE_NAME(SEOBNRv3),
125
    INITIALIZE_NAME(SEOBNRv3_pert),
126
127
    INITIALIZE_NAME(SEOBNRv3_opt),
    INITIALIZE_NAME(SEOBNRv3_opt_rk4),
128
129
    INITIALIZE_NAME(SEOBNRv4),
    INITIALIZE_NAME(SEOBNRv4_opt),
130
131
132
133
    INITIALIZE_NAME(SEOBNRv1_ROM_EffectiveSpin),
    INITIALIZE_NAME(SEOBNRv1_ROM_DoubleSpin),
    INITIALIZE_NAME(SEOBNRv2_ROM_EffectiveSpin),
    INITIALIZE_NAME(SEOBNRv2_ROM_DoubleSpin),
134
    INITIALIZE_NAME(SEOBNRv2_ROM_DoubleSpin_HI),
135
    INITIALIZE_NAME(Lackey_Tidal_2013_SEOBNRv2_ROM),
136
    INITIALIZE_NAME(SEOBNRv4_ROM),
137
138
139
140
141
142
143
144
    INITIALIZE_NAME(HGimri),
    INITIALIZE_NAME(IMRPhenomA),
    INITIALIZE_NAME(IMRPhenomB),
    INITIALIZE_NAME(IMRPhenomFA),
    INITIALIZE_NAME(IMRPhenomFB),
    INITIALIZE_NAME(IMRPhenomC),
    INITIALIZE_NAME(IMRPhenomD),
    INITIALIZE_NAME(IMRPhenomP),
145
    INITIALIZE_NAME(IMRPhenomPv2),
146
147
148
149
150
151
152
153
    INITIALIZE_NAME(IMRPhenomFC),
    INITIALIZE_NAME(TaylorEt),
    INITIALIZE_NAME(TaylorT4),
    INITIALIZE_NAME(EccentricTD),
    INITIALIZE_NAME(TaylorN),
    INITIALIZE_NAME(SpinTaylorT4Fourier),
    INITIALIZE_NAME(SpinTaylorT2Fourier),
    INITIALIZE_NAME(SpinDominatedWf),
154
    INITIALIZE_NAME(NR_hdf5),
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
};
#undef INITIALIZE_NAME

/* TODO: UPDATE ME WHENEVER A NEW PN ORDER IS ADDED */
static const char *lalSimulationPNOrderNames[] = {
    [LAL_PNORDER_NEWTONIAN]         = "newtonian",
    [LAL_PNORDER_HALF]              = "oneHalfPN",
    [LAL_PNORDER_ONE]               = "onePN",
    [LAL_PNORDER_ONE_POINT_FIVE]    = "onePointFivePN",
    [LAL_PNORDER_TWO]               = "twoPN",
    [LAL_PNORDER_TWO_POINT_FIVE]    = "twoPointFivePN",
    [LAL_PNORDER_THREE]             = "threePN",
    [LAL_PNORDER_THREE_POINT_FIVE]  = "threePointFivePN",
    [LAL_PNORDER_PSEUDO_FOUR]       = "pseudoFourPN",
};

/* TODO: UPDATE ME WHENEVER A NEW TAPER IS ADDED */
static const char *lalSimulationTaperNames[] = {
    [LAL_SIM_INSPIRAL_TAPER_NONE]       = "TAPER_NONE",
    [LAL_SIM_INSPIRAL_TAPER_START]      = "TAPER_START",
    [LAL_SIM_INSPIRAL_TAPER_END]        = "TAPER_END",
    [LAL_SIM_INSPIRAL_TAPER_STARTEND]   = "TAPER_STARTEND",
};

/* TODO: UPDATE ME WHENEVER A NEW FRAME AXIS IS ADDED */
static const char *lalSimulationFrameAxisNames[] = {
    [LAL_SIM_INSPIRAL_FRAME_AXIS_TOTAL_J]   = "TotalJ",
    [LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L] = "OrbitalL",
    [LAL_SIM_INSPIRAL_FRAME_AXIS_VIEW]      = "View",
};

/* TODO: UPDATE ME WHENEVER A NEW MODES CHOICE IS ADDED */
static const char *lalSimulationModesChoiceNames[] = {
    [LAL_SIM_INSPIRAL_MODES_CHOICE_2AND3AND4AND5L] = "L2345",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_2AND3AND4L] = "L234",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_2AND3AND5L] = "L235",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_2AND4AND5L] = "L245",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_3AND4AND5L] = "L345",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_2AND3L] = "L23",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_2AND4L] = "L24",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_3AND4L] = "L34",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_2AND5L] = "L25",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_3AND5L] = "L35",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_4AND5L] = "L45",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_RESTRICTED] = "L2",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_3L] = "L3",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_4L] = "L4",
    [LAL_SIM_INSPIRAL_MODES_CHOICE_5L] = "L5",
    /* NOTE: cannot do the "ALL" case since its value is -1 */
    // [LAL_SIM_INSPIRAL_MODES_CHOICE_ALL] = "ALL",
};

/* locates and deletes a substring in a list of substrings from a string, ignoring case;
 * if multiple substrings in the string match, delete the longest one; here, deletion
 * means replacing the substring with BEL characters */
static int delete_substring_in_list_from_string(char *string, const char *list[], size_t size)
{
    int longest_position = -1;
    int longest_offset = -1;
    int longest_length = -1;
    size_t i;

    if (string == NULL || strlen(string) == 0) // no string to search
        return -1;

    for (i = 0; i < size; ++i) {
        char *match;
        if (list[i] == NULL) // no such element in list
            continue;
        if ((match = XLALStringCaseSubstring(string, list[i]))) {
            int length = strlen(list[i]);
            if (length > longest_length) {
                longest_position = i;
                longest_offset = match - string;
                longest_length = length;
            }
        }
    }

    if (longest_position < 0) // failed to find a word
        return -1;

    /* delete word from string by replacing with BEL */
    for (i = 0; i < (size_t)longest_length; ++i)
        string[longest_offset + i] = '\b';

    return longest_position;
}

244

245
/*
246
247
248
249
250
 * certain approximants adopt the convention that f_ref=0 refers to the start
 * of the waveform while other approximants adopt the convention that f_ref=0
 * refers to the end of the waveform. in the former case, this routine will
 * return the explicit value of f_ref, which will be f_min.
 */
251
static double fixReferenceFrequency(const double f_ref, const double f_min, const Approximant approximant)
252
253
254
255
256
257
258
259
260
261
262
{
    if (f_ref == 0)
        switch (approximant) {
        case SpinTaylorT1:
        case SpinTaylorT2:
        case SpinTaylorT3:
        case SpinTaylorT4:
        case SpinTaylorT2Fourier:
        case SpinTaylorT4Fourier:
        case SpinTaylorF2:
        case IMRPhenomP:
263
        case IMRPhenomPv2:
264
265
266
267
268
269
            return f_min;
        default:
            break;
        }
    return f_ref;
}
270

271
272
273
/*
 * some helper routines for XLALSimInspiralTD
 */
274
275
static int XLALSimInspiralTDFromTD(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 m1, REAL8 m2, REAL8 S1x, REAL8 S1y, REAL8 S1z, REAL8 S2x, REAL8 S2y, REAL8 S2z, REAL8 distance, REAL8 inclination, REAL8 phiRef, REAL8 longAscNodes, REAL8 eccentricity, REAL8 meanPerAno, REAL8 deltaT, REAL8 f_min, REAL8 f_ref, LALDict *LALparams, Approximant approximant);
static int XLALSimInspiralTDFromFD(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 m1, REAL8 m2, REAL8 S1x, REAL8 S1y, REAL8 S1z, REAL8 S2x, REAL8 S2y, REAL8 S2z, REAL8 distance, REAL8 inclination, REAL8 phiRef, REAL8 longAscNodes, REAL8 eccentricity, REAL8 meanPerAno, REAL8 deltaT, REAL8 f_min, REAL8 f_ref, LALDict *LALparams, Approximant approximant);
276

277
/**
278
279
280
281
 * @addtogroup LALSimInspiral_c
 * @brief General routines for generating binary inspiral waveforms.
 *
 * @{
282
283
284
 */

/**
285
286
 * @name General Waveform Switching Generation Routines
 * @{
287
288
289
 */

/**
290
291
292
293
294
 * Chooses between different approximants when requesting a waveform to be generated
 * For spinning waveforms, all known spin effects up to given PN order are included
 * Returns the waveform in the time domain.
 *
 * The parameters passed must be in SI units.
295
 */
296
297
298
int XLALSimInspiralChooseTDWaveform(
    REAL8TimeSeries **hplus,                    /**< +-polarization waveform */
    REAL8TimeSeries **hcross,                   /**< x-polarization waveform */
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
    const REAL8 m1,                             /**< mass of companion 1 (kg) */
    const REAL8 m2,                             /**< mass of companion 2 (kg) */
    const REAL8 S1x,                            /**< x-component of the dimensionless spin of object 1 */
    const REAL8 S1y,                            /**< y-component of the dimensionless spin of object 1 */
    const REAL8 S1z,                            /**< z-component of the dimensionless spin of object 1 */
    const REAL8 S2x,                            /**< x-component of the dimensionless spin of object 2 */
    const REAL8 S2y,                            /**< y-component of the dimensionless spin of object 2 */
    const REAL8 S2z,                            /**< z-component of the dimensionless spin of object 2 */
    const REAL8 distance,                       /**< distance of source (m) */
    const REAL8 inclination,                    /**< inclination of source (rad) */
    const REAL8 phiRef,                         /**< reference orbital phase (rad) */
    const REAL8 longAscNodes,                   /**< longitude of ascending nodes, degenerate with the polarization angle, Omega in documentation */
    const REAL8 eccentricity,                   /**< eccentrocity at reference epoch */
    const REAL8 UNUSED meanPerAno,              /**< mean anomaly of periastron */
    const REAL8 deltaT,                         /**< sampling interval (s) */
    const REAL8 f_min,                          /**< starting GW frequency (Hz) */
315
    REAL8 f_ref,                                /**< reference GW frequency (Hz) */
316
317
    LALDict *LALparams,                         /**< LAL dictionary containing accessory parameters */
    const Approximant approximant               /**< post-Newtonian approximant to use for waveform production */
318
    )
319
{
320
    REAL8 LNhatx, LNhaty, LNhatz, E1x, E1y, E1z;
321
    //REAL8 tmp1, tmp2;
322
    int ret;
323
    /* N.B. the quadrupole of a spinning compact body labeled by A is
324
325
326
327
     * Q_A = - quadparam_A chi_A^2 m_A^3 (see gr-qc/9709032)
     * where quadparam = 1 for BH ~= 4-8 for NS.
     * This affects the quadrupole-monopole interaction.
     */
328
329
330
331
332
333
334
    REAL8 v0 = 1.;
    REAL8 quadparam1 = 1.+XLALSimInspiralWaveformParamsLookupdQuadMon1(LALparams);
    REAL8 quadparam2 = 1.+XLALSimInspiralWaveformParamsLookupdQuadMon2(LALparams);
    REAL8 lambda1 = XLALSimInspiralWaveformParamsLookupTidalLambda1(LALparams);
    REAL8 lambda2 = XLALSimInspiralWaveformParamsLookupTidalLambda2(LALparams);
    int amplitudeO = XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(LALparams);
    int phaseO =XLALSimInspiralWaveformParamsLookupPNPhaseOrder(LALparams);
335

336
337
338
339
340
    /* General sanity checks that will abort
     *
     * If non-GR approximants are added, include them in
     * XLALSimInspiralApproximantAcceptTestGRParams()
     */
341
    if( !XLALSimInspiralWaveformParamsNonGRAreDefault(LALparams) && XLALSimInspiralApproximantAcceptTestGRParams(approximant) != LAL_SIM_INSPIRAL_TESTGR_PARAMS ) {
342
343
        XLALPrintError("XLAL Error - %s: Passed in non-NULL pointer to LALSimInspiralTestGRParam for an approximant that does not use LALSimInspiralTestGRParam\n", __func__);
        XLAL_ERROR(XLAL_EINVAL);
344
    }
345
    /* Support variables for precessing wfs*/
346
    REAL8 incl;
347

348
349

    /* SEOBNR flag for spin aligned model version. 1 for SEOBNRv1, 2 for SEOBNRv2 */
350
    UINT4 SpinAlignedEOBversion;
351
352
353
    REAL8 spin1x,spin1y,spin1z;
    REAL8 spin2x,spin2y,spin2z;
    REAL8 polariz=longAscNodes;
354
355
356
357
358

    /* SEOBNR flag for precessing model version. 3 for SEOBNRv3, 300 for SEOBNRv3_opt */
    UINT4 PrecEOBversion;
    REAL8 spin1[3], spin2[3];

359
    //LIGOTimeGPS epoch = LIGOTIMEGPSZERO;
360

361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
    /* General sanity check the input parameters - only give warnings! */
    if( deltaT > 1. )
        XLALPrintWarning("XLAL Warning - %s: Large value of deltaT = %e requested.\nPerhaps sample rate and time step size were swapped?\n", __func__, deltaT);
    if( deltaT < 1./16385. )
        XLALPrintWarning("XLAL Warning - %s: Small value of deltaT = %e requested.\nCheck for errors, this could create very large time series.\n", __func__, deltaT);
    if( m1 < 0.09 * LAL_MSUN_SI )
        XLALPrintWarning("XLAL Warning - %s: Small value of m1 = %e (kg) = %e (Msun) requested.\nPerhaps you have a unit conversion error?\n", __func__, m1, m1/LAL_MSUN_SI);
    if( m2 < 0.09 * LAL_MSUN_SI )
        XLALPrintWarning("XLAL Warning - %s: Small value of m2 = %e (kg) = %e (Msun) requested.\nPerhaps you have a unit conversion error?\n", __func__, m2, m2/LAL_MSUN_SI);
    if( m1 + m2 > 1000. * LAL_MSUN_SI )
        XLALPrintWarning("XLAL Warning - %s: Large value of total mass m1+m2 = %e (kg) = %e (Msun) requested.\nSignal not likely to be in band of ground-based detectors.\n", __func__, m1+m2, (m1+m2)/LAL_MSUN_SI);
    if( S1x*S1x + S1y*S1y + S1z*S1z > 1.000001 )
        XLALPrintWarning("XLAL Warning - %s: S1 = (%e,%e,%e) with norm > 1 requested.\nAre you sure you want to violate the Kerr bound?\n", __func__, S1x, S1y, S1z);
    if( S2x*S2x + S2y*S2y + S2z*S2z > 1.000001 )
        XLALPrintWarning("XLAL Warning - %s: S2 = (%e,%e,%e) with norm > 1 requested.\nAre you sure you want to violate the Kerr bound?\n", __func__, S2x, S2y, S2z);
    if( f_min < 1. )
        XLALPrintWarning("XLAL Warning - %s: Small value of fmin = %e requested.\nCheck for errors, this could create a very long waveform.\n", __func__, f_min);
    if( f_min > 40.000001 )
        XLALPrintWarning("XLAL Warning - %s: Large value of fmin = %e requested.\nCheck for errors, the signal will start in band.\n", __func__, f_min);
380

381
382
383
384
    /* adjust the reference frequency for certain precessing approximants:
     * if that approximate interprets f_ref==0 to be f_min, set f_ref=f_min;
     * otherwise do nothing */
    f_ref = fixReferenceFrequency(f_ref, f_min, approximant);
385

386
387
388
389
    switch (approximant)
    {
        /* non-spinning inspiral-only models */
        case TaylorEt:
390
            /* Waveform-specific sanity checks */
391
392
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
393
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
394
                ABORT_NONZERO_SPINS(LALparams);
395
            if( !checkTidesZero(lambda1, lambda2) )
396
                ABORT_NONZERO_TIDES(LALparams);
397
398
399
400
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorEtPNGenerator(hplus, hcross, phiRef, v0,
401
                    deltaT, m1, m2, f_min, distance, inclination, amplitudeO, phaseO);
402
403
            break;

404
405
        case TaylorT1:
            /* Waveform-specific sanity checks */
406
407
408
409
410
411
	    if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);
            if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALparams) )
                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
            if( !XLALSimInspiralWaveformParamsPNSpinOrderIsDefault(LALparams) )
                ABORT_NONDEFAULT_SPIN_ORDER(LALparams);
412
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
413
                ABORT_NONZERO_SPINS(LALparams);
414
415
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorT1PNGenerator(hplus, hcross, phiRef, v0,
416
417
                    deltaT, m1, m2, f_min, f_ref, distance, inclination, lambda1, lambda2,
                    XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALparams), amplitudeO, phaseO);
418
419
            break;

420
421
        case TaylorT2:
            /* Waveform-specific sanity checks */
422
423
424
425
426
427
	    if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);
            if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALparams) )
                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
            if( !XLALSimInspiralWaveformParamsPNSpinOrderIsDefault(LALparams) )
                ABORT_NONDEFAULT_SPIN_ORDER(LALparams);
428
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
429
                ABORT_NONZERO_SPINS(LALparams);
430
431
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorT2PNGenerator(hplus, hcross, phiRef, v0,
432
433
                    deltaT, m1, m2, f_min, f_ref, distance, inclination, lambda1, lambda2,
                    XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALparams), amplitudeO, phaseO);
434
435
            break;

436
437
        case TaylorT3:
            /* Waveform-specific sanity checks */
438
439
440
441
442
443
	    if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);
            if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALparams) )
                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
            if( !XLALSimInspiralWaveformParamsPNSpinOrderIsDefault(LALparams) )
                ABORT_NONDEFAULT_SPIN_ORDER(LALparams);
444
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
445
                ABORT_NONZERO_SPINS(LALparams);
446
447
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorT3PNGenerator(hplus, hcross, phiRef, v0,
448
449
                    deltaT, m1, m2, f_min, f_ref, distance, inclination, lambda1, lambda2,
                    XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALparams), amplitudeO, phaseO);
450
            break;
451

452
453
        case TaylorT4:
            /* Waveform-specific sanity checks */
454
455
456
457
458
459
	    if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);
            if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALparams) )
                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
            if( !XLALSimInspiralWaveformParamsPNSpinOrderIsDefault(LALparams) )
                ABORT_NONDEFAULT_SPIN_ORDER(LALparams);
460
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
461
                ABORT_NONZERO_SPINS(LALparams);
462
463
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorT4PNGenerator(hplus, hcross, phiRef, v0,
464
465
                    deltaT, m1, m2, f_min, f_ref, distance, inclination, lambda1, lambda2,
                    XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALparams), amplitudeO, phaseO);
466
467
            break;

468
469
	case EccentricTD:
	    /* Waveform-specific sanity checks */
470
471
472
473
474
475
476
477
	    if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);
            if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALparams) )
                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
            if( !XLALSimInspiralWaveformParamsPNSpinOrderIsDefault(LALparams) )
                ABORT_NONDEFAULT_SPIN_ORDER(LALparams);
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
                ABORT_NONZERO_SPINS(LALparams);
478
	    if( !checkTidesZero(lambda1, lambda2) )
479
		ABORT_NONZERO_TIDES(LALparams);
480
481
	    /* Call the waveform driver routine */
	    ret = XLALSimInspiralEccentricTDPNGenerator(hplus, hcross, phiRef,
482
		    deltaT, m1, m2, f_min, f_ref, distance, inclination, eccentricity,
483
484
485
486
487
		    amplitudeO, phaseO);
	    if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
	    break;

        /* non-spinning inspiral-merger-ringdown models */
488
        case IMRPhenomA:
489
            /* Waveform-specific sanity checks */
490
491
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
492
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
493
                ABORT_NONZERO_SPINS(LALparams);
494
            if( !checkTidesZero(lambda1, lambda2) )
495
                ABORT_NONZERO_TIDES(LALparams);
496
497
498
499
500
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
            // NB: f_max = 0 will generate up to the ringdown cut-off frequency
            ret = XLALSimIMRPhenomAGenerateTD(hplus, hcross, phiRef, deltaT,
501
                    m1, m2, f_min, 0., distance, inclination);
502
503
            break;

504
505
        case EOBNRv2HM:
            /* Waveform-specific sanity checks */
506
507
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
508
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
509
	      ABORT_NONZERO_SPINS(LALparams);
510
            if( !checkTidesZero(lambda1, lambda2) )
511
                ABORT_NONZERO_TIDES(LALparams);
512
513
514
515
516
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
            // FIXME: need to create a function to take in different modes or produce an error if all modes not given
            ret = XLALSimIMREOBNRv2AllModes(hplus, hcross, phiRef, deltaT,
517
                    m1, m2, f_min, distance, inclination);
518
519
            break;

520
521
        case EOBNRv2:
            /* Waveform-specific sanity checks */
522
523
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
524
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
525
                ABORT_NONZERO_SPINS(LALparams);
526
            if( !checkTidesZero(lambda1, lambda2) )
527
                ABORT_NONZERO_TIDES(LALparams);
528
529
530
531
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
            ret = XLALSimIMREOBNRv2DominantMode(hplus, hcross, phiRef, deltaT,
532
                    m1, m2, f_min, distance, inclination);
533
534
            break;

535
        /* spinning inspiral-only models */
536
        case SpinTaylorT2:
537
            /* Waveform-specific sanity checks */
538
539
	    if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) && (XLALSimInspiralWaveformParamsLookupPNSpinOrder(LALparams)>5) )
	      ABORT_NONZERO_TRANSVERSE_SPINS_HIGH_SPINO(LALparams);
540
541
	    XLALSimInspiralInitialConditionsPrecessingApproxs(&incl,&spin1x,&spin1y,&spin1z,&spin2x,&spin2y,&spin2z,inclination,S1x,S1y,S1z,S2x,S2y,S2z,m1,m2,f_ref,phiRef,XLALSimInspiralWaveformParamsLookupFrameAxis(LALparams));
            LNhatx = sin(incl);
542
            LNhaty = 0.;
543
544
545
546
547
            LNhatz = cos(incl);
            E1x = 0.;
            E1y = 1.;
            E1z = 0.;
	    polariz+=LAL_PI/2.;
548
            /* Maximum PN amplitude order for precessing waveforms is
549
             * MAX_PRECESSING_AMP_PN_ORDER */
550
            amplitudeO = amplitudeO <= MAX_PRECESSING_AMP_PN_ORDER ?
551
552
553
                    amplitudeO : MAX_PRECESSING_AMP_PN_ORDER;
            /* Call the waveform driver routine */
            ret = XLALSimInspiralSpinTaylorT2(hplus, hcross, phiRef, v0, deltaT,
554
555
556
557
558
					      m1, m2, f_min, f_ref, distance, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
					      LNhatx, LNhaty, LNhatz, E1x, E1y, E1z, lambda1, lambda2,
					      quadparam1, quadparam2,
					      LALparams,
					      phaseO, amplitudeO);
559
560
            break;

561
562
        // need to make a consistent choice for SpinTaylorT4 and PSpinInspiralRD waveform inputs
        // proposal: TotalJ frame of PSpinInspiralRD
563
564
565
        // inclination denotes the angle between the view direction
        // and J (J is constant during the evolution, J//z, both N and initial
        // L are in the x-z plane) and the spin coordinates are given wrt
566
567
568
        // initial ** L **.
        case SpinTaylorT4:
            /* Waveform-specific sanity checks */
569
570
	    if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) && (XLALSimInspiralWaveformParamsLookupPNSpinOrder(LALparams)>5) )
	      ABORT_NONZERO_TRANSVERSE_SPINS_HIGH_SPINO(LALparams);
571
572
            XLALSimInspiralInitialConditionsPrecessingApproxs(&incl,&spin1x,&spin1y,&spin1z,&spin2x,&spin2y,&spin2z,inclination,S1x,S1y,S1z,S2x,S2y,S2z,m1,m2,f_ref,phiRef,XLALSimInspiralWaveformParamsLookupFrameAxis(LALparams));
            LNhatx = sin(incl);
573
            LNhaty = 0.;
574
575
576
577
578
            LNhatz = cos(incl);
            E1x = 0.;
            E1y = 1.;
            E1z = 0.;
	    polariz+=LAL_PI/2.;
579
            /* Maximum PN amplitude order for precessing waveforms is
580
             * MAX_PRECESSING_AMP_PN_ORDER */
581
            amplitudeO = amplitudeO <= MAX_PRECESSING_AMP_PN_ORDER ?
582
583
584
                    amplitudeO : MAX_PRECESSING_AMP_PN_ORDER;
            /* Call the waveform driver routine */
            ret = XLALSimInspiralSpinTaylorT4(hplus, hcross, phiRef, v0, deltaT,
585
                    m1, m2, f_min, f_ref, distance, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
586
                    LNhatx, LNhaty, LNhatz, E1x, E1y, E1z, lambda1, lambda2,
587
	            quadparam1, quadparam2, LALparams,
588
589
                    phaseO, amplitudeO);
            break;
590

591
592
	 case SpinTaylorT1:
            /* Waveform-specific sanity checks */
593
594
595
            /* Waveform-specific sanity checks */
	    if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) && (XLALSimInspiralWaveformParamsLookupPNSpinOrder(LALparams)>5) )
	      ABORT_NONZERO_TRANSVERSE_SPINS_HIGH_SPINO(LALparams);
596
597
	    XLALSimInspiralInitialConditionsPrecessingApproxs(&incl,&spin1x,&spin1y,&spin1z,&spin2x,&spin2y,&spin2z,inclination,S1x,S1y,S1z,S2x,S2y,S2z,m1,m2,f_ref,phiRef,XLALSimInspiralWaveformParamsLookupFrameAxis(LALparams));
            LNhatx = sin(incl);
598
            LNhaty = 0.;
599
600
601
602
603
            LNhatz = cos(incl);
            E1x = 0.;
            E1y = 1.;
            E1z = 0.;
	    polariz+=LAL_PI/2.;
604
605
606
607
608
609
            /* Maximum PN amplitude order for precessing waveforms is
             * MAX_PRECESSING_AMP_PN_ORDER */
            amplitudeO = amplitudeO <= MAX_PRECESSING_AMP_PN_ORDER ?
                    amplitudeO : MAX_PRECESSING_AMP_PN_ORDER;
            /* Call the waveform driver routine */
            ret = XLALSimInspiralSpinTaylorT1(hplus, hcross, phiRef, v0, deltaT,
610
611
					      m1, m2, f_min, f_ref, distance, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z, LNhatx, LNhaty, LNhatz, E1x, E1y, E1z, lambda1, lambda2,
                    quadparam1, quadparam2, LALparams,
612
613
                    phaseO, amplitudeO);
            break;
614

615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
        case SpinDominatedWf:
                // waveform specific sanity checks
                if (S2x != 0. || S2y != 0. || S2z != 0.){
                XLALPrintError("XLAL Error : The spindominatedwf approximant is only for 1 spin case.\n");
                XLAL_ERROR(XLAL_EDOM);
                }
                /*Maximal PN amplitude order is 1.5, maximal phase order is 2 PN*/
                if (amplitudeO > 3) {
                XLALPrintError("XLAL Error : Foe the spindominatedwf approximant maximal amplitude correction is 1.5 PN\n");
                XLAL_ERROR(XLAL_EDOM);
                }
                if (phaseO > 4){
                XLALPrintError("XLAL Error : For the spindominatedwf approximant maximal phase correction is 2 PN\n");
                XLAL_ERROR(XLAL_EDOM);
                }
630
631
	            incl=inclination;
                LNhatx = 0.;
632
                LNhaty = 0.;
633
                LNhatz = 1.;
634
                /* Call the waveform driver routine */
635
                ret = XLALSimInspiralSpinDominatedWaveformInterfaceTD(hplus, hcross, deltaT, m1, m2, f_min, f_ref, distance, S1x, S1y, S1z, LNhatx, LNhaty, LNhatz, incl, phaseO, amplitudeO, phiRef);
636
                break;
637

638
639
640
        /* spin aligned inspiral-merger-ringdown models */
        case IMRPhenomB:
            /* Waveform-specific sanity checks */
641
642
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
643
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
644
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
645
            if( !checkTidesZero(lambda1, lambda2) )
646
                ABORT_NONZERO_TIDES(LALparams);
647
648
649
650
651
652
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
            // NB: f_max = 0 will generate up to the ringdown cut-off frequency
            ret = XLALSimIMRPhenomBGenerateTD(hplus, hcross, phiRef, deltaT,
                    m1, m2, XLALSimIMRPhenomBComputeChi(m1, m2, S1z, S2z),
653
                    f_min, 0., distance, inclination);
654
            break;
655

656
657
        case PhenSpinTaylor:
            /* Waveform-specific sanity checks */
658

659
            /* Call the waveform driver routine */
660
661
662
	    XLALSimInspiralInitialConditionsPrecessingApproxs(&incl,&spin1x,&spin1y,&spin1z,&spin2x,&spin2y,&spin2z,inclination,S1x,S1y,S1z,S2x,S2y,S2z,m1,m2,f_ref,phiRef,XLALSimInspiralWaveformParamsLookupFrameAxis(LALparams));
	    polariz+=LAL_PI/2.;

663
            ret = XLALSimSpinInspiralGenerator(hplus, hcross, phiRef,
664
665
					       deltaT, m1, m2, f_min, f_ref, distance, incl, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
					       phaseO, amplitudeO, lambda1, lambda2, quadparam1, quadparam2, LALparams);
666
	    break;
667

668
669
        case IMRPhenomC:
            /* Waveform-specific sanity checks */
670
671
	  if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
672
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
673
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
674
            if( !checkTidesZero(lambda1, lambda2) )
675
                ABORT_NONZERO_TIDES(LALparams);
676
677
678
679
680
681
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
            // NB: f_max = 0 will generate up to the ringdown cut-off frequency
            ret = XLALSimIMRPhenomCGenerateTD(hplus, hcross, phiRef, deltaT,
                    m1, m2, XLALSimIMRPhenomBComputeChi(m1, m2, S1z, S2z),
682
                    f_min, 0., distance, inclination, LALparams);
683
            break;
684

685
	case IMRPhenomD:
686
687
	    if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
		    ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
688
	    if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
689
		    ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
690
	    if( !checkTidesZero(lambda1, lambda2) )
691
		    ABORT_NONZERO_TIDES(LALparams);
692
693
694
	    // generate TD waveforms with zero inclincation so that amplitude can be
	    // calculated from hplus and hcross, apply inclination-dependent factors
	    // in loop below
695
	    ret = XLALSimInspiralTDFromFD(hplus, hcross, m1, m2, S1x, S1y, S1z, S2x, S2y, S2z, distance, 0., phiRef, longAscNodes, eccentricity, meanPerAno, deltaT, f_min, f_ref, LALparams, approximant);
696
697
698
699
700
	    REAL8 maxamp=0;
	    REAL8TimeSeries *hp = *hplus;
	    REAL8TimeSeries *hc = *hcross;
	    INT4 maxind=hp->data->length - 1;
	    INT4 loopi;
701
	    const REAL8 cfac=cos(inclination);
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
	    const REAL8 pfac = 0.5 * (1. + cfac*cfac);

	    for (loopi=hp->data->length - 1; loopi > -1; loopi--)
	    {
		    REAL8 ampsqr = (hp->data->data[loopi])*(hp->data->data[loopi]) +
			   (hc->data->data[loopi])*(hc->data->data[loopi]);
		    if (ampsqr > maxamp)
		    {
			    maxind=loopi;
			    maxamp=ampsqr;
		    }
		    hp->data->data[loopi] *= pfac;
		    hc->data->data[loopi] *= cfac;
	    }
	    XLALGPSSetREAL8(&(hp->epoch), (-1.) * deltaT * maxind);
	    XLALGPSSetREAL8(&(hc->epoch), (-1.) * deltaT * maxind);
	    break;

720
	case IMRPhenomPv2:
721
	    ret = XLALSimInspiralTDFromFD(hplus, hcross, m1, m2, S1x, S1y, S1z, S2x, S2y, S2z, distance, inclination, phiRef, longAscNodes, eccentricity, meanPerAno, deltaT, f_min, f_ref, LALparams, approximant);
722
723
	    break;

724
725
726
        case PhenSpinTaylorRD:
            /* Waveform-specific sanity checks */
            if( !checkTidesZero(lambda1, lambda2) )
727
                ABORT_NONZERO_TIDES(LALparams);
728
729
730
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does use f_ref. The reference phase will be defined at the start.\n", __func__);
            /* Call the waveform driver routine */
731
732
	    XLALSimInspiralInitialConditionsPrecessingApproxs(&incl,&spin1x,&spin1y,&spin1z,&spin2x,&spin2y,&spin2z,inclination,S1x,S1y,S1z,S2x,S2y,S2z,m1,m2,f_ref,phiRef,XLALSimInspiralWaveformParamsLookupFrameAxis(LALparams));
	    polariz+=LAL_PI/2.;
733
            ret = XLALSimIMRPhenSpinInspiralRDGenerator(hplus, hcross, phiRef,
734
735
							deltaT, m1, m2, f_min, f_ref, distance, incl, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
							phaseO, amplitudeO, lambda1, lambda2, quadparam1, quadparam2, LALparams);
736
            break;
737

738
739
        case SEOBNRv1:
            /* Waveform-specific sanity checks */
740
741
	    if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
	        ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
742
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
743
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
744
            if( !checkTidesZero(lambda1, lambda2) )
745
	      ABORT_NONZERO_TIDES(LALparams);
746
747
748
749
750
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does not use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
            SpinAlignedEOBversion = 1;
            ret = XLALSimIMRSpinAlignedEOBWaveform(hplus, hcross, phiRef,
751
                    deltaT, m1, m2, f_min, distance, inclination, S1z, S2z, SpinAlignedEOBversion);
752
            break;
753

754
        case SEOBNRv2_opt:
755
756
        case SEOBNRv2:
            /* Waveform-specific sanity checks */
757
758
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
759
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
760
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
761
            if( !checkTidesZero(lambda1, lambda2) )
762
                ABORT_NONZERO_TIDES(LALparams);
763
764
765
766
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does not use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
            SpinAlignedEOBversion = 2;
767
            if(approximant==SEOBNRv2_opt) SpinAlignedEOBversion = 200;
768
            ret = XLALSimIMRSpinAlignedEOBWaveform(hplus, hcross, phiRef,
769
                    deltaT, m1, m2, f_min, distance, inclination, S1z, S2z, SpinAlignedEOBversion);
770
            break;
771

772
        case SEOBNRv4_opt:
773
774
        case SEOBNRv4:
            /* Waveform-specific sanity checks */
775
776
	    if(!XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams))
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
777
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
778
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
779
            if( !checkTidesZero(lambda1, lambda2) )
780
                ABORT_NONZERO_TIDES(LALparams);
781
782
783
784
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does not use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
            SpinAlignedEOBversion = 4;
785
            if(approximant==SEOBNRv4_opt) SpinAlignedEOBversion = 400;
786
            ret = XLALSimIMRSpinAlignedEOBWaveform(hplus, hcross, phiRef,
787
                                                   deltaT, m1, m2, f_min, distance, inclination, S1z, S2z, SpinAlignedEOBversion);
788
789
            break;

790
791
792
        case SEOBNRv3_opt_rk4:
        case SEOBNRv3_opt:
        case SEOBNRv3_pert:
793
794
        case SEOBNRv3:
            /* Waveform-specific sanity checks */
795
796
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
797
            if( !checkTidesZero(lambda1, lambda2) )
798
                ABORT_NONZERO_TIDES(LALparams);
799
800
801
            if( f_ref != 0.)
                XLALPrintWarning("XLAL Warning - %s: This approximant does use f_ref. The reference phase will be defined at coalescence.\n", __func__);
            /* Call the waveform driver routine */
802
803
804
805
806
	    //REAL8 spin1[3], spin2[3];
            spin1[0] = S1x; spin1[1] = S1y; spin1[2] = S1z;
            spin2[0] = S2x; spin2[1] = S2y; spin2[2] = S2z;
	    polariz+=-LAL_PI/2.;
            PrecEOBversion = 3;
807
808
809
810
811
812
            if(approximant==SEOBNRv3_pert) {
              const double m1pert = m1*(1.0 + 1e-15);
              ret = XLALSimIMRSpinEOBWaveform(hplus, hcross, /*&epoch,*/ phiRef,
                                              deltaT, m1pert, m2, f_min, distance, inclination, spin1, spin2, PrecEOBversion);

            } else {
813
814
              if(approximant==SEOBNRv3_opt) PrecEOBversion = 300;
              if(approximant==SEOBNRv3_opt_rk4) PrecEOBversion = 304;
815
816
817
              ret = XLALSimIMRSpinEOBWaveform(hplus, hcross, /*&epoch,*/ phiRef,
                                              deltaT, m1, m2, f_min, distance, inclination, spin1, spin2, PrecEOBversion);
            }
818
            break;
819

820
821
822
	case HGimri:
	     /* Waveform-specific sanity checks */
	     if( !checkTidesZero(lambda1, lambda2) )
823
		 ABORT_NONZERO_TIDES(LALparams);
824
	     if( !checkCOSpinZero(S2x, S2y, S2z) )
825
	         ABORT_NONZERO_SPIN2(LALparams);
826
	     /* Call the waveform driver */
827
	     ret = XLALHGimriGenerator(hplus, hcross, phiRef, deltaT, m1, m2, f_min, distance, inclination, S1z);
828
	     break;
829

830
831
832
833
        case NR_hdf5:
            /* Waveform-specific sanity checks */
            /* Call the waveform driver routine */
            ret = XLALSimInspiralNRWaveformGetHplusHcross(hplus, hcross,
834
835
                    phiRef, inclination, deltaT, m1, m2, distance, f_min, f_ref, S1x, S1y, S1z,
                    S2x, S2y, S2z, XLALSimInspiralWaveformParamsLookupNumRelData(LALparams));
836
837
838
            break;


839
840
841
842
        default:
            XLALPrintError("TD version of approximant not implemented in lalsimulation\n");
            XLAL_ERROR(XLAL_EINVAL);
    }
843

844
845
846
847
848
849
850
851
852
853
854
855
    if (polariz) {
      REAL8 tmpP,tmpC;
      REAL8 cp=cos(2.*polariz);
      REAL8 sp=sin(2.*polariz);
      for (UINT4 idx=0;idx<(*hplus)->data->length;idx++) {
	tmpP=(*hplus)->data->data[idx];
	tmpC=(*hcross)->data->data[idx];
	(*hplus)->data->data[idx] =cp*tmpP+sp*tmpC;
	(*hcross)->data->data[idx]=cp*tmpC-sp*tmpP;
      }
    }

856
    if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
857

858
    return ret;
859
860
}

861
/**
862
863
864
 * Chooses between different approximants when requesting a waveform to be generated
 * For spinning waveforms, all known spin effects up to given PN order are included
 * Returns the waveform in the frequency domain.
865
 */
866
867
868
int XLALSimInspiralChooseFDWaveform(
    COMPLEX16FrequencySeries **hptilde,     /**< FD plus polarization */
    COMPLEX16FrequencySeries **hctilde,     /**< FD cross polarization */
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
    const REAL8 m1,                         /**< mass of companion 1 (kg) */
    const REAL8 m2,                         /**< mass of companion 2 (kg) */
    const REAL8 S1x,                        /**< x-component of the dimensionless spin of object 1 */
    const REAL8 S1y,                        /**< y-component of the dimensionless spin of object 1 */
    const REAL8 S1z,                        /**< z-component of the dimensionless spin of object 1 */
    const REAL8 S2x,                        /**< x-component of the dimensionless spin of object 2 */
    const REAL8 S2y,                        /**< y-component of the dimensionless spin of object 2 */
    const REAL8 S2z,                        /**< z-component of the dimensionless spin of object 2 */
    const REAL8 distance,                   /**< distance of source (m) */
    const REAL8 inclination,                /**< inclination of source (rad) */
    const REAL8 phiRef,                     /**< reference orbital phase (rad) */
    const REAL8 longAscNodes,               /**< longitude of ascending nodes, degenerate with the polarization angle, Omega in documentation */
    const REAL8 eccentricity,               /**< eccentricity at reference epoch */
    const REAL8 UNUSED meanPerAno,          /**< mean anomaly of periastron */
    // frequency sampling parameters, no default value
    const REAL8 deltaF,                     /**< sampling interval (Hz) */
    const REAL8 f_min,                      /**< starting GW frequency (Hz) */
    const REAL8 f_max,                      /**< ending GW frequency (Hz) */
887
    REAL8 f_ref,                            /**< Reference frequency (Hz) */
888
889
    LALDict *LALparams,                     /**< LAL dictionary containing accessory parameters */
    const Approximant approximant           /**< post-Newtonian approximant to use for waveform production */
890
    )
891
{
892
893
894
895
896
897
898
899
900
    REAL8 LNhatx, LNhaty, LNhatz;
    REAL8 tmp1, tmp2;
    REAL8 E1x, E1y, E1z;
    REAL8 kMax;
    REAL8 v0, fStart;
    int ret;
    unsigned int j;
    REAL8 pfac, cfac;
    INT4 phiRefAtEnd;
901
902
903
904
905
906
    REAL8 quadparam1 = 1.+XLALSimInspiralWaveformParamsLookupdQuadMon1(LALparams);
    REAL8 quadparam2 = 1.+XLALSimInspiralWaveformParamsLookupdQuadMon2(LALparams);
    REAL8 lambda1 = XLALSimInspiralWaveformParamsLookupTidalLambda1(LALparams);
    REAL8 lambda2 = XLALSimInspiralWaveformParamsLookupTidalLambda2(LALparams);
    int amplitudeO = XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(LALparams);
    int phaseO =XLALSimInspiralWaveformParamsLookupPNPhaseOrder(LALparams);
907

908
    /* Support variables for precessing wfs*/
909
910
    REAL8 spin1x,spin1y,spin1z;
    REAL8 spin2x,spin2y,spin2z;
911

912
    /* Variables for IMRPhenomP and IMRPhenomPv2 */
913
914
    REAL8 chi1_l, chi2_l, chip, thetaJN, alpha0, phi_aligned, zeta_polariz;
    COMPLEX16 PhPpolp,PhPpolc;
915

916
917
918
919
920
    /* General sanity checks that will abort
     *
     * If non-GR approximants are added, include them in
     * XLALSimInspiralApproximantAcceptTestGRParams()
     */
921
    if( !XLALSimInspiralWaveformParamsNonGRAreDefault(LALparams) && XLALSimInspiralApproximantAcceptTestGRParams(approximant) != LAL_SIM_INSPIRAL_TESTGR_PARAMS ) {
922
923
924
        XLALPrintError("XLAL Error - %s: Passed in non-NULL pointer to LALSimInspiralTestGRParam for an approximant that does not use LALSimInspiralTestGRParam\n", __func__);
        XLAL_ERROR(XLAL_EINVAL);
    }
925

926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
    /* General sanity check the input parameters - only give warnings! */
    if( deltaF > 1. )
        XLALPrintWarning("XLAL Warning - %s: Large value of deltaF = %e requested...This corresponds to a very short TD signal (with padding). Consider a smaller value.\n", __func__, deltaF);
    if( deltaF < 1./4096. )
        XLALPrintWarning("XLAL Warning - %s: Small value of deltaF = %e requested...This corresponds to a very long TD signal. Consider a larger value.\n", __func__, deltaF);
    if( m1 < 0.09 * LAL_MSUN_SI )
        XLALPrintWarning("XLAL Warning - %s: Small value of m1 = %e (kg) = %e (Msun) requested...Perhaps you have a unit conversion error?\n", __func__, m1, m1/LAL_MSUN_SI);
    if( m2 < 0.09 * LAL_MSUN_SI )
        XLALPrintWarning("XLAL Warning - %s: Small value of m2 = %e (kg) = %e (Msun) requested...Perhaps you have a unit conversion error?\n", __func__, m2, m2/LAL_MSUN_SI);
    if( m1 + m2 > 1000. * LAL_MSUN_SI )
        XLALPrintWarning("XLAL Warning - %s: Large value of total mass m1+m2 = %e (kg) = %e (Msun) requested...Signal not likely to be in band of ground-based detectors.\n", __func__, m1+m2, (m1+m2)/LAL_MSUN_SI);
    if( S1x*S1x + S1y*S1y + S1z*S1z > 1.000001 )
        XLALPrintWarning("XLAL Warning - %s: S1 = (%e,%e,%e) with norm > 1 requested...Are you sure you want to violate the Kerr bound?\n", __func__, S1x, S1y, S1z);
    if( S2x*S2x + S2y*S2y + S2z*S2z > 1.000001 )
        XLALPrintWarning("XLAL Warning - %s: S2 = (%e,%e,%e) with norm > 1 requested...Are you sure you want to violate the Kerr bound?\n", __func__, S2x, S2y, S2z);
    if( f_min < 1. )
        XLALPrintWarning("XLAL Warning - %s: Small value of fmin = %e requested...Check for errors, this could create a very long waveform.\n", __func__, f_min);
    if( f_min > 40.000001 )
        XLALPrintWarning("XLAL Warning - %s: Large value of fmin = %e requested...Check for errors, the signal will start in band.\n", __func__, f_min);

    /* adjust the reference frequency for certain precessing approximants:
     * if that approximate interprets f_ref==0 to be f_min, set f_ref=f_min;
     * otherwise do nothing */
    f_ref = fixReferenceFrequency(f_ref, f_min, approximant);

    /* The non-precessing waveforms return h(f) for optimal orientation
     * (i=0, Fp=1, Fc=0; Lhat pointed toward the observer)
     * To get generic polarizations we multiply by inclination dependence
     * and note hc(f) \propto -I * hp(f)
     * Non-precessing waveforms multiply hp by pfac, hc by -I*cfac
956
     */
957
    cfac = cos(inclination);
958
    pfac = 0.5 * (1. + cfac*cfac);
959

960
    switch (approximant)
961
    {
962
963
964
        /* inspiral-only models */
	case EccentricFD:
            /* Waveform-specific sanity checks */
965
966
967
968
	    if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);
            if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALparams) )
                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
969
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
970
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
971
972
973
974
975
976
977
            /* Call the waveform driver routine */
            /* Note that for generic inclined eccentric waveforms
 *             * it is not possible to decompose hc(f) \propto I * hp(f)
 *                         * we call both polarizations independently
 *                                     */
            /*ret = XLALSimInspiralEFD(hptilde, hctilde, phiRef, deltaF, m1, m2,
 *             f_min, f_max, i, r, lambda1, lambda2, phaseO);*/
978
            // IMPORTANT CHECK: verify that inclination_azimuth is the longitude of ascending nodes
979
            ret = XLALSimInspiralEFD(hptilde, hctilde, phiRef, deltaF, m1, m2,
980
981
            f_min, f_max, inclination, distance,  longAscNodes,
             eccentricity, phaseO);
982
983
            if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
            break;
984

985
986
        case TaylorF2:
            /* Waveform-specific sanity checks */
987
988
989
990
            if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);
            if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALparams) )
                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
991
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
992
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
993

994
995
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorF2(hptilde, phiRef, deltaF, m1, m2,
996
                    S1z, S2z, f_min, f_max, f_ref, distance,
997
                    LALparams);
998
999
1000
            if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
            /* Produce both polarizations */
            *hctilde = XLALCreateCOMPLEX16FrequencySeries("FD hcross",
For faster browsing, not all history