LALInspiralEccentricity.c

Go to the documentation of this file.

/*
*  Copyright (C) 2007
*
*  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., 51 Franklin Street, Fifth Floor, Boston,
*  MA  02110-1301  USA
*/
 
/**
 * \author Devanka Pathak and Thomas Cokelaer
 * \file
 *
 * \brief The code \c LALInspiralEccentricity generates a time-domain inspiral waveform corresponding to the
 * \c approximant \c Eccentricity as outlined PRD 60 for the Newtonian case.
 *
 * ### Prototypes ###
 *
 * <tt>LALInspiralEccentricity()</tt>
 * <ul>
 * <li> \c signalvec: Output containing the inspiral waveform.</li>
 * <li> \c params: Input containing binary chirp parameters and eccentricity.</li>
 * </ul>
 *
 * <tt>LALInspiralEccentricityTemplates()</tt>
 * <ul>
 * <li> \c signalvec1: Output containing the 0-phase inspiral waveform.</li>
 * <li> \c signalvec2: Output containing the \f$\pi/2\f$-phase inspiral waveform.</li>
 * <li> \c params: Input containing binary chirp parameters.</li>
 * </ul>
 *
 * ### Description ###
 *
 * \c LALInspiralEccentricity is called if the user has specified the
 * \c enum \c approximant to be
 * either \c TaylorT1 or \c PadeT1.
 * \c LALInspiralEccentricityTemplates is exactly the same as <tt>LALInspiralEccentricity,</tt> except that
 * it generates two templates one for which the starting phase is
 * <tt>params.startPhase</tt> and the other for which the phase is
 * <tt>params.startPhase + \f$\pi/2\f$</tt>.
 *
 * ### Algorithm ###
 *
 * This code uses a fourth-order Runge-Kutta algorithm to solve the ODEs
 * in \eqref{eq_ode2}.
 *
 * ### Uses ###
 *
 * \code
 * LALInspiralSetup()
 * LALInspiralChooseModel()
 * LALInspiralVelocity()
 * LALInspiralPhasing1()
 * LALInspiralDerivatives()
 * LALRungeKutta4()
 * \endcode
 *
 * ### Notes ###
 *
 */
 
/*
   Interface routine needed to generate time-domain T- or a P-approximant
   waveforms by solving the ODEs using a 4th order Runge-Kutta; April 5, 00.
*/
#include <math.h>
#include <lal/LALInspiral.h>
#include <lal/LALStdlib.h>
#include <lal/Units.h>
#include <lal/SeqFactories.h>
 
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif
 
/* structure to provide M and eta. */
typedef struct
tagecc_CBC_ODE_Struct
{
  double totalMass; /* M=m1+m2 */
  double eta;       /* eta=m1 m2/M^2 */
}
ecc_CBC_ODE_Input;
 
 
void
LALInspiralEccentricityDerivatives (
   REAL8Vector *values,
   REAL8Vector *dvalues,
   void        *params
   );
 
static void
LALInspiralEccentricityEngine(
   LALStatus        *status,
   REAL4Vector      *signalvec1,
   REAL4Vector      *signalvec2,
   REAL4Vector      *a,
   REAL4Vector      *ff,
   REAL8Vector      *phi,
   INT4             *countback,
   InspiralTemplate *params
   );
 
void
LALInspiralEccentricity(
   LALStatus        *status,
   REAL4Vector      *signalvec,
   InspiralTemplate *params
   )
 {
 
   INT4 count;
 
   INITSTATUS(status);
   ATTATCHSTATUSPTR(status);
 
   ASSERT(signalvec, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signalvec->data, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
 
 
   /* Initially the waveform is empty*/
   memset(signalvec->data, 0, signalvec->length*sizeof(REAL4));
 
   /*Call the engine function*/
   LALInspiralEccentricityEngine(status->statusPtr, signalvec, NULL, NULL, NULL, NULL, &count, params);
   CHECKSTATUSPTR(status);
 
   DETATCHSTATUSPTR(status);
   RETURN (status);
 
}
 
 
 
/*
   Interface routine needed to generate time-domain T- or a P-approximant
   waveforms by solving the ODEs using a 4th order Runge-Kutta; April 5, 00.
*/
 
void
LALInspiralEccentricityTemplates(
   LALStatus        *status,
   REAL4Vector      *signalvec1,
   REAL4Vector      *signalvec2,
   InspiralTemplate *params
   )
 {
 
   INT4 count;
 
   INITSTATUS(status);
   ATTATCHSTATUSPTR(status);
 
   ASSERT(signalvec1, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signalvec2, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signalvec1->data, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signalvec2->data, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
 
   /* Initially the waveforms are empty */
   memset(signalvec1->data, 0, signalvec1->length * sizeof(REAL4));
   memset(signalvec2->data, 0, signalvec2->length * sizeof(REAL4));
 
   /* Call the engine function */
   LALInspiralEccentricityEngine(status->statusPtr, signalvec1, signalvec2, NULL, NULL, NULL, &count, params);
   CHECKSTATUSPTR(status);
 
   DETATCHSTATUSPTR(status);
   RETURN (status);
 
}
 
/*
 *  Engine function for use by other LALInspiralEccentricity* functions
 *  Craig Robinson April 2005
 */
 
void
LALInspiralEccentricityEngine(
                LALStatus        *status,
                REAL4Vector      *signalvec1,
                REAL4Vector      *signalvec2,
                REAL4Vector      *a,
                REAL4Vector      UNUSED *ff,
                REAL8Vector      UNUSED *phi,
                INT4             *countback,
                InspiralTemplate *params)
{
   INT4 number_of_diff_equations = 3;
   INT4 count = 0;
   ecc_CBC_ODE_Input in3;
   REAL8 phase;
   REAL8 orbital_element_p,orbital_element_e_squared;
   REAL8 orbital_element_e;
   REAL8 twoPhim2Beta = 0;
   REAL8 threePhim2Beta = 0;
   REAL8 phim2Beta = 0;
   REAL8 twoBeta;
   REAL8 rbyM=1e6, rbyMFlso=6.;
   REAL8 sin2Beta,cos2Beta,iota,onepCosSqI, SinSqI, cosI, e0, f_min, beta, p0;
 
 
 
   REAL8 amp, m, dt, t,  h1, h2, f,  fHigh, piM, fu;
   REAL8Vector dummy, values, dvalues, valuesNew, yt, dym, dyt;
   INT4 done=0;
   rk4In in4;
   rk4GSLIntegrator *integrator;
   void *funcParams;
   expnCoeffs ak;
   expnFunc func;
 
#if 0
   REAL8 mTot = 0;
   REAL8 unitHz = 0;
   REAL8 f2a = 0;
   REAL8 mu = 0;
   REAL8 etab = 0;
   REAL8 fFac = 0; /* SI normalization for f and t */
   REAL8 f2aFac = 0;/* factor multiplying f in amplitude function */
   REAL8 apFac = 0, acFac = 0;/* extra factor in plus and cross amplitudes */
#endif
 
   INITSTATUS(status);
   ATTATCHSTATUSPTR(status);
 
   ASSERT (params,  status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT (params->nStartPad >= 0, status, LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);
   ASSERT (params->nEndPad >= 0, status, LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);
   ASSERT (params->fLower > 0, status, LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);
   ASSERT (params->tSampling > 0, status, LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);
 
   LALInspiralSetup (status->statusPtr, &ak, params);
   CHECKSTATUSPTR(status);
   LALInspiralChooseModel(status->statusPtr, &func, &ak, params);
   CHECKSTATUSPTR(status);
 
   m = ak.totalmass = params->mass1+params->mass2;
 
   values.length = dvalues.length = valuesNew.length =
   yt.length = dym.length = dyt.length = number_of_diff_equations;
   dummy.length = number_of_diff_equations * 6;
   if (!(dummy.data = (REAL8 * ) LALMalloc(sizeof(REAL8) * number_of_diff_equations * 6))) {
      ABORT(status, LALINSPIRALH_EMEM, LALINSPIRALH_MSGEMEM);
   }
 
   values.data = &dummy.data[0];
   dvalues.data = &dummy.data[number_of_diff_equations];
   valuesNew.data = &dummy.data[2*number_of_diff_equations];
   yt.data = &dummy.data[3*number_of_diff_equations];
   dym.data = &dummy.data[4*number_of_diff_equations];
   dyt.data = &dummy.data[5*number_of_diff_equations];
 
/*   m = ak.totalmass;*/
   dt = 1./params->tSampling;
 
   if (a)
   {
     /*to be implemented. */
 
 
/*      mTot   = params->mass1 + params->mass2;
      etab   = params->mass1 * params->mass2;
      etab  /= mTot;
      etab  /= mTot;
      unitHz = mTot *LAL_MTSUN_SI*(REAL8)LAL_PI;
      cosI   = cos( params->inclination );
      mu     = etab * mTot;
      fFac   = 1.0 / ( 4.0*LAL_TWOPI*LAL_MTSUN_SI*mTot );
      f2aFac = LAL_PI*LAL_MTSUN_SI*mTot*fFac;
      apFac  = acFac = -2.0 * mu * LAL_MRSUN_SI/params->distance;
      apFac *= 1.0 + cosI*cosI;
      acFac *= 2.0*cosI;
      params->nStartPad = 0;
      */
   }
 
   ASSERT(ak.totalmass > 0, status, LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);
 
   t = 0.0;
 
 
 
   in3.totalMass = (params->mass1 + params->mass2) * LAL_MTSUN_SI ;
   in3.eta = (params->mass1 * params->mass2) /(params->mass1 + params->mass2) / (params->mass1 + params->mass2);
   funcParams = (void *) &in3;
 
 
   piM = LAL_PI * m * LAL_MTSUN_SI;
/*   f = (v*v*v)/piM;
 
   fu = params->fCutoff;
   if (fu)
      fHigh = (fu < ak.flso) ? fu : ak.flso;
   else
      fHigh = ak.flso;
   f = (v*v*v)/(LAL_PI*m);
 
   ASSERT(fHigh < 0.5/dt, status, LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);
   ASSERT(fHigh > params->fLower, status, LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);
*/
 
 
   /* e0 is set at f_min */
   e0 = params->eccentricity;
 
   /* the second harmonic will start at fLower*2/3 */
   f_min = params->fLower;
   iota = params->inclination; /*overwritten later */
 
   beta = 0.;
   twoBeta = 2.* beta;
   cos2Beta = cos(twoBeta);
   sin2Beta = sin(twoBeta);
   iota = LAL_PI/4.;
   onepCosSqI = 1. + cos(iota) * cos(iota);
   SinSqI = sin(iota) * sin(iota);
   cosI = cos(iota);
 
   p0 = (1. - e0*e0)/pow(2. * LAL_PI * m * LAL_MTSUN_SI* f_min/3. , 2./3.);
 
   *(values.data) = orbital_element_p = p0;
   *(values.data+1) = phase = params->startPhase;
   *(values.data+2) = orbital_element_e = e0;
 
 
 
 
   in4.function = LALInspiralEccentricityDerivatives;
   in4.x = t;
   in4.y = &values;
   in4.h = dt;
   in4.n = number_of_diff_equations;
   in4.yt = &yt;
   in4.dym = &dym;
   in4.dyt = &dyt;
 
   xlalErrno = 0;
   /* Initialize GSL integrator */
   if (!(integrator = XLALRungeKutta4Init(number_of_diff_equations, &in4)))
   {
     INT4 errNum = XLALClearErrno();
     LALFree(dummy.data);
 
     if (errNum == XLAL_ENOMEM)
       ABORT(status, LALINSPIRALH_EMEM, LALINSPIRALH_MSGEMEM);
     else
       ABORTXLAL( status );
   }
 
   count = 0;
   if (signalvec2) {
   params->nStartPad = 0;} /* for template generation, that value must be zero*/
 
   else if (signalvec1) {
     count = params->nStartPad;
   }
 
   t = 0.0;
 
 
   fu = params->fCutoff;
   if (fu)
      fHigh = (fu < ak.flso) ? fu : ak.flso;
   else
      fHigh = ak.flso;
 
   f = 1./(pow(orbital_element_p, 3./2.))/piM;
 
 
   /*fprintf(stderr, "fFinal = %f %f %f %f\n", fu,fHigh,f,ak.flso);*/
 done = 0;
   do {
      /* Free up memory and abort if writing beyond the end of vector*/
      /*if ((signalvec1 && (UINT4)count >= signalvec1->length) || (ff && (UINT4)count >= ff->length))*/
      if ((signalvec1 && (UINT4)count >= signalvec1->length))
      {
          XLALRungeKutta4Free( integrator );
          LALFree(dummy.data);
          ABORT(status, LALINSPIRALH_EVECTOR, LALINSPIRALH_MSGEVECTOR);
      }
 
      /* Non-injection case */
      if (signalvec1)
      {
        twoPhim2Beta = 2.* phase - twoBeta;
        phim2Beta = phase - twoBeta;
        threePhim2Beta = 3.* phase - twoBeta;
        orbital_element_e_squared = orbital_element_e * orbital_element_e;
        amp = params->signalAmplitude / orbital_element_p;
 
/*        fprintf(stderr, "%e %e %e %e %e\n", twoBeta, twoPhim2Beta, phim2Beta, threePhim2Beta, orbital_element_e_squared);*/
 
 
        h1 = amp * ( ( 2. * cos(twoPhim2Beta) + 2.5 * orbital_element_e * cos(phim2Beta)
          + 0.5 * orbital_element_e * cos(threePhim2Beta) + orbital_element_e_squared * cos2Beta) * onepCosSqI +
          + ( orbital_element_e * cos(orbital_element_p) + orbital_element_e_squared) * SinSqI);
        if ((f >= params->fLower) && (done == 0))
        {
        /*fprintf(stderr, "freq=%e p=%e, e=%e, phase = %e\n", f,orbital_element_p, orbital_element_e, phase);fflush(stderr);
*/
        params->alpha1 =  orbital_element_e;
        done = 1;
         }
         /*if (f>=params->fLower)*/
        {
          *(signalvec1->data + count) = (REAL4) h1;
         }
 
         if (signalvec2)
         {
            h2 = amp * ( ( 4. * sin(twoPhim2Beta) + 5 * orbital_element_e * sin(phim2Beta)
          + orbital_element_e * sin(threePhim2Beta) - 2. * orbital_element_e_squared * sin2Beta) * cosI);
/*           if (f>=params->fLower)*/
           {
              *(signalvec2->data + count) = (REAL4) h2;
            }
         }
      }
 
      /* Injection case */
      else if (a)
      {
        /*to be done*/
        /*
        omega = v*v*v;
 
          ff->data[count]       = (REAL4)(omega/unitHz);
          f2a                   = pow (f2aFac * omega, 2./3.);
          a->data[2*count]      = (REAL4)(4.*apFac * f2a);
          a->data[2*count+1]    = (REAL4)(4.*acFac * f2a);
          phi->data[count]      = (REAL8)(p);
          */
      }
 
      LALInspiralEccentricityDerivatives(&values, &dvalues,funcParams);
      CHECKSTATUSPTR(status);
 
      in4.dydx = &dvalues;
      in4.x=t;
 
      LALRungeKutta4(status->statusPtr, &valuesNew, integrator, funcParams);
      CHECKSTATUSPTR(status);
 
      *(values.data) = orbital_element_p = *(valuesNew.data);
      *(values.data+1) = phase = *(valuesNew.data+1);
      *(values.data+2) = orbital_element_e = *(valuesNew.data+2);
 
      t = (++count-params->nStartPad) * dt;
      /* v^3 is equal to p^(3/2)*/
      f = 1./(pow(orbital_element_p, 3./2.))/piM;
/*      fprintf(stderr, "p=%e, e=%e, phase = %e\n", orbital_element_p, orbital_element_e, phase);
      fflush(stderr);*/
      rbyM = orbital_element_p/(1.+orbital_element_e * cos(phase));
      /*fprintf(stderr, "rbyM=%e rbyMFlso=%e t=%e, ak.tn=%e f=%e e=%e\n",rbyM, rbyMFlso, t, ak.tn,f,orbital_element_e);
      fflush(stderr);*/
 
   } while ( (t < ak.tn) && (rbyM>rbyMFlso) && (f<fHigh));
 
 
   /*fprintf(stderr, "t=%e ak.tn=%e rbyM=%e rbyMFlso=%e f=%f fHigh=%f", t, ak.tn, rbyM, rbyMFlso, f, fHigh);*/
 
   /*also need to add for the fupper nyquist frequency**/
   params->vFinal = orbital_element_p;
   params->fFinal = f;
   params->tC = t;
   *countback = count;
 
   XLALRungeKutta4Free( integrator );
   LALFree(dummy.data);
 
   DETATCHSTATUSPTR(status);
   RETURN (status);
 
}
 
 
 
 
void
LALInspiralEccentricityDerivatives (
   REAL8Vector *values,
   REAL8Vector *dvalues,
   void        *params
   )
 {
 
  ecc_CBC_ODE_Input *par;
  double M, eta, mu, c1, e0, e2, UNUSED p0, p2, p3, p4, phi;
  par = (ecc_CBC_ODE_Input *) params;
 
  /* affectation */
  M = par->totalMass; /* in MTSUN*/
  eta = par->eta;  /*dimensionless*/
  mu = eta * M; /*therefore in MTSUN*/
  phi = values->data[1];
  e0 = values->data[2];
  e2 = values->data[2] * values->data[2];
  p0 = values->data[0];
  p2 = values->data[0] * values->data[0];
  p3 = p2 * values->data[0];
  p4 = p3 * values->data[0];
 
  c1 = pow(1. - e2, 1.5);
 
  /* y[0] is p
     y[1] is phase
   * y[2] is e
   */
 
/*  fprintf(stderr, "before=%e %e %e\n", p0, e0, phi);*/
  /* Eq 6 */
  dvalues->data[1] = 1./(M * sqrt (p3)) * ( 1. + e0 * cos(phi) ) * ( 1. + e0 * cos(phi) );
 
  /* Eq 8 */
  dvalues->data[0] = (-64./5.) * mu / M / M * c1 / p3 * (1 + 7./8. * e2);
 
  /* Eq 9*/
  dvalues->data[2] = (-304./15.) * mu / M / M / p4 * c1 * e0 * (1. + 121./304. * e2);
 
/*  fprintf(stderr, "new values p=%e, e=%e, phase = %e\n",
      dvalues->data[0],
      dvalues->data[2],
      dvalues->data[1]);*/
 
  return;
}