LALSimIMRPSpinInspiralRD.c

/*
 * Copyright (C) 2011 Riccardo Sturani, John Veitch, Drew Keppel
 *
 *  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 <stdlib.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <lal/LALStdlib.h>
#include <lal/AVFactories.h>
#include <lal/SeqFactories.h>
#include <lal/Units.h>
#include <lal/TimeSeries.h>
#include <lal/LALConstants.h>
#include <lal/SeqFactories.h>
#include <lal/RealFFT.h>
#include <lal/SphericalHarmonics.h>
#include <lal/LALAdaptiveRungeKutta4.h>
#include <lal/LALSimInspiral.h>
#include <lal/LALSimIMR.h>
#include <lal/LALSimSphHarmMode.h>
#include <lal/Date.h>

//#include "LALSimIMRPhenSpin.h"
#include "LALSimPhenSpinRingDown.c"
#include "LALSimInspiralPNCoefficients.c"

#define LALSIMINSPIRAL_PST4_TEST_ENERGY 		1025
#define LALSIMINSPIRAL_PST4_TEST_OMEGADOT 		1026
#define LALSIMINSPIRAL_PST4_TEST_COORDINATE 		1027
#define LALSIMINSPIRAL_PST4_TEST_OMEGANAN 		1028
#define LALSIMINSPIRAL_PST4_TEST_FREQBOUND 		1029
#define LALSIMINSPIRAL_PST4_DERIVATIVE_OMEGANONPOS 	1030
#define LALSIMINSPIRAL_PST4_TEST_OMEGAMATCH             1031

#define LAL_NUM_PST4_VARIABLES 12

#define LAL_PST4_ABSOLUTE_TOLERANCE 1.e-12
#define LAL_PST4_RELATIVE_TOLERANCE 1.e-12

/* Minimum integration length */
#define minIntLen    16
/* For turning on debugging messages*/
#define DEBUG_RD  0

#define nModes 8
#define RD_EFOLDS 10

static REAL8 dsign(INT4 i){
  if (i>=0) return 1.;
  else return -1.;
}

typedef struct tagLALSimInspiralPhenSpinTaylorT4Coeffs {
  REAL8 M; ///< total mass in seconds
  REAL8 eta; ///< symmetric mass ratio
  REAL8 m1ByM; ///< m1 / M
  REAL8 m2ByM; ///< m2 / M
  REAL8 dt; ///< UNDOCUMENTED
  REAL8 wdotnewt; ///< leading order coefficient of wdot = \f$\dot{\omega}\f$
  REAL8 wdotcoeff[LAL_MAX_PN_ORDER]; ///< coeffs. of PN corrections to wdot
  REAL8 wdotlogcoeff; ///< coefficient of log term in wdot
  REAL8 Enewt; ///< coeffs. of PN corrections to energy
  REAL8 Ecoeff[LAL_MAX_PN_ORDER]; ///< coeffs. of PN corrections to energy
  REAL8 wdot3S1O, wdot3S2O; ///< non-dynamical 1.5PN SO corrections
  REAL8 wdot4S1S2; ///< non-dynamical 2PN SS correction
  REAL8 wdot4S1OS2O; ///< non-dynamical 2PN SS correction
  REAL8 wdot4S1S1,wdot4S2S2; ///< non-dynamical 2PN SS correction
  REAL8 wdot4S1OS1O,wdot4S2OS2O; ///< non-dynamical 2PN SS correction
  REAL8 wdot4QMS1; ///< non-dynamical S1^2 2PN quadrupole-monopole correction
  REAL8 wdot4QMS1O; ///< non-dynamical (S1.L)^2 2PN quadrupole-monopole correction
  REAL8 wdot4QMS2; ///< non-dynamical S2^2 2PN quadrupole-monopole correction
  REAL8 wdot4QMS2O; ///< non-dynamical (S2.L)^2 2PN quadrupole-monopole correction
  REAL8 wdot5S1O, wdot5S2O; ///< non-dynamical 2.5PN SO corrections
  REAL8 wdot6S1O, wdot6S2O; ///< non-dynamical 3PN SO corrections
  REAL8 E3S1O, E3S2O; ///< non-dynamical 1.5PN SO corrections
  REAL8 E4S1S2,E4S1OS2O; ///< non-dynamical 2PN SS correction
  REAL8 E4QMS1; ///< non-dynamical S1^2 2PN quadrupole-monopole correction
  REAL8 E4QMS1O;///< non-dynamical (S1.L)^2 2PN quadrupole-monopole correction
  REAL8 E4QMS2; ///< non-dynamical S2^2 2PN quadrupole-monopole correction
  REAL8 E4QMS2O;///< non-dynamical (S2.L)^2 2PN quadrupole-monopole correction
  REAL8 E5S1O, E5S2O; ///< non-dynamical 2.5PN SO corrections
  REAL8 E6S1O, E6S2O; ///< non-dynamical 2.5PN SO corrections
  REAL8 LNhat3S1O, LNhat3S2O; ///< non-dynamical 1.5PN SO corrections
  REAL8 LNhat4SS; ///< non-dynamical 2PN SS correction
  REAL8 wdottidal5pn;	///< leading order tidal correction
  REAL8 wdottidal6pn;	///< next to leading order tidal correction
  REAL8 Etidal5pn; ///< leading order tidal correction to energy
  REAL8 Etidal6pn; ///< next to leading order tidal correction to energy
  REAL8 S1dot3, S2dot3;  ///< UNDOCUMENTED
  REAL8 Sdot4S2,Sdot4S2O;  ///< UNDOCUMENTED
  REAL8 S1dot4QMS1O,S2dot4QMS2O;  ///< UNDOCUMENTED
  REAL8 S1dot5, S2dot5;  ///< UNDOCUMENTED
  REAL8 S1dot6, S2dot6;  ///< UNDOCUMENTED
  REAL8 fStart; ///< starting GW frequency of integration
  REAL8 fEnd; ///< ending GW frequency of integration
} LALSimInspiralPhenSpinTaylorT4Coeffs;

static REAL8 OmMatch(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2) {

  const REAL8 omM       = 0.05;
  const REAL8 omMsz12   =    9.97e-4;
  const REAL8 omMs1d2   =  -2.032e-3;
  const REAL8 omMssq    =   5.629e-3;
  const REAL8 omMsz1d2  =   8.646e-3;
  const REAL8 omMszsq   =  -5.909e-3;
  const REAL8 omM3s1d2  =   1.801e-3;
  const REAL8 omM3ssq   = -1.4059e-2;
  const REAL8 omM3sz1d2 =  1.5483e-2;
  const REAL8 omM3szsq  =   8.922e-3;

  return omM + /*6.05e-3 * sqrtOneMinus4Eta +*/
    omMsz12   * (LNhS1 + LNhS2) +
    omMs1d2   * (S1S2) +
    omMssq    * (S1S1 + S2S2) +
    omMsz1d2  * (LNhS1 * LNhS2) +
    omMszsq   * (LNhS1 * LNhS1 + LNhS2 * LNhS2) +
    omM3s1d2  * (LNhS1 + LNhS2) * (S1S2) +
    omM3ssq   * (LNhS1 + LNhS2) * (S1S1+S2S2) +
    omM3sz1d2 * (LNhS1 + LNhS2) * (LNhS1*LNhS2) +
    omM3szsq  * (LNhS1 + LNhS2) * (LNhS1*LNhS1+LNhS2*LNhS2);
} /* End of OmMatch */

static REAL8 fracRD(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2) {

  const double frac0      = 0.840;
  const double fracsz12   = -2.145e-2;
  const double fracs1d2   = -4.421e-2;
  const double fracssq    = -2.643e-2;
  const double fracsz1d2  = -5.876e-2;
  const double fracszsq   = -2.215e-2;

  return frac0 +
    fracsz12   * (LNhS1 + LNhS2) +
    fracs1d2   * (S1S2) +
    fracssq    * (S1S1 + S2S2) +
    fracsz1d2  * (LNhS1 * LNhS2) +
    fracszsq   * (LNhS1 * LNhS1 + LNhS2 * LNhS2);
} /* End of fracRD */

/**
 * Convenience function to set up XLALSimInspiralSpinTaylotT4Coeffs struct
 */
static INT4 XLALSimIMRPhenSpinParamsSetup(LALSimInspiralPhenSpinTaylorT4Coeffs  *params, /**< PN params [returned] */
                                         REAL8 dt,                                      /**< Sampling in secs */
                                         REAL8 fStart,                                  /**< Starting frequency of integration*/
                                         REAL8 fEnd,                                    /**< Ending frequency of integration*/
                                         REAL8 mass1,                                   /**< Mass 1 in solar mass units */
                                         REAL8 mass2,                                   /**< Mass 2 in solar mass units */
					 REAL8 lambda1,      /**< Tidal par1*/
                                         REAL8 lambda2,      /**< Tidal par2*/
                                         REAL8 quadparam1,   /**< Quad-monop par1*/
					 REAL8 quadparam2,   /**< Quad-monop par2*/
					 int spinO,                 /**< Spin interaction */
					 int tideO,                /**< Tidal iteraction interaction */
					 LALDict *LALparams,    /**< Extra parameters */
                                         UINT4 order                                    /**< twice PN Order in Phase */
)
{
  /* Zero the coefficients */
  memset(params, 0, sizeof(LALSimInspiralPhenSpinTaylorT4Coeffs));
  params->eta    = mass1*mass2/(mass1+mass2)/(mass1+mass2);
  params->m1ByM  = mass1 / (mass1+mass2);
  params->m2ByM  = mass2 / (mass1+mass2);
  params->M     = (mass1+mass2)*LAL_MTSUN_SI;
  REAL8 unitHz   = params->M *((REAL8) LAL_PI);

  params->fEnd   = fEnd*unitHz;    /*On the left side there is actually an omega*/
  params->fStart = fStart*unitHz;  /*On the left side there is actually an omega*/
  if (fEnd>0.)
    params->dt     = dt*(fEnd-fStart)/fabs(fEnd-fStart);
  else
    params->dt     = dt;

  REAL8 phi1 = XLALSimInspiralWaveformParamsLookupNonGRPhi1(LALparams);

  params->wdotnewt = XLALSimInspiralTaylorT4wdot_0PNCoeff(params->eta);
  params->Enewt    = XLALSimInspiralPNEnergy_0PNCoeff(params->eta);

  switch (order) {

    case -1: // Use the highest PN order available.
    case 7:
      params->wdotcoeff[7]  = XLALSimInspiralTaylorT4wdot_7PNCoeff(params->eta);

    case 6:
      params->Ecoeff[6]     = XLALSimInspiralPNEnergy_6PNCoeff(params->eta);
      params->wdotcoeff[6]  = XLALSimInspiralTaylorT4wdot_6PNCoeff(params->eta);
      params->wdotlogcoeff  = XLALSimInspiralTaylorT4wdot_6PNLogCoeff(params->eta);
      params->wdot6S1O   = XLALSimInspiralTaylorT4wdot_6PNSOCoeff(params->m1ByM);
      params->wdot6S2O   = XLALSimInspiralTaylorT4wdot_6PNSOCoeff(params->m2ByM);

    case 5:
      params->Ecoeff[5]     = 0.;
      params->wdotcoeff[5]  = XLALSimInspiralTaylorT4wdot_5PNCoeff(params->eta);
      params->E5S1O     = XLALSimInspiralPNEnergy_5PNSOCoeff(params->m1ByM);
      params->E5S2O     = XLALSimInspiralPNEnergy_5PNSOCoeff(params->m1ByM);
      params->wdot5S1O  = XLALSimInspiralTaylorT4wdot_5PNSOCoeff(params->m1ByM);
      params->wdot5S2O  = XLALSimInspiralTaylorT4wdot_5PNSOCoeff(params->m2ByM);
      params->S1dot5    = XLALSimInspiralSpinDot_5PNCoeff(params->m1ByM);
      params->S2dot5    = XLALSimInspiralSpinDot_5PNCoeff(params->m2ByM);

    case 4:
      params->wdotcoeff[4]  = XLALSimInspiralTaylorT4wdot_4PNCoeff(params->eta);
      params->Ecoeff[4]     = XLALSimInspiralPNEnergy_4PNCoeff(params->eta);
      params->wdot4S1S2     = XLALSimInspiralTaylorT4wdot_4PNS1S2CoeffAvg(params->eta);
      params->wdot4S1OS2O   = XLALSimInspiralTaylorT4wdot_4PNS1S2OCoeffAvg(params->eta);
      params->E4S1S2      = XLALSimInspiralPNEnergy_4PNS1S2CoeffAvg(params->eta);
      params->E4S1OS2O    = XLALSimInspiralPNEnergy_4PNS1S2OCoeffAvg(params->eta);
      params->wdot4S1S1     = XLALSimInspiralTaylorT4wdot_4PNSelf2SCoeff(params->m1ByM);
      params->wdot4S1OS1O   = XLALSimInspiralTaylorT4wdot_4PNSelf2SOCoeff(params->m1ByM);
      params->wdot4S2S2     = XLALSimInspiralTaylorT4wdot_4PNSelf2SCoeff(params->m2ByM);
      params->wdot4S2OS2O   = XLALSimInspiralTaylorT4wdot_4PNSelf2SOCoeff(params->m2ByM);
      params->wdot4QMS1     = quadparam1*XLALSimInspiralTaylorT4wdot_4PNQM2SCoeff(params->m1ByM);
      params->wdot4QMS1O    = quadparam1*XLALSimInspiralTaylorT4wdot_4PNQM2SOCoeff(params->m1ByM);
      params->wdot4QMS2     = quadparam2*XLALSimInspiralTaylorT4wdot_4PNQM2SCoeff(params->m2ByM);
      params->wdot4QMS2O    = quadparam2*XLALSimInspiralTaylorT4wdot_4PNQM2SOCoeff(params->m2ByM);
      params->E4QMS1      = quadparam1*XLALSimInspiralPNEnergy_4PNQM2SCoeffAvg(params->m1ByM);
      params->E4QMS2      = quadparam2*XLALSimInspiralPNEnergy_4PNQM2SCoeffAvg(params->m2ByM);
      params->E4QMS1O     = quadparam1*XLALSimInspiralPNEnergy_4PNQM2SOCoeffAvg(params->m1ByM);
      params->E4QMS2O     = quadparam2*XLALSimInspiralPNEnergy_4PNQM2SOCoeffAvg(params->m2ByM);
      params->Sdot4S2     = XLALSimInspiralSpinDot_4PNS2CoeffAvg;
      params->Sdot4S2O    = XLALSimInspiralSpinDot_4PNS2OCoeffAvg;
      params->S1dot4QMS1O  = quadparam1*XLALSimInspiralSpinDot_4PNQMSOCoeffAvg(params->m1ByM);
      params->S2dot4QMS2O  = quadparam2*XLALSimInspiralSpinDot_4PNQMSOCoeffAvg(params->m2ByM);

    case 3:
      params->Ecoeff[3]      = 0.;
      params->wdotcoeff[3]   = XLALSimInspiralTaylorT4wdot_3PNCoeff(params->eta);
      params->wdot3S1O  = XLALSimInspiralTaylorT4wdot_3PNSOCoeff(params->m1ByM);
      params->wdot3S2O  = XLALSimInspiralTaylorT4wdot_3PNSOCoeff(params->m1ByM);
      params->E3S1O     = XLALSimInspiralPNEnergy_3PNSOCoeff(params->m1ByM);
      params->E3S2O     = XLALSimInspiralPNEnergy_3PNSOCoeff(params->m1ByM);
      params->S1dot3    = XLALSimInspiralSpinDot_3PNCoeff(params->m1ByM);
      params->S2dot3    = XLALSimInspiralSpinDot_3PNCoeff(params->m2ByM);

    case 2:
      params->Ecoeff[2]  = XLALSimInspiralPNEnergy_2PNCoeff(params->eta);
      params->wdotcoeff[2] = XLALSimInspiralTaylorT4wdot_2PNCoeff(params->eta);

    case 1:
      params->Ecoeff[1]  = 0.;
      params->wdotcoeff[1] = phi1;

    case 0:
      params->Ecoeff[0]  = 1.;
      params->wdotcoeff[0] = 1.;
      break;

    case 8:
      XLALPrintError("*** LALSimIMRPhenSpinInspiralRD ERROR: PhenSpin approximant not available at pseudo4PN order\n");
                        XLAL_ERROR(XLAL_EDOM);
      break;

    default:
      XLALPrintError("*** LALSimIMRPhenSpinInspiralRD ERROR: Impossible to create waveform with %d order\n",order);
                        XLAL_ERROR(XLAL_EFAILED);
      break;
  }

  switch (spinO) {

    case LAL_SIM_INSPIRAL_SPIN_ORDER_0PN:
    case LAL_SIM_INSPIRAL_SPIN_ORDER_05PN:
    case LAL_SIM_INSPIRAL_SPIN_ORDER_1PN:
      /*This kills all spin effects in the phase. Still there are spin effects
        in the waveform due to orbital plane precession*/
      params->E6S1O      = 0.;
      params->E6S2O      = 0.;
      params->wdot6S1O   = 0.;
      params->wdot6S1O   = 0.;
      params->S1dot6     = 0.;
      params->S2dot6     = 0.;

    case LAL_SIM_INSPIRAL_SPIN_ORDER_15PN:
      /* This keeps only the leading spin-orbit interactions*/
      params->wdot4S1S2   = 0.;
      params->wdot4S1OS2O = 0.;
      params->E4QMS1  = 0.;
      params->E4QMS2  = 0.;
      params->E4QMS1O = 0.;
      params->E4QMS2O = 0.;

    case LAL_SIM_INSPIRAL_SPIN_ORDER_2PN:
      /* This kills all spin interaction intervening at 2.5PN order or higher*/
      params->E5S1O       = 0.;
      params->E5S2O       = 0.;
      params->wdot5S1O    = 0.;
      params->wdot5S2O    = 0.;
      params->S1dot5      = 0.;
      params->S2dot5      = 0.;

    case LAL_SIM_INSPIRAL_SPIN_ORDER_25PN:
      /* This kills all spin interaction intervening at 3PN order or higher*/
      params->wdot6S1O   = 0.;
      params->wdot6S2O   = 0.;

    case LAL_SIM_INSPIRAL_SPIN_ORDER_3PN:
    case LAL_SIM_INSPIRAL_SPIN_ORDER_ALL:
    default:
      break;
  }

  switch( tideO ) {
    case LAL_SIM_INSPIRAL_TIDAL_ORDER_ALL:
    case LAL_SIM_INSPIRAL_TIDAL_ORDER_6PN:
      params->wdottidal6pn = lambda1 * XLALSimInspiralTaylorT4wdot_12PNTidalCoeff(params->m1ByM) + lambda2 * XLALSimInspiralTaylorT4wdot_12PNTidalCoeff(params->m2ByM);
      params->Etidal6pn =  lambda1*XLALSimInspiralPNEnergy_12PNTidalCoeff(params->m1ByM) + lambda2*XLALSimInspiralPNEnergy_12PNTidalCoeff(params->m2ByM);
    case LAL_SIM_INSPIRAL_TIDAL_ORDER_5PN:
      params->wdottidal5pn = lambda1 * XLALSimInspiralTaylorT4wdot_10PNTidalCoeff(params->m1ByM) + lambda2 * XLALSimInspiralTaylorT4wdot_10PNTidalCoeff(params->m2ByM);
      params->Etidal5pn = lambda1*XLALSimInspiralPNEnergy_10PNTidalCoeff(params->m1ByM) + lambda2*XLALSimInspiralPNEnergy_10PNTidalCoeff(params->m2ByM);
    case LAL_SIM_INSPIRAL_TIDAL_ORDER_0PN:
      break;
    default:
      XLALPrintError("XLAL Error - %s: Invalid tidal PN order %d\n",
		     __func__, tideO );
      XLAL_ERROR(XLAL_EINVAL);
      break;
  }

  return XLAL_SUCCESS;
} /* End of XLALSimIMRPhenSpinParamsSetup */

static INT4 XLALSpinInspiralDerivatives(UNUSED double t,
                                       const double values[],
                                       double dvalues[],
                                       void *mparams)
{
  REAL8 omega;                // time-derivative of the orbital phase
  REAL8 LNhx, LNhy, LNhz;     // orbital angular momentum unit vector
  REAL8 S1x, S1y, S1z;        // dimension-less spin variable S/M^2
  REAL8 S2x, S2y, S2z;
  REAL8 LNhS1, LNhS2;         // scalar products
  REAL8 domega;               // derivative of omega
  REAL8 dLNhx, dLNhy, dLNhz;  // derivatives of \f$\hat L_N\f$ components
  REAL8 dS1x, dS1y, dS1z;     // derivative of \f$S_i\f$
  REAL8 dS2x, dS2y, dS2z;
  REAL8 energy,energyold;

  /* auxiliary variables*/
  REAL8 S1S2, S1S1, S2S2;     // Scalar products
  REAL8 alphadotcosi;         // alpha is the right ascension of L, i(iota) the angle between L and J
  REAL8 v, v2, v4, v5, v6, v7;
  REAL8 tmpx, tmpy, tmpz, cross1x, cross1y, cross1z, cross2x, cross2y, cross2z, LNhxy;

  LALSimInspiralPhenSpinTaylorT4Coeffs *params = (LALSimInspiralPhenSpinTaylorT4Coeffs *) mparams;

  /* --- computation start here --- */
  omega = values[1];

  LNhx = values[2];
  LNhy = values[3];
  LNhz = values[4];

  S1x = values[5];
  S1y = values[6];
  S1z = values[7];

  S2x = values[8];
  S2y = values[9];
  S2z = values[10];

  energyold = values[11];

  v = cbrt(omega);
  v2 = v * v;
  v4 = omega * v;
  v5 = omega * v2;
  v6 = omega * omega;
  v7 = omega * omega * v;

  // Omega derivative without spin effects up to 3.5 PN
  // this formula does not include the 1.5PN shift mentioned in arXiv:0810.5336, five lines below (3.11)
  domega = params->wdotcoeff[0]
          + v * (params->wdotcoeff[1]
                 + v * (params->wdotcoeff[2]
                        + v * (params->wdotcoeff[3]
                               + v * (params->wdotcoeff[4]
                                      + v * (params->wdotcoeff[5]
                                             + v * (params->wdotcoeff[6] + params->wdotlogcoeff * log(v)
                                                    + v * params->wdotcoeff[7]))))));

  energy = (params->Ecoeff[0] + v2 * (params->Ecoeff[2] +
                                      v2 * (params->Ecoeff[4] +
                                            v2 * params->Ecoeff[6])));

  domega+= omega*v7* ( params->wdottidal5pn + v2 * ( params->wdottidal6pn ) );
  energy+= omega*v7* ( params->Etidal5pn + v2 * params->Etidal6pn);

  // Adding spin effects
  // L dot S1,2
  LNhS1 = (LNhx * S1x + LNhy * S1y + LNhz * S1z);
  LNhS2 = (LNhx * S2x + LNhy * S2y + LNhz * S2z);

  // wdotSO15si = -1/12 (...)
  domega += omega * (params->wdot3S1O * LNhS1 + params->wdot3S2O * LNhS2); // see e.g. Buonanno et al. gr-qc/0211087

  energy += omega * (params->E3S1O * LNhS1 + params->E3S2O * LNhS2);  // see e.g. Blanchet et al. gr-qc/0605140

  // wdotSS2 = -1/48 eta ...
  S1S1 = (S1x * S1x + S1y * S1y + S1z * S1z);
  S2S2 = (S2x * S2x + S2y * S2y + S2z * S2z);
  S1S2 = (S1x * S2x + S1y * S2y + S1z * S2z);
  domega += v4 * ( params->wdot4S1S2 * S1S2 + params->wdot4S1OS2O * LNhS1 * LNhS2);  // see e.g. Buonanno et al. arXiv:0810.5336
  domega += v4 * ( params->wdot4QMS1O * LNhS1*LNhS1 + params->wdot4QMS2O * LNhS2*LNhS2 + params->wdot4QMS1 * S1S1 + params->wdot4QMS2 * S2S2 );
  // see Racine et al. arXiv:0812.4413

  energy += v4 * (params->E4S1S2 * S1S2 + params->E4S1OS2O * LNhS1 * LNhS2);    // see e.g. Buonanno et al. as above
  energy += v4 * (params->E4QMS1 * S1S1 + params->E4QMS2 * S2S2 + params->E4QMS1O * LNhS1 * LNhS1 + params->E4QMS2O * LNhS2 * LNhS2);   // see Racine et al. as above

  // wdotspin25SiLNh = see below
  domega += v5 * (params->wdot5S1O * LNhS1 + params->wdot5S2O * LNhS2);   //see (8.3) of Blanchet et al.
  energy += v5 * (params->E5S1O * LNhS1 + params->E5S2O * LNhS2);    //see (7.9) of Blanchet et al.

  domega += omega*omega * (params->wdot6S1O * LNhS1 + params->wdot6S2O * LNhS2); // see (6.5) of arXiv:1104.5659

  // Setting the right pre-factor
  domega *= params->wdotnewt * v5 * v6;
  energy *= params->Enewt * v2;

  /*Derivative of the angular momentum and spins */

  /* dS1, 1.5PN */
  cross1x = (LNhy * S1z - LNhz * S1y);
  cross1y = (LNhz * S1x - LNhx * S1z);
  cross1z = (LNhx * S1y - LNhy * S1x);

  dS1x = params->S1dot3 * v5 * cross1x;
  dS1y = params->S1dot3 * v5 * cross1y;
  dS1z = params->S1dot3 * v5 * cross1z;

  /* dS1, 2PN */
  tmpx = S1z * S2y - S1y * S2z;
  tmpy = S1x * S2z - S1z * S2x;
  tmpz = S1y * S2x - S1x * S2y;

  // S1S2 contribution see. eq. 2.23 of arXiv:0812.4413
  dS1x += v6 * (params->Sdot4S2*tmpx + params->Sdot4S2O * LNhS2 * cross1x);
  dS1y += v6 * (params->Sdot4S2*tmpy + params->Sdot4S2O * LNhS2 * cross1y);
  dS1z += v6 * (params->Sdot4S2*tmpz + params->Sdot4S2O * LNhS2 * cross1z);
  // S1S1 contribution
  dS1x += v6 * LNhS1 * cross1x * params->S1dot4QMS1O;
  dS1y += v6 * LNhS1 * cross1y * params->S1dot4QMS1O;
  dS1z += v6 * LNhS1 * cross1z * params->S1dot4QMS1O;

  // dS1, 2.5PN, eq. 7.8 of Blanchet et al. gr-qc/0605140
  dS1x += params->S1dot5 * v7 * cross1x;
  dS1y += params->S1dot5 * v7 * cross1y;
  dS1z += params->S1dot5 * v7 * cross1z;

  /* dS2, 1.5PN */
  cross2x = (LNhy * S2z - LNhz * S2y);
  cross2y = (LNhz * S2x - LNhx * S2z);
  cross2z = (LNhx * S2y - LNhy * S2x);

  dS2x = params->S2dot3 * v5 * cross2x;
  dS2y = params->S2dot3 * v5 * cross2y;
  dS2z = params->S2dot3 * v5 * cross2z;

  /* dS2, 2PN */
  dS2x += v6 * (-params->Sdot4S2*tmpx + params->Sdot4S2O * LNhS1 * cross2x);
  dS2y += v6 * (-params->Sdot4S2*tmpy + params->Sdot4S2O * LNhS1 * cross2y);
  dS2z += v6 * (-params->Sdot4S2*tmpz + params->Sdot4S2O * LNhS1 * cross2z);
  // S2S2 contribution
  dS2x += v6 * LNhS2 * cross2x * params->S2dot4QMS2O;
  dS2y += v6 * LNhS2 * cross2y * params->S2dot4QMS2O;
  dS2z += v6 * LNhS2 * cross2z * params->S2dot4QMS2O;

  // dS2, 2.5PN, eq. 7.8 of Blanchet et al. gr-qc/0605140
  dS2x += params->S2dot5 * v7 * cross2x;
  dS2y += params->S2dot5 * v7 * cross2y;
  dS2z += params->S2dot5 * v7 * cross2z;

  dLNhx = -(dS1x + dS2x) * v / params->eta;
  dLNhy = -(dS1y + dS2y) * v / params->eta;
  dLNhz = -(dS1z + dS2z) * v / params->eta;

  /* dphi */
  LNhxy = LNhx * LNhx + LNhy * LNhy;

  if (LNhxy > 0.0)
    alphadotcosi = LNhz * (LNhx * dLNhy - LNhy * dLNhx) / LNhxy;
  else
  {
    //XLALPrintWarning("*** LALSimIMRPSpinInspiralRD WARNING ***: alphadot set to 0, LNh:(%12.4e %12.4e %12.4e)\n",LNhx,LNhy,LNhz);
    alphadotcosi = 0.;
  }

  /* dvalues->data[0] is the phase derivative */
  /* omega is the derivative of the orbital phase omega \neq dvalues->data[0] */
  dvalues[0] = omega - alphadotcosi;
  dvalues[1] = domega;

  dvalues[2] = dLNhx;
  dvalues[3] = dLNhy;
  dvalues[4] = dLNhz;

  dvalues[5] = dS1x;
  dvalues[6] = dS1y;
  dvalues[7] = dS1z;

  dvalues[8] = dS2x;
  dvalues[9] = dS2y;
  dvalues[10] = dS2z;

  dvalues[11] = (energy-energyold)/params->dt*params->M;

  return GSL_SUCCESS;
} /* end of XLALSpinInspiralDerivatives */

static INT4 XLALGenerateWaveDerivative (REAL8Vector *dwave,
                                       REAL8Vector *wave,
                                       REAL8 dt
)
{
  /* XLAL error handling */
  INT4 errcode = XLAL_SUCCESS;

  /* For checking GSL return codes */
  INT4 gslStatus;

  UINT4 j;
  double *x, *y;
  double dy;
  gsl_interp_accel *acc;
  gsl_spline *spline;

  if (wave->length!=dwave->length)
    XLAL_ERROR( XLAL_EFUNC );

  /* Getting interpolation and derivatives of the waveform using gsl spline routine */
  /* Initialize arrays and supporting variables for gsl */

  x = (double *) LALMalloc(wave->length * sizeof(double));
  y = (double *) LALMalloc(wave->length * sizeof(double));

  if ( !x || !y )
  {
    if ( x ) LALFree (x);
    if ( y ) LALFree (y);
    XLAL_ERROR( XLAL_ENOMEM );
  }

  for (j = 0; j < wave->length; ++j)
  {
                x[j] = j;
                y[j] = wave->data[j];
  }

  XLAL_CALLGSL( acc = (gsl_interp_accel*) gsl_interp_accel_alloc() );
  XLAL_CALLGSL( spline = (gsl_spline*) gsl_spline_alloc(gsl_interp_cspline, wave->length) );
  if ( !acc || !spline )
  {
    if ( acc )    gsl_interp_accel_free(acc);
    if ( spline ) gsl_spline_free(spline);
    LALFree( x );
    LALFree( y );
    XLAL_ERROR( XLAL_ENOMEM );
  }

  /* Gall gsl spline interpolation */
  XLAL_CALLGSL( gslStatus = gsl_spline_init(spline, x, y, wave->length) );
  if ( gslStatus != GSL_SUCCESS )
  {
    gsl_spline_free(spline);
    gsl_interp_accel_free(acc);
    LALFree( x );
    LALFree( y );
    XLAL_ERROR( XLAL_EFUNC );
  }

  /* Getting first and second order time derivatives from gsl interpolations */
  for (j = 0; j < wave->length; ++j)
  {
    XLAL_CALLGSL(gslStatus = gsl_spline_eval_deriv_e( spline, j, acc, &dy ) );
    if (gslStatus != GSL_SUCCESS )
    {
      gsl_spline_free(spline);
      gsl_interp_accel_free(acc);
      LALFree( x );
      LALFree( y );
      XLAL_ERROR( XLAL_EFUNC );
    }
    dwave->data[j]  = (REAL8)(dy / dt);
  }

  /* Free gsl variables */
  gsl_spline_free(spline);
  gsl_interp_accel_free(acc);
  LALFree(x);
  LALFree(y);

  return errcode;
} /* End of XLALGenerateWaveDerivative */

static INT4 XLALSimSpinInspiralTest(UNUSED double t, const double values[], double dvalues[], void *mparams) {

  LALSimInspiralPhenSpinTaylorT4Coeffs *params = (LALSimInspiralPhenSpinTaylorT4Coeffs *) mparams;

  REAL8 omega   =   values[1];
  REAL8 energy  =  values[11];
  REAL8 denergy = dvalues[11];

  if ( (energy > 0.0) || (( denergy*params->dt/params->M > - 0.001*energy )&&(energy<0.) ) ) {
    if (energy>0.) XLALPrintWarning("*** Test: LALSimIMRPSpinInspiral WARNING **: Bounding energy >ve!\n");
    else
      XLALPrintWarning("*** Test: LALSimIMRPSpinInspiral WARNING **:  Energy increases dE %12.6e dE*dt %12.6e 1pMEn %12.4e M: %12.4e, eta: %12.4e  om %12.6e \n", denergy, denergy*params->dt/params->M, - 0.001*energy, params->M/LAL_MTSUN_SI, params->eta, omega);
    return LALSIMINSPIRAL_PST4_TEST_ENERGY;
  }
  else if (omega < 0.0) {
    XLALPrintWarning("** LALSimIMRPSpinInspiral WARNING **: omega < 0  M: %12.4e, eta: %12.4e  om %12.6e\n",params->M, params->eta, omega);
    return LALSIMINSPIRAL_PST4_DERIVATIVE_OMEGANONPOS;
  }
  else if (dvalues[1] < 0.0) {
    /* omegadot < 0 */
    return LALSIMINSPIRAL_PST4_TEST_OMEGADOT;
  }
  else if (isnan(omega)) {
    /* omega is nan */
    return LALSIMINSPIRAL_PST4_TEST_OMEGANAN;
  } 
  else if ( params->fEnd > 0. && params->fStart > params->fEnd && omega < params->fEnd) {
    /* freq. below bound in backward integration */
    return LALSIMINSPIRAL_PST4_TEST_FREQBOUND;
  }
  else if ( params->fEnd > params->fStart && omega > params->fEnd) {
    /* freq. above bound in forward integration */
    return LALSIMINSPIRAL_PST4_TEST_FREQBOUND;
  }
  else
    return GSL_SUCCESS;
} /* End of XLALSimSpinInspiralTest */


static INT4 XLALSimIMRPhenSpinTest(UNUSED double t, const double values[], double dvalues[], void *mparams) {

  LALSimInspiralPhenSpinTaylorT4Coeffs *params = (LALSimInspiralPhenSpinTaylorT4Coeffs *) mparams;

  REAL8 omega   =   values[1];
  REAL8 energy  =  values[11];
  REAL8 denergy = dvalues[11];

  REAL8 LNhS1=(values[2]*values[5]+values[3]*values[6]+values[4]*values[7])/params->m1ByM/params->m1ByM;
  REAL8 LNhS2=(values[2]*values[8]+values[3]*values[9]+values[4]*values[10])/params->m2ByM/params->m2ByM;
  REAL8 S1sq =(values[5]*values[5]+values[6]*values[6]+values[7]*values[7])/pow(params->m1ByM,4);
  REAL8 S2sq =(values[8]*values[8]+values[9]*values[9]+values[10]*values[10])/pow(params->m2ByM,4);
  REAL8 S1S2 =(values[5]*values[8]+values[6]*values[9]+values[7]*values[10])/pow(params->m1ByM*params->m2ByM,2);

  REAL8 omegaMatch=OmMatch(LNhS1,LNhS2,S1sq,S1S2,S2sq)+0.0005;

  if ( (energy > 0.0) || (( denergy*params->dt/params->M > - 0.001*energy )&&(energy<0.) ) ) {
    if (energy>0.) XLALPrintWarning("*** Test: LALSimIMRPSpinInspiralRD WARNING **: Bounding energy >ve!\n");
    else
      XLALPrintWarning("*** Test: LALSimIMRPSpinInspiralRD WARNING **:  Energy increases dE %12.6e dE*dt %12.6e 1pMEn %12.4e M: %12.4e, eta: %12.4e  om %12.6e \n", denergy, denergy*params->dt/params->M, - 0.001*energy, params->M/LAL_MTSUN_SI, params->eta, omega);
    return LALSIMINSPIRAL_PST4_TEST_ENERGY;
  }
  else if (omega < 0.0) {
    XLALPrintWarning("** LALSimIMRPSpinInspiralRD WARNING **: omega < 0  M: %12.4e, eta: %12.4e  om %12.6e\n",params->M, params->eta, omega);
    return LALSIMINSPIRAL_PST4_DERIVATIVE_OMEGANONPOS;
  }
  else if (dvalues[1] < 0.0) {
    /* omegadot < 0 */
    return LALSIMINSPIRAL_PST4_TEST_OMEGADOT;
  }
  else if (isnan(omega)) {
    /* omega is nan */
    return LALSIMINSPIRAL_PST4_TEST_OMEGANAN;
  }
  else if ( params->fEnd > 0. && params->fStart > params->fEnd && omega < params->fEnd) {
    /* freq. below bound in backward integration */
    return LALSIMINSPIRAL_PST4_TEST_FREQBOUND;
  }
  else if ( params->fEnd > params->fStart && omega > params->fEnd) {
    /* freq. above bound in forward integration */
    return LALSIMINSPIRAL_PST4_TEST_FREQBOUND;
  }
  else if (omega>omegaMatch) {
    return LALSIMINSPIRAL_PST4_TEST_OMEGAMATCH;
  }
  else
    return GSL_SUCCESS;
} /* End of XLALSimIMRPhenSpinTest */

typedef struct tagLALSimInspiralInclAngle {
  REAL8 cHi;
  REAL8 sHi;
  REAL8 ci;
  REAL8 si;
  REAL8 ci2;
  REAL8 si2;
  REAL8 cHi2;
  REAL8 sHi2;
  REAL8 cHi3;
  REAL8 sHi3;
  REAL8 cHi4;
  REAL8 sHi4;
  REAL8 cHi5;
  REAL8 sHi5;
  REAL8 cHi6;
  REAL8 sHi6;
  REAL8 cHi8;
  REAL8 sHi8;
  REAL8 cDi;
  REAL8 sDi;
} LALSimInspiralInclAngle;

static INT4 XLALSimSpinInspiralFillL2Modes(COMPLEX16Vector *hL2,
                                   REAL8 v,
                                   REAL8 eta,
                                   REAL8 dm,
                                   REAL8 Psi,
                                   REAL8 alpha,
                                   LALSimInspiralInclAngle *an
                                   )
{
  const INT4 os=2;
  REAL8 amp20 = sqrt(1.5);
  REAL8 v2    = v*v;
  REAL8 damp  = 1.;

  hL2->data[2+os] = ( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
                      ( cos(2.*(Psi+alpha)) * an->cHi4 + cos(2.*(Psi-alpha)) * an->sHi4 ) +
                      v * dm/3.*an->si * ( cos(Psi-2.*alpha) * an->sHi2 + cos(Psi + 2.*alpha) * an->cHi2 ) );

  hL2->data[2+os]+=I*( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
                       (-sin(2.*(Psi+alpha)) * an->cHi4 + sin(2.*(Psi-alpha)) * an->sHi4 ) +
                       v * dm/3.*an->si * ( sin(Psi-2.*alpha) * an->sHi2 - sin(Psi + 2. * alpha) * an->cHi2 ) );

  hL2->data[-2+os] = ( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
                       ( cos(2. * (Psi + alpha)) * an->cHi4 + cos(2. * (Psi - alpha)) * an->sHi4 ) -
                       v * dm / 3. * an->si * ( cos(Psi - 2. * alpha) * an->sHi2 + cos(Psi + 2. * alpha) * an->cHi2 ) );

  hL2->data[-2+os]+=I*( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
                        ( sin(2.*(Psi + alpha))*an->cHi4 - sin(2.*(Psi-alpha)) * an->sHi4 ) +
                         v*dm/3.*an->si * ( sin(Psi-2.*alpha) * an->sHi2 - sin(Psi+2.*alpha) * an->cHi2 ) );

  hL2->data[1+os] = an->si * ( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
                                ( -cos(2. * Psi - alpha) * an->sHi2 + cos(2. * Psi + alpha) * an->cHi2 ) +
                                v * dm / 3. * ( -cos(Psi + alpha) * (an->ci + an->cDi)/2. - cos(Psi - alpha) * an->sHi2 * (1. + 2. * an->ci) ) );

  hL2->data[1+os]+= an->si *I*( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
                                ( -sin(2.*Psi-alpha ) * an->sHi2 - sin(2.*Psi + alpha) * an->cHi2 ) +
                                v * dm / 3. * (sin(Psi + alpha) * (an->ci + an->cDi)/2. - sin(Psi - alpha) * an->sHi2 * (1.+2.*an->ci) ) );

  hL2->data[-1+os] = an->si * ( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
                                ( cos(2.*Psi-alpha) * an->sHi2 - cos(2.*Psi+alpha)*an->cHi2) +
                                v * dm / 3. * ( -cos(Psi + alpha) * (an->ci + an->cDi)/2. - cos(Psi - alpha) * an->sHi2 * (1. + 2. * an->ci) ) );

  hL2->data[-1+os]+= an->si *I*( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
                                 ( -sin(2. * Psi - alpha) * an->sHi2 - sin(2. * Psi + alpha) * an->cHi2 ) -
                                 v * dm / 3. * ( sin(Psi + alpha) * (an->ci + an->cDi)/2. - sin(Psi - alpha) * an->sHi2 * (1. + 2. * an->ci) ) );

  hL2->data[os] = amp20 * ( an->si2/( 1. + damp *v2/42.*(107.-55.*eta) )*cos(2.*Psi) + I*v*dm/3.*an->sDi*sin(Psi) );

  return XLAL_SUCCESS;
} /* End of XLALSimSpinInspiralFillL2Modes*/

static INT4 XLALSimSpinInspiralFillL3Modes(COMPLEX16Vector *hL3,
                                   REAL8 v,
                                   REAL8 eta,
                                   REAL8 dm,
                                   REAL8 Psi,
                                   REAL8 alpha,
                                   LALSimInspiralInclAngle *an)
{
  const INT4 os=3;
  REAL8 amp32 = sqrt(1.5);
  REAL8 amp31 = sqrt(0.15);
  REAL8 amp30 = 1. / sqrt(5)/2.;
  REAL8 v2    = v*v;

  hL3->data[3+os] = (v * dm * (-9.*cos(3.*(Psi-alpha))*an->sHi6 - cos(Psi-3.*alpha)*an->sHi4*an->cHi2 + cos(Psi+3.*alpha)*an->sHi2*an->cHi4 + 9.*cos(3.*(Psi+alpha))*an->cHi6) +
                     v2 * 4. * an->si *(1.-3.*eta)* ( -cos(2.*Psi-3.*alpha)*an->sHi4 + cos(2.*Psi+3.*alpha)*an->cHi4) );

  hL3->data[3+os]+= I*(v * dm * (-9.*sin(3.*(Psi-alpha))*an->sHi6 - sin(Psi-3.*alpha)*an->sHi4*an->cHi2 - sin(Psi+3.*alpha)*an->sHi2*an->cHi4 - 9.*sin(3.*(Psi+alpha))* an->cHi6) +
                       v2 * 4. * an->si *(1.-3.*eta)* ( -sin(2.*Psi-3.*alpha)*an->sHi4 -sin(2.*Psi+3.*alpha)*an->cHi4) );

  hL3->data[-3+os] = (-v * dm * (-9.*cos(3.*(Psi-alpha))*an->sHi6 - cos(Psi-3.*alpha)*an->sHi4*an->cHi2 + cos(Psi+3.*alpha)*an->sHi2*an->cHi4 + 9.*cos(3.*(Psi+alpha))*an->cHi6) +
                      v2 * 4. * an->si *(1.-3.*eta)*( -cos(2.*Psi-3.*alpha)*an->sHi4 + cos(2.*Psi+3.*alpha)*an->cHi4) );

  hL3->data[-3+os]+=I*(v * dm *(-9.*sin(3.*(Psi-alpha))*an->sHi6 - sin(Psi-3.*alpha)*an->sHi4*an->cHi2 - sin(Psi+3.*alpha)*an->sHi2*an->cHi4 - 9.*sin(3.*(Psi+alpha))* an->cHi6) -
                       v2 * 4. * an->si * (1.-3.*eta)*( -sin(2.*Psi-3.*alpha)*an->sHi4 - sin(2.*Psi+3.*alpha)*an->cHi4 ) );

  hL3->data[2+os] = amp32 * ( v * dm/3. * (27.*cos(3.*Psi-2.*alpha)*an->si*an->sHi4 + 27.*cos(3.*Psi+2.*alpha)*an->si*an->cHi4 + cos(Psi+2.*alpha)*an->cHi3*(5.*an->sHi-3.*an->si*an->cHi-3.*an->ci*an->sHi) /2. + cos(Psi-2.*alpha)*an->sHi3*(5.*an->cHi+3.*an->ci*an->cHi-3.*an->si*an->sHi) /2. ) +
                              v2*(1./3.-eta) * (-8.*an->cHi4*(3.*an->ci-2.)*cos(2.*(Psi+alpha)) + 8.*an->sHi4*(3.*an->ci+2.)*cos(2.*(Psi-alpha)) ) );

  hL3->data[2+os]+= amp32*I*( v * dm/3. * ( 27.*sin(3.*Psi-2.*alpha)*an->si*an->sHi4 - 27.*cos(3.*Psi+2.*alpha)*an->si*an->cHi4 - sin(Psi+2.*alpha)*an->cHi3*(5.*an->sHi-3.*an->si*an->cHi-3.*an->ci*an->sHi) /2. + sin(Psi-2.*alpha)*an->sHi3*(5.*an->cHi+3.*an->ci*an->cHi-3.*an->si*an->sHi)/2. ) +
                              v2*(1./3.-eta) * ( 8.*an->cHi4*(3.*an->ci-2.)*sin(2.*(Psi+alpha)) + 8.*an->sHi4*(3.*an->ci+2.)*sin(2.*(Psi-alpha)) ) );

  hL3->data[-2+os] = amp32 * ( v * dm/3. * (27.*cos(3.*Psi-2.*alpha)*an->si*an->sHi4 + 27.*cos(3.*Psi+2.*alpha)*an->si*an->cHi4 + cos(Psi+2.*alpha)*an->cHi3*(5.*an->sHi-3.*an->si*an->cHi-3.*an->ci*an->sHi) /2. + cos(Psi-2.*alpha)*an->sHi3*(5.*an->cHi+3.*an->ci*an->cHi-3.*an->si*an->sHi) /2. ) -
                               v2*(1./3.-eta) * ( 8.*an->cHi4*(3.*an->ci-2.)*cos(2.*(Psi+alpha)) - 8.*an->sHi4*(3.*an->ci+2.)*cos(2.*(Psi-alpha)) ) );

  hL3->data[-2+os]+= amp32*I*(-v * dm/3. * (27.*sin(3.*Psi-2.*alpha)*an->si*an->sHi4 - 27.*cos(3.*Psi+2.*alpha)*an->si*an->cHi4 - sin(Psi+2.*alpha)*an->cHi3*(5.*an->sHi-3.*an->si*an->cHi-3.*an->ci*an->sHi) /2.+ sin(Psi-2.*alpha)*an->sHi3*(5.*an->cHi+3.*an->ci*an->cHi-3.*an->si*an->sHi) /2.) +
                             v2*(1./3.-eta) * (8.*an->cHi4*(3.*an->ci-2.)*sin(2.*(Psi+alpha)) + 8.*an->sHi4*(3.*an->ci+2.)*sin(2.*(Psi-alpha)) ) );

  hL3->data[1+os] = amp31 * ( v * dm/6. * ( -135.*cos(3.*Psi-alpha)*an->sHi*an->sHi2 + 135.*cos(3.*Psi+alpha)*an->sHi*an->cHi2 + cos(Psi+alpha)*an->cHi2*(15.*an->cDi-20.*an->ci+13.)/2. - cos(Psi-alpha)*an->sHi2*(15.*an->cDi+20.*an->ci+13.)/2.)
                            + v2*(1./3.-eta) * ( 20.*an->cHi3*cos(2.*Psi+alpha)*(3.*(an->sHi*an->ci+an->cHi*an->si)-5.*an->sHi) + 20.*an->sHi3*cos(2.*Psi-alpha)*(3.*(an->cHi2*an->ci-an->sHi*an->si)+5.*an->cHi) ) );

  hL3->data[1+os]+= amp31*I*(-v * dm/6. * ( -135.*cos(3.*Psi-alpha)*an->si2*an->sHi2 + 135.*cos(3.*Psi+alpha)*an->si2*an->cHi2 + cos(Psi+alpha)*an->cHi2*(15.*an->cDi-20.*an->ci+13.)/2. - cos(Psi-alpha)*an->sHi2*(15.*an->cDi+20.*an->ci+13.)/2. )
                           - v2*(1./3.-eta) * ( 20.*an->cHi3*cos(2.*Psi+alpha)*(3.*(an->sHi*an->ci+an->cHi*an->si)-5.*an->sHi) + 20.*an->sHi3*cos(2.*Psi-alpha)*(3.*(an->cHi*an->ci-an->sHi*an->si)+5.*an->cHi) ) );

  hL3->data[-1+os] = amp31 * (-v * dm/6. * ( -135.*cos(3.*Psi-alpha)*an->si2*an->sHi2 + 135.*cos(3.*Psi+alpha)*an->si2*an->cHi2 + cos(Psi+alpha)*an->cHi2*(15.*an->cDi-20.*an->ci+13.)/2.- cos(Psi-alpha) * an->sHi2*(15.*an->cDi+20.*an->ci+13.)/2. ) -
                               v2 * (1./3.-eta)* ( 20.*an->cHi3*cos(2.*Psi+alpha)*(3.*(an->sHi*an->ci+an->cHi*an->si)-5.*an->sHi) + 20.*an->sHi3*cos(2.*Psi-alpha)*(3.*(an->cHi*an->ci-an->sHi*an->si)+5.*an->cHi) ) );

  hL3->data[-1+os]+= amp31*I*(v * dm/6. * ( -135.*sin(3.*Psi-alpha)*an->si2*an->sHi2 - 135.*sin(3.*Psi+alpha)*an->si2*an->cHi2 - sin(Psi+alpha)*an->cHi2*(15.*an->cDi-20.*an->ci+13.)/2. - sin(Psi-alpha)*an->sHi2*(15.*an->cDi+20.*an->ci+13.)/2.)
                              -v2 * (1./3.-eta)* ( 20.*an->cHi3*sin(2.*Psi+alpha)*(3.*(an->sHi*an->ci+an->ci2*an->si)-5.*an->si2) - 20.*an->sHi3*sin(2.*Psi-alpha)*(3.*(an->ci2*an->ci-an->si2*an->si)+5.*an->ci2) ) );

  hL3->data[os] = amp30 * I * ( v * dm * ( cos(Psi)*an->si*(cos(2.*Psi)*(45.*an->si2)-(25.*an->cDi-21.) ) ) +
                                v2*(1.-3.*eta) * (80.*an->si2*an->cHi*sin(2.*Psi) ) );

  return XLAL_SUCCESS;

} /*End of XLALSimSpinInspiralFillL3Modes*/

static INT4 XLALSimSpinInspiralFillL4Modes(COMPLEX16Vector *hL4,
                                   UNUSED REAL8 v,
                                   REAL8 eta,
                                   UNUSED REAL8 dm,
                                   REAL8 Psi,
                                   REAL8 alpha,
                                   LALSimInspiralInclAngle *an
                                   )
{
  const INT4 os=4;
  REAL8 amp43 = - sqrt(2.);
  REAL8 amp42 = sqrt(7.)/2.;
  REAL8 amp41 = sqrt(3.5)/4.;
  REAL8 amp40 = sqrt(17.5)/16.;

  hL4->data[4+os] = (1. - 3.*eta) * ( 4.*an->sHi8*cos(4.*(Psi-alpha)) + cos(2.*Psi-4.*alpha)*an->sHi6*an->cHi2 + an->sHi2*an->cHi6*cos(2.*Psi+4.*alpha) + 4.*an->cHi8*cos(4.*(Psi+alpha)) );

  hL4->data[4+os]+= (1. - 3.*eta)*I*( 4.*an->sHi8*sin(4.*(Psi-alpha)) + sin(2.*Psi-4.*alpha)*an->sHi6*an->cHi2 - an->sHi2*an->cHi6*sin(2.*Psi+4.*alpha) - 4.*an->cHi8*sin(4.*(Psi+alpha)) );

  hL4->data[-4+os] = (1. - 3.*eta) * (4.*an->sHi8*cos(4.*(Psi-alpha)) + cos(2.*Psi-4.*alpha)*an->sHi6*an->cHi2 + an->sHi2*an->cHi6*cos(2.*Psi+4.*alpha) + 4.*an->cHi8*cos(4.*(Psi+alpha) ) );

  hL4->data[-4+os]+=-(1. - 3.*eta) *I*(4.*an->sHi8*sin(4.*(Psi-alpha)) + sin(2*Psi-4.*alpha)*an->sHi6*an->cHi2 - an->sHi2*an->cHi6*sin(2.*Psi+4.*alpha) - 4.*an->cHi8*sin(4.*(Psi+alpha)) );

  hL4->data[3+os] = amp43 * (1. - 3.*eta) * an->si * ( 4.*an->sHi6*cos(4.*Psi-3.*alpha) - 4.*an->cHi6*cos(4.*Psi+3.*alpha) - an->sHi4*(an->ci+0.5)/2.*cos(2.*Psi-3.*alpha) + an->cHi4*(an->ci-0.5)*cos(2.*Psi+3.*alpha) ); /****/

  hL4->data[3+os]+= amp43*I*(1. - 3.*eta) * an->si * ( 4.*an->sHi6*sin(4.*Psi-3.*alpha) + 4.*an->cHi6*sin(4.*Psi+3.*alpha) - an->sHi4*(an->ci+0.5)/2.*sin(2.*Psi-3.*alpha) + an->cHi4*(an->ci-0.5)*sin(2.*Psi+3.*alpha) ); /****/

  hL4->data[-3+os] = -amp43 * (1. - 3.*eta) * an->si * ( 4.*an->sHi6*cos(4.*Psi-3.*alpha) - 4.*an->cHi6*cos(4.*Psi+3.*alpha) - an->sHi4*(an->ci+0.5)/2.*cos(2.*Psi-3.*alpha) + an->cHi4*(an->ci-0.5)*cos(2.*Psi+3.*alpha) ); /****/

  hL4->data[-3+os]+= amp43*I*(1. - 3.*eta) * an->si * ( 4.*an->sHi6*sin(4.*Psi-3.*alpha) + 4.*an->cHi6*sin(4.*Psi+3.*alpha) - an->sHi4*(an->ci+0.5)/2.*sin(2.*Psi-3.*alpha) + an->cHi4*(an->ci-0.5)*sin(2.*Psi+3.*alpha) ); /****/

  hL4->data[2+os] = amp42 * (1. - 3.*eta) * ( 16.*an->sHi6*an->cHi2*cos(4.*Psi-2.*alpha) + 16.*an->cHi6*an->sHi2*cos(4.*Psi+2.*alpha) - an->cHi4*cos(2.*(Psi+alpha))*(an->cDi-2.*an->ci+9./7.)/2. - an->sHi4*cos(2.*(Psi-alpha))*(an->cDi+2.*an->ci+9./7.)/2. );

  hL4->data[2+os]+= amp42 *I*(1. - 3.*eta) * ( 16.*an->sHi6*an->cHi2 * sin(4.*Psi-2.*alpha) - 16.*an->cHi6*an->sHi2*sin(4.*Psi+2.*alpha) + an->cHi4*sin(2.*(Psi+alpha))*(an->cDi-2.*an->ci+9./7.)/2. - an->sHi4*sin(2.*(Psi-alpha))*(an->cDi+2.*an->ci+9./7.)/2. );

  hL4->data[-2+os] = amp42 * (1. - 3.*eta) * ( 16.*an->sHi6*an->cHi2*cos(4.*Psi-2.*alpha) + 16.*an->cHi6*an->sHi2*cos(4.*Psi+2.*alpha) - an->cHi4*cos(2.*(Psi+alpha))*(an->cDi-2.*an->ci+9./7.)/2. - an->sHi4*cos(2.*(Psi-alpha))*(an->cDi+2.*an->ci+9./7.)/2. );

  hL4->data[-2+os]+=-amp42 *I*(1. - 3.*eta) * ( 16.*an->sHi6*an->cHi2*sin(4.*Psi-2.*alpha) - 16.*an->cHi6*an->sHi2*sin(4.*Psi+2.*alpha) + an->cHi4*sin(2.*(Psi+alpha))*(an->cDi-2.*an->ci+9./7.)/2. - an->sHi4*sin(2.*(Psi-alpha))*(an->cDi+2.*an->ci+9./7.)/2. );

  hL4->data[1+os] = amp41 * (1. - 3.*eta) * ( -64.*an->sHi5*an->cHi3*cos(4.*Psi-alpha) + 64.*an->sHi3*an->cHi5*cos(4.*Psi+alpha) - an->sHi3*cos(2.*Psi-alpha)*((an->cDi*an->cHi-an->sDi*an->sHi) + 2.*(an->cHi*an->ci-an->sHi*an->si) + 19./7.*an->cHi) + an->cHi3*cos(2.*Psi+alpha)*((an->cDi*an->sHi+an->sDi*an->cHi) - 2.*(an->si*an->cHi+an->ci*an->si2) +19./7.*an->cHi) );

  hL4->data[1+os]+= amp41*I*(1. - 3.*eta) * ( -64.*an->sHi5*an->cHi3 * sin(4.*Psi-alpha) - 64.*an->sHi3*an->cHi5 * sin(4.*Psi+alpha) - an->sHi3*sin(2.*Psi-alpha)*((an->cDi*an->cHi-an->sDi*an->sHi) + 2.*(an->cHi*an->ci-an->sHi*an->si) + 19./7.*an->cHi) - an->cHi3*sin(2.*Psi+alpha)*((an->cDi*an->sHi+an->sDi*an->cHi) - 2.*(an->si*an->cHi+an->ci*an->sHi) + 19./7.*an->cHi) );

  hL4->data[-1+os] = -amp41 * (1. - 3.*eta) * ( -64*an->sHi5*an->cHi3 * cos(4.*Psi-alpha) + 64.*an->sHi3*an->cHi5*cos(4.*Psi+alpha) - an->sHi3*cos(2.*Psi-alpha)*((an->cDi*an->cHi-an->sDi*an->sHi) + 2.*(an->cHi*an->ci-an->sHi*an->si) + 19./7.*an->ci2) + an->cHi3*cos(2.*Psi+alpha)*((an->cDi*an->sHi+an->sDi*an->cHi) - 2.*(an->si*an->cHi+an->ci*an->sHi) + 19./7.*an->ci2) );

  hL4->data[-1+os]+= amp41 *I*(1. - 3.*eta) *I*( -64.*an->sHi5*an->cHi3 * sin(4.*Psi-alpha) - 64.*an->sHi3*an->cHi5 * sin(4.*Psi+alpha) - an->sHi3*sin(2.*Psi-alpha)*((an->cDi*an->cHi-an->sDi*an->sHi) + 2.*(an->cHi*an->ci-an->sHi*an->si) + 19./7.*an->ci2) - an->cHi3*sin(2.*Psi+alpha)*((an->cDi*an->sHi+an->sDi*an->cHi) - 2.*(an->si*an->cHi+an->ci*an->sHi) + 19./7.*an->cHi) );

  hL4->data[os] = amp40 * (1.-3.*eta) * an->si2 * (8.*an->si2*cos(4.*Psi) + cos(2.*Psi)*(an->cDi+5./7.) );

  return XLAL_SUCCESS;
} /* End of XLALSimSpinInspiralFillL4Modes*/

static INT4 XLALSimInspiralSpinTaylorT4Engine(REAL8TimeSeries **omega,      /**< post-Newtonian parameter [returned]*/
                                             REAL8TimeSeries **Phi,        /**< orbital phase            [returned]*/
                                             REAL8TimeSeries **LNhatx,     /**< LNhat vector x component [returned]*/
                                             REAL8TimeSeries **LNhaty,     /**< "    "    "  y component [returned]*/
                                             REAL8TimeSeries **LNhatz,     /**< "    "    "  z component [returned]*/
                                             REAL8TimeSeries **S1x,        /**< Spin1 vector x component [returned]*/
                                             REAL8TimeSeries **S1y,        /**< "    "    "  y component [returned]*/
                                             REAL8TimeSeries **S1z,        /**< "    "    "  z component [returned]*/
                                             REAL8TimeSeries **S2x,        /**< Spin2 vector x component [returned]*/
                                             REAL8TimeSeries **S2y,        /**< "    "    "  y component [returned]*/
                                             REAL8TimeSeries **S2z,        /**< "    "    "  z component [returned]*/
                                             REAL8TimeSeries **Energy,     /**< Energy                   [returned]*/
                                             const REAL8 yinit[],          /**< UNDOCUMENTED */
                                             const INT4  lengthH,          /**< UNDOCUMENTED */
                                             const Approximant approx,     /**< Allow to choose w/o ringdown */
                                             LALSimInspiralPhenSpinTaylorT4Coeffs *params /**< UNDOCUMENTED */
                                             )
{
  UINT4 idx;
  INT4 jdx;
  UINT4 intLen;
  INT4 intReturn;

  REAL8 S1x0,S1y0,S1z0,S2x0,S2y0,S2z0;  /** Used to store initial spin values */
  REAL8Array *yout;                     /** Used to store integration output */

  LALAdaptiveRungeKutta4Integrator *integrator;

  /* allocate the integrator */
  if (approx == PhenSpinTaylor)
    integrator = XLALAdaptiveRungeKutta4Init(LAL_NUM_PST4_VARIABLES,XLALSpinInspiralDerivatives,XLALSimSpinInspiralTest,LAL_PST4_ABSOLUTE_TOLERANCE,LAL_PST4_RELATIVE_TOLERANCE);
  else
    integrator = XLALAdaptiveRungeKutta4Init(LAL_NUM_PST4_VARIABLES,XLALSpinInspiralDerivatives,XLALSimIMRPhenSpinTest,LAL_PST4_ABSOLUTE_TOLERANCE,LAL_PST4_RELATIVE_TOLERANCE);

  if (!integrator) {
    XLALPrintError("XLAL Error - %s: Cannot allocate integrator\n", __func__);
    XLAL_ERROR(XLAL_EFUNC);
  }

  /* stop the integration only when the test is true */
  integrator->stopontestonly = 1;

  REAL8 *yin = (REAL8 *) LALMalloc(sizeof(REAL8) * LAL_NUM_PST4_VARIABLES);
  for (idx=0; idx<LAL_NUM_PST4_VARIABLES; idx++) yin[idx]=yinit[idx];
  S1x0=yinit[5];
  S1y0=yinit[6];
  S1z0=yinit[7];
  S2x0=yinit[8];
  S2y0=yinit[9];
  S2z0=yinit[10];

  //REAL8 dtInt=1./OmMatch(0,0,0,0,0)/50.*fabs(params->dt)/params->dt;
  REAL8 length=((REAL8)lengthH)*fabs(params->dt)/params->M;
  intLen    = XLALAdaptiveRungeKutta4Hermite(integrator,(void *)params,yin,0.0,length,params->dt/params->M,&yout);

  intReturn = integrator->returncode;
  XLALAdaptiveRungeKutta4Free(integrator);

  if (intReturn == XLAL_FAILURE) {
    XLALPrintError("** LALSimIMRPSpinInspiralRD Error **: Adaptive Integrator\n");
    XLALPrintError("             m:  %12.4e  %12.4e  Mom  %12.4e\n",params->m1ByM*params->M,params->m2ByM*params->M,params->fStart);
    XLALPrintError("             S1: %12.4e  %12.4e  %12.4e\n",S1x0,S1y0,S1z0);
    XLALPrintError("             S2: %12.4e  %12.4e  %12.4e\n",S2x0,S2y0,S2z0);
    XLAL_ERROR(XLAL_EFUNC);
  }
  /* End integration*/

  /* Start of the integration checks*/
  if (intLen<minIntLen) {
    XLALPrintError("** LALSimIMRPSpinInspiralRD ERROR **: integration too short! intReturnCode %d, integration length %d, at least %d required\n",intReturn,intLen,minIntLen);
    if (XLALClearErrno() == XLAL_ENOMEM) {
      XLAL_ERROR(  XLAL_ENOMEM);
    } else {
      XLAL_ERROR( XLAL_EFAILED);
    }
  }

  const LIGOTimeGPS tStart=LIGOTIMEGPSZERO;
  *omega  = XLALCreateREAL8TimeSeries( "OMEGA", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *Phi    = XLALCreateREAL8TimeSeries( "ORBITAL_PHASE", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *LNhatx = XLALCreateREAL8TimeSeries( "LNHAT_X_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *LNhaty = XLALCreateREAL8TimeSeries( "LNHAT_Y_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *LNhatz = XLALCreateREAL8TimeSeries( "LNHAT_Z_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *S1x    = XLALCreateREAL8TimeSeries( "SPIN1_X_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *S1y    = XLALCreateREAL8TimeSeries( "SPIN1_Y_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *S1z    = XLALCreateREAL8TimeSeries( "SPIN1_Z_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *S2x    = XLALCreateREAL8TimeSeries( "SPIN2_X_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *S2y    = XLALCreateREAL8TimeSeries( "SPIN2_Y_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *S2z    = XLALCreateREAL8TimeSeries( "SPIN2_Z_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  *Energy = XLALCreateREAL8TimeSeries( "LNHAT_Z_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
  if ( !omega || !Phi || !S1x || !S1y || !S1z || !S2x || !S2y || !S2z || !LNhatx || !LNhaty || !LNhatz || !Energy ) {
    XLALDestroyREAL8Array(yout);
    XLAL_ERROR(XLAL_EFUNC);
  }

  /* Copy dynamical variables from yout array to output time series.
   * Note the first 'len' members of yout are the time steps.
   */
  INT4 sign=params->dt > 0. ? 1 : -1;
  jdx= (intLen-1)*(-sign+1)/2;

  for (idx=0;idx<intLen;idx++) {
    (*Phi)->data->data[idx]    = yout->data[intLen+jdx];
    (*omega)->data->data[idx]  = yout->data[2*intLen+jdx];
    (*LNhatx)->data->data[idx] = yout->data[3*intLen+jdx];
    (*LNhaty)->data->data[idx] = yout->data[4*intLen+jdx];
    (*LNhatz)->data->data[idx] = yout->data[5*intLen+jdx];
    (*S1x)->data->data[idx]    = yout->data[6*intLen+jdx];
    (*S1y)->data->data[idx]    = yout->data[7*intLen+jdx];
    (*S1z)->data->data[idx]    = yout->data[8*intLen+jdx];
    (*S2x)->data->data[idx]    = yout->data[9*intLen+jdx];
    (*S2y)->data->data[idx]    = yout->data[10*intLen+jdx];
    (*S2z)->data->data[idx]    = yout->data[11*intLen+jdx];
    (*Energy)->data->data[idx] = yout->data[12*intLen+jdx];
    jdx+=sign;
  }

  XLALDestroyREAL8Array(yout);