LALSimInspiral.c

/*
 * Copyright (C) 2008 J. Creighton, S. Fairhurst, B. Krishnan, L. Santamaria, D. Keppel, Evan Ochsner, C. Pankow, 2014 A. Klein
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with with program; see the file COPYING. If not, write to the
 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
 *  MA  02111-1307  USA
 */

#include <complex.h>
#include <math.h>

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

#include <lal/SphericalHarmonics.h>
#include <lal/LALSimInspiral.h>
#include <lal/LALSimIMR.h>
#include <lal/LALSimSphHarmMode.h>
#include <lal/LALConstants.h>
#include <lal/LALStdlib.h>
#include <lal/LALString.h>
#include <lal/Sequence.h>
#include <lal/TimeSeries.h>
#include <lal/FrequencySeries.h>
#include <lal/TimeFreqFFT.h>
#include <lal/BandPassTimeSeries.h>
#include <lal/Units.h>
#include <lal/LALSimBlackHoleRingdown.h>

#include "LALSimInspiralPNCoefficients.c"
#include "check_series_macros.h"
#include "check_waveform_macros.h"

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

/**
 * (Twice) the highest known PN order of amplitude correction for
 * non-precessing binaries.
 */
#define MAX_NONPRECESSING_AMP_PN_ORDER 6

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

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

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


#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),
    INITIALIZE_NAME(EOBNRv2_ROM),
    INITIALIZE_NAME(EOBNRv2HM_ROM),
    INITIALIZE_NAME(SEOBNRv1),
    INITIALIZE_NAME(SEOBNRv2),
    INITIALIZE_NAME(SEOBNRv2_opt),
    INITIALIZE_NAME(SEOBNRv3),
    INITIALIZE_NAME(SEOBNRv3_pert),
    INITIALIZE_NAME(SEOBNRv3_opt),
    INITIALIZE_NAME(SEOBNRv3_opt_rk4),
    INITIALIZE_NAME(SEOBNRv4),
    INITIALIZE_NAME(SEOBNRv4_opt),
    INITIALIZE_NAME(SEOBNRv1_ROM_EffectiveSpin),
    INITIALIZE_NAME(SEOBNRv1_ROM_DoubleSpin),
    INITIALIZE_NAME(SEOBNRv2_ROM_EffectiveSpin),
    INITIALIZE_NAME(SEOBNRv2_ROM_DoubleSpin),
    INITIALIZE_NAME(SEOBNRv2_ROM_DoubleSpin_HI),
    INITIALIZE_NAME(Lackey_Tidal_2013_SEOBNRv2_ROM),
    INITIALIZE_NAME(SEOBNRv4_ROM),
    INITIALIZE_NAME(HGimri),
    INITIALIZE_NAME(IMRPhenomA),
    INITIALIZE_NAME(IMRPhenomB),
    INITIALIZE_NAME(IMRPhenomFA),
    INITIALIZE_NAME(IMRPhenomFB),
    INITIALIZE_NAME(IMRPhenomC),
    INITIALIZE_NAME(IMRPhenomD),
    INITIALIZE_NAME(IMRPhenomP),
    INITIALIZE_NAME(IMRPhenomPv2),
    INITIALIZE_NAME(IMRPhenomFC),
    INITIALIZE_NAME(TaylorEt),
    INITIALIZE_NAME(TaylorT4),
    INITIALIZE_NAME(EccentricTD),
    INITIALIZE_NAME(TaylorN),
    INITIALIZE_NAME(SpinTaylorT4Fourier),
    INITIALIZE_NAME(SpinTaylorT2Fourier),
    INITIALIZE_NAME(SpinDominatedWf),
    INITIALIZE_NAME(NR_hdf5),
};
#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;
}


/*
 * 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.
 */
static double fixReferenceFrequency(const double f_ref, const double f_min, const Approximant approximant)
{
    if (f_ref == 0)
        switch (approximant) {
        case SpinTaylorT1:
        case SpinTaylorT2:
        case SpinTaylorT3:
        case SpinTaylorT4:
        case SpinTaylorT2Fourier:
        case SpinTaylorT4Fourier:
        case SpinTaylorF2:
        case IMRPhenomP:
        case IMRPhenomPv2:
            return f_min;
        default:
            break;
        }
    return f_ref;
}

/*
 * some helper routines for XLALSimInspiralTD
 */
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);

/**
 * @addtogroup LALSimInspiral_c
 * @brief General routines for generating binary inspiral waveforms.
 *
 * @{
 */

/**
 * @name General Waveform Switching Generation Routines
 * @{
 */

/**
 * 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.
 */
int XLALSimInspiralChooseTDWaveform(
    REAL8TimeSeries **hplus,                    /**< +-polarization waveform */
    REAL8TimeSeries **hcross,                   /**< x-polarization waveform */
    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) */
    REAL8 f_ref,                                /**< reference GW frequency (Hz) */
    LALDict *LALparams,                         /**< LAL dictionary containing accessory parameters */
    const Approximant approximant               /**< post-Newtonian approximant to use for waveform production */
    )
{
    REAL8 LNhatx, LNhaty, LNhatz, E1x, E1y, E1z;
    //REAL8 tmp1, tmp2;
    int ret;
    /* N.B. the quadrupole of a spinning compact body labeled by A is
     * 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.
     */
    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);

    /* General sanity checks that will abort
     *
     * If non-GR approximants are added, include them in
     * XLALSimInspiralApproximantAcceptTestGRParams()
     */
    if( !XLALSimInspiralWaveformParamsNonGRAreDefault(LALparams) && XLALSimInspiralApproximantAcceptTestGRParams(approximant) != LAL_SIM_INSPIRAL_TESTGR_PARAMS ) {
        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);
    }
    /* Support variables for precessing wfs*/
    REAL8 incl;


    /* SEOBNR flag for spin aligned model version. 1 for SEOBNRv1, 2 for SEOBNRv2 */
    UINT4 SpinAlignedEOBversion;
    REAL8 spin1x,spin1y,spin1z;
    REAL8 spin2x,spin2y,spin2z;
    REAL8 polariz=longAscNodes;

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

    //LIGOTimeGPS epoch = LIGOTIMEGPSZERO;

    /* 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);

    /* 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);

    switch (approximant)
    {
        /* non-spinning inspiral-only models */
        case TaylorEt:
            /* Waveform-specific sanity checks */
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
                ABORT_NONZERO_SPINS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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,
                    deltaT, m1, m2, f_min, distance, inclination, amplitudeO, phaseO);
            break;

        case TaylorT1:
            /* Waveform-specific sanity checks */
	    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);
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorT1PNGenerator(hplus, hcross, phiRef, v0,
                    deltaT, m1, m2, f_min, f_ref, distance, inclination, lambda1, lambda2,
                    XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALparams), amplitudeO, phaseO);
            break;

        case TaylorT2:
            /* Waveform-specific sanity checks */
	    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);
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorT2PNGenerator(hplus, hcross, phiRef, v0,
                    deltaT, m1, m2, f_min, f_ref, distance, inclination, lambda1, lambda2,
                    XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALparams), amplitudeO, phaseO);
            break;

        case TaylorT3:
            /* Waveform-specific sanity checks */
	    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);
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorT3PNGenerator(hplus, hcross, phiRef, v0,
                    deltaT, m1, m2, f_min, f_ref, distance, inclination, lambda1, lambda2,
                    XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALparams), amplitudeO, phaseO);
            break;

        case TaylorT4:
            /* Waveform-specific sanity checks */
	    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);
            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorT4PNGenerator(hplus, hcross, phiRef, v0,
                    deltaT, m1, m2, f_min, f_ref, distance, inclination, lambda1, lambda2,
                    XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALparams), amplitudeO, phaseO);
            break;

	case EccentricTD:
	    /* Waveform-specific sanity checks */
	    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);
	    if( !checkTidesZero(lambda1, lambda2) )
		ABORT_NONZERO_TIDES(LALparams);
	    /* Call the waveform driver routine */
	    ret = XLALSimInspiralEccentricTDPNGenerator(hplus, hcross, phiRef,
		    deltaT, m1, m2, f_min, f_ref, distance, inclination, eccentricity,
		    amplitudeO, phaseO);
	    if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
	    break;

        /* non-spinning inspiral-merger-ringdown models */
        case IMRPhenomA:
            /* Waveform-specific sanity checks */
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
                ABORT_NONZERO_SPINS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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,
                    m1, m2, f_min, 0., distance, inclination);
            break;

        case EOBNRv2HM:
            /* Waveform-specific sanity checks */
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
	      ABORT_NONZERO_SPINS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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,
                    m1, m2, f_min, distance, inclination);
            break;

        case EOBNRv2:
            /* Waveform-specific sanity checks */
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkSpinsZero(S1x, S1y, S1z, S2x, S2y, S2z) )
                ABORT_NONZERO_SPINS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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,
                    m1, m2, f_min, distance, inclination);
            break;

        /* spinning inspiral-only models */
        case SpinTaylorT2:
            /* Waveform-specific sanity checks */
	    if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) && (XLALSimInspiralWaveformParamsLookupPNSpinOrder(LALparams)>5) )
	      ABORT_NONZERO_TRANSVERSE_SPINS_HIGH_SPINO(LALparams);
	    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);
            LNhaty = 0.;
            LNhatz = cos(incl);
            E1x = 0.;
            E1y = 1.;
            E1z = 0.;
	    polariz+=LAL_PI/2.;
            /* 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 = XLALSimInspiralSpinTaylorT2(hplus, hcross, phiRef, v0, deltaT,
					      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);
            break;

        // need to make a consistent choice for SpinTaylorT4 and PSpinInspiralRD waveform inputs
        // proposal: TotalJ frame of PSpinInspiralRD
        // 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
        // initial ** L **.
        case SpinTaylorT4:
            /* Waveform-specific sanity checks */
	    if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) && (XLALSimInspiralWaveformParamsLookupPNSpinOrder(LALparams)>5) )
	      ABORT_NONZERO_TRANSVERSE_SPINS_HIGH_SPINO(LALparams);
            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);
            LNhaty = 0.;
            LNhatz = cos(incl);
            E1x = 0.;
            E1y = 1.;
            E1z = 0.;
	    polariz+=LAL_PI/2.;
            /* 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 = XLALSimInspiralSpinTaylorT4(hplus, hcross, phiRef, v0, deltaT,
                    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);
            break;

	 case SpinTaylorT1:
            /* Waveform-specific sanity checks */
            /* Waveform-specific sanity checks */
	    if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) && (XLALSimInspiralWaveformParamsLookupPNSpinOrder(LALparams)>5) )
	      ABORT_NONZERO_TRANSVERSE_SPINS_HIGH_SPINO(LALparams);
	    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);
            LNhaty = 0.;
            LNhatz = cos(incl);
            E1x = 0.;
            E1y = 1.;
            E1z = 0.;
	    polariz+=LAL_PI/2.;
            /* 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,
					      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);
            break;

        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);
                }
	            incl=inclination;
                LNhatx = 0.;
                LNhaty = 0.;
                LNhatz = 1.;
                /* Call the waveform driver routine */
                ret = XLALSimInspiralSpinDominatedWaveformInterfaceTD(hplus, hcross, deltaT, m1, m2, f_min, f_ref, distance, S1x, S1y, S1z, LNhatx, LNhaty, LNhatz, incl, phaseO, amplitudeO, phiRef);
                break;

        /* spin aligned inspiral-merger-ringdown models */
        case IMRPhenomB:
            /* Waveform-specific sanity checks */
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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),
                    f_min, 0., distance, inclination);
            break;

        case PhenSpinTaylor:
            /* Waveform-specific sanity checks */

            /* Call the waveform driver routine */
	    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.;

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

        case IMRPhenomC:
            /* Waveform-specific sanity checks */
	  if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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),
                    f_min, 0., distance, inclination, LALparams);
            break;

	case IMRPhenomD:
	    if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
		    ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
	    if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
		    ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
	    if( !checkTidesZero(lambda1, lambda2) )
		    ABORT_NONZERO_TIDES(LALparams);
	    // generate TD waveforms with zero inclincation so that amplitude can be
	    // calculated from hplus and hcross, apply inclination-dependent factors
	    // in loop below
	    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);
	    REAL8 maxamp=0;
	    REAL8TimeSeries *hp = *hplus;
	    REAL8TimeSeries *hc = *hcross;
	    INT4 maxind=hp->data->length - 1;
	    INT4 loopi;
	    const REAL8 cfac=cos(inclination);
	    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;

	case IMRPhenomPv2:
	    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);
	    break;

        case PhenSpinTaylorRD:
            /* Waveform-specific sanity checks */
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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 */
	    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.;
            ret = XLALSimIMRPhenSpinInspiralRDGenerator(hplus, hcross, phiRef,
							deltaT, m1, m2, f_min, f_ref, distance, incl, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
							phaseO, amplitudeO, lambda1, lambda2, quadparam1, quadparam2, LALparams);
            break;

        case SEOBNRv1:
            /* Waveform-specific sanity checks */
	    if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
	        ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
	      ABORT_NONZERO_TIDES(LALparams);
            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,
                    deltaT, m1, m2, f_min, distance, inclination, S1z, S2z, SpinAlignedEOBversion);
            break;

        case SEOBNRv2_opt:
        case SEOBNRv2:
            /* Waveform-specific sanity checks */
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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;
            if(approximant==SEOBNRv2_opt) SpinAlignedEOBversion = 200;
            ret = XLALSimIMRSpinAlignedEOBWaveform(hplus, hcross, phiRef,
                    deltaT, m1, m2, f_min, distance, inclination, S1z, S2z, SpinAlignedEOBversion);
            break;

        case SEOBNRv4_opt:
        case SEOBNRv4:
            /* Waveform-specific sanity checks */
	    if(!XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams))
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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;
            if(approximant==SEOBNRv4_opt) SpinAlignedEOBversion = 400;
            ret = XLALSimIMRSpinAlignedEOBWaveform(hplus, hcross, phiRef,
                                                   deltaT, m1, m2, f_min, distance, inclination, S1z, S2z, SpinAlignedEOBversion);
            break;

        case SEOBNRv3_opt_rk4:
        case SEOBNRv3_opt:
        case SEOBNRv3_pert:
        case SEOBNRv3:
            /* Waveform-specific sanity checks */
            if( !XLALSimInspiralWaveformParamsFlagsAreDefault(LALparams) )
                ABORT_NONDEFAULT_LALDICT_FLAGS(LALparams);
            if( !checkTidesZero(lambda1, lambda2) )
                ABORT_NONZERO_TIDES(LALparams);
            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 */
	    //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;
            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 {
              if(approximant==SEOBNRv3_opt) PrecEOBversion = 300;
              if(approximant==SEOBNRv3_opt_rk4) PrecEOBversion = 304;
              ret = XLALSimIMRSpinEOBWaveform(hplus, hcross, /*&epoch,*/ phiRef,
                                              deltaT, m1, m2, f_min, distance, inclination, spin1, spin2, PrecEOBversion);
            }
            break;

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

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


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

    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;
      }
    }

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

    return ret;
}

/**
 * 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.
 */
int XLALSimInspiralChooseFDWaveform(
    COMPLEX16FrequencySeries **hptilde,     /**< FD plus polarization */
    COMPLEX16FrequencySeries **hctilde,     /**< FD cross polarization */
    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) */
    REAL8 f_ref,                            /**< Reference frequency (Hz) */
    LALDict *LALparams,                     /**< LAL dictionary containing accessory parameters */
    const Approximant approximant           /**< post-Newtonian approximant to use for waveform production */
    )
{
    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;
    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);

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

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

    /* General sanity checks that will abort
     *
     * If non-GR approximants are added, include them in
     * XLALSimInspiralApproximantAcceptTestGRParams()
     */
    if( !XLALSimInspiralWaveformParamsNonGRAreDefault(LALparams) && XLALSimInspiralApproximantAcceptTestGRParams(approximant) != LAL_SIM_INSPIRAL_TESTGR_PARAMS ) {
        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);
    }

    /* 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
     */
    cfac = cos(inclination);
    pfac = 0.5 * (1. + cfac*cfac);

    switch (approximant)
    {
        /* inspiral-only models */
	case EccentricFD:
            /* Waveform-specific sanity checks */
	    if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);
            if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALparams) )
                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);
            /* 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);*/
            // IMPORTANT CHECK: verify that inclination_azimuth is the longitude of ascending nodes
            ret = XLALSimInspiralEFD(hptilde, hctilde, phiRef, deltaF, m1, m2,
            f_min, f_max, inclination, distance,  longAscNodes,
             eccentricity, phaseO);
            if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
            break;

        case TaylorF2:
            /* Waveform-specific sanity checks */
            if( !XLALSimInspiralWaveformParamsFrameAxisIsDefault(LALparams) )
                ABORT_NONDEFAULT_FRAME_AXIS(LALparams);
            if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALparams) )
                ABORT_NONDEFAULT_MODES_CHOICE(LALparams);
            if( !checkTransverseSpinsZero(S1x, S1y, S2x, S2y) )
                ABORT_NONZERO_TRANSVERSE_SPINS(LALparams);

            /* Call the waveform driver routine */
            ret = XLALSimInspiralTaylorF2(hptilde, phiRef, deltaF, m1, m2,
                    S1z, S2z, f_min, f_max, f_ref, distance,
                    LALparams);
            if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
            /* Produce both polarizations */
            *hctilde = XLALCreateCOMPLEX16FrequencySeries("FD hcross",