LALEOBWaveform.c

Go to the documentation of this file.

/*
*  Copyright (C) 2007 Stas Babak, David Churches, Duncan Brown, David Chin, Jolien Creighton, 
*                     B.S. Sathyaprakash, Craig Robinson , Thomas Cokelaer, Evan Ochsner
*
*  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
*/

/*  <lalVerbatim file="LALEOBWaveformCV">
Author: Sathyaprakash, B. S., Cokelaer T.
$Id: LALEOBWaveform.c,v 1.55 2008/05/29 17:08:13 spxcar Exp $
</lalVerbatim>  */

/*  <lalLaTeX>

\subsection{Module \texttt{LALEOBWaveform.c} and 
\texttt{LALEOBWaveformTemplates.c}}

Module to generate effective-one-body waveforms.

\subsubsection*{Prototypes}
\vspace{0.1in}
\input{LALEOBWaveformCP}
\index{\verb&LALEOBWaveform()&}
\begin{itemize}
\item {\tt signal:} Output containing the inspiral waveform.
\item {\tt params:} Input containing binary chirp parameters.
\end{itemize}

\input{LALEOBWaveformTemplatesCP}
\index{\verb&LALEOBWaveformTemplates()&}
\begin{itemize}
\item {\tt signal1:} Output containing the 0-phase inspiral waveform.
\item {\tt signal2:} Output containing the $\pi/2$-phase inspiral waveform.
\item {\tt params:} Input containing binary chirp parameters.
\end{itemize}

\input{LALEOBWaveformForInjectionCP}
\index{\verb&LALEOBWaveformForInjection()&}
\begin{itemize}
\item {\tt inject\_hc:} Output containing the 0-phase inspiral waveform.
\item {\tt inject\_hp:} Output containing the $\pi/2$-phase inspiral waveform.
\item {\tt inject\_phase:} Output containing the phase of inspiral waveform.
\item {\tt inject\_freq:} Output containing the frequency of inspiral waveform.
\item {\tt params:} Input containing binary chirp parameters.
\end{itemize}


\subsubsection*{Description}
By solving four coupled ordinary differential equations in
Eq.~(\ref{eq:3.28})-(\ref{3.31}) this module computes the
waveform in Eq.~(\ref{4.1}) (see discussion in Sec.~\ref{sec:EOB}
for details on how the initial conditions are chosen, when the
waveform is terminated and so on). 
No quasi-normal mode oscillations are added to the plunge signal
so the waveform is terminated around $2.8\,M$.
\subsection*{3PN vs 2PN}
At 3PN, two additional parameters exist namely OmegaS and Zeta2. 
The first parameters should be set to zero. If the  second parameter
is also set to zero then the waveform correponds to the standard 
waveforms. 
\subsubsection*{Algorithm}
A fourth order Runge-Kutta is used to solve the differential equations.

\subsubsection*{Uses}
\begin{verbatim}
   LALInspiralSetup 
   LALInspiralChooseModel
   LALInspiralVelocity
   LALInspiralPhasing1
   LALDBisectionFindRoot
   LALRungeKutta4
   LALHCapDerivatives
   LALHCapDerivatives3PN
   LALHCapDerivativesP4PN
   LALlightRingRadius
   LALlightRingRadius3PN
   LALlightRingRadiusP4PN
   LALpphiInit
   LALpphiInit3PN
   LALpphiInitP4PN
   LALprInit
   LALprInit3PN
   LALprInitP4PN
   LALrOfOmega
   LALrOfOmega3PN
   LALrOfOmegaP4PN
\end{verbatim}

\subsubsection*{Notes}
The length of the waveform returned by {\tt LALInspiralWaveLength} is 
occassionally smaller than what is required to hold an EOB waveform.
This is because EOB goes beyond the last stable orbit up to the light
ring while {\tt LALInspiralWaveLength} assumes that the waveform terminates
at the last stable orbit. It is recommended that a rather generous 
{\tt params->nEndPad} be used to prevent the code from crashing.

\vfill{\footnotesize\input{LALEOBWaveformCV}}

</lalLaTeX>  */
#include <lal/Units.h>
#include <lal/LALInspiral.h>
#include <lal/FindRoot.h>
#include <lal/SeqFactories.h>
#include <lal/NRWaveInject.h>

typedef struct tagrOfOmegaIn {
   REAL8 eta, omega;
} rOfOmegaIn;

typedef struct tagPr3In {
  REAL8 eta, zeta2, omegaS, omega, vr,r,q;
  InspiralDerivativesIn in3copy; 
} pr3In;

static void 
omegaofr3PN (
             REAL8 *x,
             REAL8 r, 
             void *params) ;

static void
omegaofrP4PN (
             REAL8 *x,
             REAL8 r,
             void *params) ;

static 
void LALHCapDerivatives(        REAL8Vector *values, 
                                REAL8Vector *dvalues, 
                                void            *funcParams);

static 
void LALprInit( REAL8           *pr, 
                                REAL8           r, 
                                InspiralDerivativesIn   *ak);

static 
void LALpphiInit(       REAL8 *phase, 
                        REAL8 r, 
                        REAL8 eta);

static 
void LALlightRingRadius(        LALStatus       *status,
                                REAL8           *x, 
                                REAL8           r, 
                                void    *params);
                                                        
static void LALrOfOmega (       LALStatus       *status, 
                                REAL8           *x,
                                REAL8           r, 
                                void            *params);

static
void LALHCapDerivatives3PN(     REAL8Vector     *values,
                                REAL8Vector     *dvalues, 
                                void                    *funcParams);
                                                        
static
void LALprInit3PN(LALStatus *status, REAL8 *pr , REAL8 , void  *params);

static
void LALpphiInit3PN(REAL8 *phase, REAL8 r, REAL8 eta, REAL8 omegaS);

static
void LALlightRingRadius3PN(LALStatus *status, REAL8 *x, REAL8 r, void *params);

static
void LALrOfOmega3PN (LALStatus *status, REAL8 *x, REAL8 r, void *params);

static
void LALvr3PN(REAL8 *vr, void *params);

static
void LALHCapDerivativesP4PN(     REAL8Vector     *values,
                                REAL8Vector     *dvalues,
                                void                    *funcParams);

static
void LALprInitP4PN(LALStatus *status, REAL8 *pr , REAL8 , void  *params);

static
void LALpphiInitP4PN(REAL8 *phase, REAL8 r, REAL8 eta, REAL8 omegaS);

static
void LALlightRingRadiusP4PN(LALStatus *status, REAL8 *x, REAL8 r, void *params);

static
void LALrOfOmegaP4PN (LALStatus *status, REAL8 *x, REAL8 r, void *params);

static
void LALvrP4PN(REAL8 *vr, void *params);

static void
LALEOBWaveformEngine (
                LALStatus        *status,
                REAL4Vector      *signal1,
                REAL4Vector      *signal2,
                REAL4Vector      *h,
                REAL4Vector      *a,
                REAL4Vector      *ff,
                REAL8Vector      *phi,
                UINT4            *countback,
                InspiralTemplate *params,
                InspiralInit     *paramsInit
                );

NRCSID (LALEOBWAVEFORMC, 
"$Id: LALEOBWaveform.c,v 1.55 2008/05/29 17:08:13 spxcar Exp $");

/*--------------------------------------------------------------------*/

static void 
LALprInit(
   REAL8 *pr, 
   REAL8 r, 
   InspiralDerivativesIn *ak
   ) 
{ 

/* 
        This combines Eq. (4.13), (4.14) (4.16) of BD2 to get 
        the initial value of pr 
*/

   REAL8        eta,
                        z,
                        omega,
                        jadiab,
                        djadiabdr, 
                        H0cap, 
                        v,
                        FDIS,
                        A,
                        r2,  /* temp variable */
                        r3,  /* temp variable */
                        cr;

   eta = ak->coeffs->eta;
   r2 = r*r;
   r3 = r2*r;

   A = 1. - 2./r + 2.*eta/r3;
   z = r3 * A * A/(r3 - 3.*r2 + 5.*eta);
   omega = sqrt((1.-3*eta/r2)/(r3 * (1. + 2*eta*(sqrt(z)-1.))));
   jadiab = sqrt (r2 * (r2 - 3.*eta)/(r3 - 3.*r2 + 5.*eta));
   djadiabdr = (r3*r3 - 6.*r3*r2 + 3.*eta*r2*r2 +20.*eta*r3-30.*eta*eta*r)/
               (pow(r3 - 3.*r2 + 5.*eta, 2.)*2.*jadiab);
   H0cap = sqrt(1. + 2.*eta*(-1. + sqrt(z)))/eta;
   cr = A*A/((1.-6.*eta/r2) * eta * H0cap * sqrt(z));
   v = pow(omega,oneby3);
   FDIS = -ak->flux(v, ak->coeffs)/(eta * omega); 
   *pr = FDIS/(djadiabdr*cr);
}

/*--------------------------------------------------------------------*/
static void 
LALpphiInit(
   REAL8 *phase, 
   REAL8 r, 
   REAL8 eta
   )
{ 
   REAL8 r2, r3;
   r2 = r*r;
   r3 = r2*r;
   *phase = pow(r2 * (r2 - 3.*eta) / (r3 - 3.* r2 + 5.*eta), 0.5);
}

/*--------------------------------------------------------------------*/
NRCSID (LALROFOMEGAC,
"$Id: LALEOBWaveform.c,v 1.55 2008/05/29 17:08:13 spxcar Exp $");

static void 
LALrOfOmega (
   LALStatus *status, 
   REAL8 *x, 
   REAL8 r, 
   void *params
   ) 
{ 

   REAL8 r2, r3, A, z, eta, omega;
   rOfOmegaIn *rofomegain;

   INITSTATUS(status, "LALrOfOmega", LALROFOMEGAC);
   ATTATCHSTATUSPTR(status);
   rofomegain = (rOfOmegaIn *) params;
   eta = rofomegain->eta;
   omega = rofomegain->omega;
   r2 = r*r;
   r3 = r2*r;

   A = 1. - 2./r + 2.*eta/r3;
   z = r3 * A * A/(r3 - 3.*r2 + 5.*eta);
   *x = omega - sqrt((1.-3*eta/r2)/(r3 * (1. + 2*eta*(sqrt(z)-1.))));
   DETATCHSTATUSPTR(status);
   RETURN(status);
}

/*--------------------------------------------------------------------*/
NRCSID (LALLIGHTRINGRADIUSC, 
"$Id: LALEOBWaveform.c,v 1.55 2008/05/29 17:08:13 spxcar Exp $");
static void 
LALlightRingRadius(
   LALStatus *status, 
   REAL8 *x, 
   REAL8 r, 
   void *params
   ) 
{ 
   REAL8 eta;
   rOfOmegaIn *rofomegain;

   INITSTATUS(status, "LALlightRingRadius", LALLIGHTRINGRADIUSC);
   ATTATCHSTATUSPTR(status);
   rofomegain = (rOfOmegaIn *) params;
   eta = rofomegain->eta;
   *x = pow(r,3.) - 3.*r*r + 5.* eta;
   DETATCHSTATUSPTR(status);
   RETURN(status);
}


static void 
LALHCapDerivatives(
   REAL8Vector *values, 
   REAL8Vector *dvalues, 
   void *funcParams
   ) 
{ 
   REAL8 r, s, p, q, r2, r3, p2, q2, A, B, dA, dB, hcap, Hcap, etahH;
   REAL8 omega, v, eta;
   InspiralDerivativesIn *ak;

   ak = (InspiralDerivativesIn *) funcParams;

   eta = ak->coeffs->eta;

   r = values->data[0];
   s = values->data[1];
   p = values->data[2];
   q = values->data[3];

   r2 = r*r;
   r3 = r2*r;
   p2 = p*p;
   q2 = q*q;
   A = 1. - 2./r + 2.*eta/r3;
   B = (1. - 6.*eta/r2)/A;
   dA = 2./r2 - 6.*eta/pow(r,4.);
   dB = (-dA * B + 12.*eta/r3)/A;
   hcap = pow (A*(1. + p2/B + q2/r2), 0.5);
   Hcap = pow (1. + 2.*eta*(hcap - 1.), 0.5) / eta;
   etahH = eta*hcap*Hcap;

   dvalues->data[0] = A*A*p/((1. - 6.*eta/r2) * etahH);
   dvalues->data[1] = omega = A * q / (r2 * etahH);

   v = pow(omega,oneby3);

   dvalues->data[2] = -0.5 * (dA * hcap * hcap/A - p2 * A * dB/(B*B) 
                    - 2. * A * q2/r3) / etahH;
   dvalues->data[3] = -ak->flux(v, ak->coeffs)/(eta * omega); 
   /*
   printf("%e %e %e %e %e %e %e %e\n", r, s, p, q, A, B, hcap, Hcap);
   */
}



/*--------------------------------------------------------------------*/

  void
LALpphiInit3PN(
            REAL8 *phase,
            REAL8 r,
            REAL8 eta,
            REAL8 omegaS
            )
{
  REAL8 u, u2, u3,  a4, a4p4eta, a4peta2, NA, DA, A, dA;
  

  u = 1./r;
  u2 = u*u;
  u3 = u2*u;
  a4 = (ninty4by3etc - 2. * omegaS) * eta;
  a4p4eta = a4 + 4. * eta;
  a4peta2 = a4 + eta * eta;
  NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
  DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
  A = NA/DA;
  dA = ( (a4 - 16. + 8. * eta) * DA - NA * (a4p4eta + 4. * a4p4eta * u 
                        + 12. * a4peta2  * u2))/(DA*DA);

  *phase = sqrt(-dA/(2.*u*A + u2 * dA));

}
/*---------------------------------------------------------------*/

NRCSID (LALPRINIT3PN, 
"$Id: LALEOBWaveform.c,v 1.55 2008/05/29 17:08:13 spxcar Exp $");
 void 
LALprInit3PN(
             LALStatus *status, 
             REAL8 *pr,
             REAL8 p,
             void *params
             ) 
{
  REAL8   u, u2, u3, u4, p2, p3, p4, q2, A, DA, NA;
  REAL8  onebyD, AbyD, Heff, HReal, etahH;
  REAL8 omegaS, eta, a4, a4p4eta, a4peta2, z3, r, vr, q;
  pr3In *ak;
  
  INITSTATUS(status, "LALprInit3PN", LALPRINIT3PN);
  ATTATCHSTATUSPTR(status);
  ak = (pr3In *) params;
  
  eta = ak->eta;
  vr = ak->vr;
  r = ak->r;
  q = ak->q;
  omegaS = ak->omegaS;
  
     
   p2 = p*p;
   p3 = p2*p;
   p4 = p2*p2;
   q2 = q*q;
   u = 1./ r;
   u2 = u*u;
   u3 = u2 * u;
   u4 = u2 * u2;
   z3 = 2. * (4. - 3. * eta) * eta;
   a4 = (ninty4by3etc - 2. * omegaS) * eta;
   a4p4eta = a4 + 4. * eta;
   a4peta2 = a4 + eta * eta;  

/* From DJS 14, 15 */
   NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
   DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
   A = NA/DA;
   onebyD = 1. + 6.*eta*u2 + 2. * ( 26. - 3. * eta) * eta * u3;
   AbyD = A * onebyD;   

   Heff = pow (A*(1. + AbyD * p2 + q*q * u2 + z3 * p4 * u2), 0.5);
   HReal = pow (1. + 2.*eta*(Heff - 1.), 0.5) / eta;
   etahH = eta*Heff*HReal;
   
   *pr = -vr +  A*(AbyD*p + 2. * z3 * u2 * p3)/etahH; 
   
   DETATCHSTATUSPTR(status);
   RETURN(status);
}

static void 
omegaofr3PN (
             REAL8 *x,
             REAL8 r, 
             void *params) 
{
   REAL8 u, u2, u3, a4, a4p4eta, a4peta2, eta, NA, DA, A, dA;
   REAL8   omegaS;

   /*include a status here ?*/
   pr3In *ak;
   ak = (pr3In *) params;
   omegaS = ak->omegaS;
   eta = ak->eta;

   u = 1./r;
   u2 = u*u;
   u3 = u2*u;
   a4 = (ninty4by3etc - 2. * omegaS) * eta;
   
   a4p4eta = a4 + 4. * eta;
   a4peta2 = a4 + eta * eta;
   NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
   DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
   A = NA/DA;
   dA = ( (a4 - 16. + 8. * eta) * DA - NA 
      * (a4p4eta + 4. * a4p4eta * u + 12. * a4peta2  * u2))/(DA*DA);
   *x = pow(u,1.5) * sqrt ( -0.5 * dA /(1. + 2.*eta * (A/sqrt(A+0.5 * u*dA)-1.)));

}

void 
LALrOfOmega3PN(
            LALStatus *status, 
            REAL8 *x, 
            REAL8 r, 
            void *params)
{
  REAL8  omega1,omega2,eta ;
  pr3In *pr3in;
       
  status = NULL;
  pr3in = (pr3In *) params;
  eta = pr3in->eta;

  omega1 = pr3in->omega;
  omegaofr3PN(&omega2,r, params);
  *x = -omega1 + omega2;

}
/*--------------------------------------------------------------------*/
NRCSID (LALLIGHTRINGRADIUS3PNC,
"$Id: LALEOBWaveform.c,v 1.55 2008/05/29 17:08:13 spxcar Exp $");
 void 
LALlightRingRadius3PN(
                      LALStatus *status, 
                      REAL8 *x, 
                      REAL8 r, 
                      void *params
                      ) 
{ 
  REAL8 eta, u, u2, u3, a4, a4p4eta,a4peta2, NA, DA, A, dA;
  rOfOmegaIn *rofomegain;
  status = NULL;
  rofomegain = (rOfOmegaIn *) params;
  eta = rofomegain->eta;


  u = 1./r;
  u2 = u*u;
  u3 = u2*u;
  a4 = ninty4by3etc * eta;
  a4p4eta = a4 + 4. * eta;
  a4peta2 = a4 + eta * eta;
  NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
  DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
  A = NA/DA;
  dA = ( (a4 - 16. + 8. * eta) * DA - NA * (a4p4eta + 4. 
        * a4p4eta * u + 12. * a4peta2  * u2))/(DA*DA);
  *x = 2 * A + dA * u;
}
/*--------------------------------------------------------------------*/

 void 
LALHCapDerivatives3PN(
                  REAL8Vector *values,
                  REAL8Vector *dvalues,
                  void *funcParams
                  )
{
   REAL8 r, s, p, q, u, u2, u3, u4, p2, p3, p4, q2, Apot, DA, NA;
   REAL8  dA, onebyD, DonebyD, AbyD, Heff, HReal, etahH;
   REAL8 omega, v, eta, a4, a4p4eta, a4peta2, z2, z30, z3, zeta2;
   REAL8 n1, c1, d1, d2, d3, oneby4meta;
   REAL8    flexNonAdiab = 0;
   REAL8    flexNonCirc = 0;

   InspiralDerivativesIn *ak;

   ak = (InspiralDerivativesIn *) funcParams;
   eta = ak->coeffs->eta;
   zeta2 = ak->coeffs->zeta2;

   r = values->data[0];
   s = values->data[1];
   p = values->data[2];
   q = values->data[3];

   p2 = p*p;
   p3 = p2*p;
   p4 = p2*p2;
   q2 = q*q;
   u = 1./r;
   u2 = u*u;
   u3 = u2 * u;
   u4 = u2 * u2;
   z30 = 2.L * (4.L - 3.L * eta) * eta;
   z2 = 0.75L * z30 * zeta2,
   z3 = z30 * (1.L - zeta2);

   a4 = ninty4by3etc * eta;
   a4p4eta = a4 + 4. * eta;
   a4peta2 = a4 + eta * eta;
   
   /* From DJS 3PN Hamiltonian */
   oneby4meta = 1./(4.-eta);
   n1 = 0.5 * (a4 - 16. + 8.*eta) * oneby4meta;
   d1 = 0.5 * a4p4eta * oneby4meta;
   d2 = a4p4eta * oneby4meta;
   d3 = 2. * a4peta2 * oneby4meta;
   NA = 1. + n1 * u;
   DA = 1 + d1*u + d2*u2 + d3*u3;
   Apot = NA/DA;
   
   onebyD = 1. + 6.*eta*u2 + (2. * ( 26. - 3. * eta) * eta - z2)* u3;
   AbyD = Apot * onebyD;
   Heff = pow (Apot*(1. + AbyD * p2 + q*q * u2 + z30 * (p4 + zeta2*(-0.25*p4
        + 0.75  * p2 * q2 * u2 )) * u2), 0.5);
   HReal = pow (1. + 2.*eta*(Heff - 1.), 0.5) / eta;

   dA = -u2/(DA*DA) * (n1*DA - NA * (d1 + 2.*d2*u + 3.*d3*u2));

   DonebyD = -12.*eta*u3 - (6.*(26. - 3.*eta)*eta - z2)*u4;
   etahH = eta*Heff*HReal;

   dvalues->data[0] = Apot*(AbyD*p +  z30 * u2 *(2.* p3 
              + zeta2*(-0.5*p3 + 0.75*p*q2*u2)))/etahH;
   dvalues->data[1] = omega = Apot * q * u2 * (1. + 0.75*z30*zeta2*p2*u2)/ etahH;
   v = pow(omega,oneby3);

   dvalues->data[2] = -0.5 * Apot * (dA*Heff*Heff/(Apot*Apot) - 2.*q2*u3 
              + (dA * onebyD + Apot * DonebyD) * p2
              + z30 * u3 *(-2.* p4+zeta2*(0.5*p4 - 3.0*p2*q2*u2))) / etahH;

   c1 = 1.+(u2 - 2.*u3*Apot/dA) * q2;
   dvalues->data[3] = -(1. - flexNonAdiab*c1) * (1. + flexNonCirc*p2/(q2*u2)) 
                                        * ak->flux(v,ak->coeffs)/(eta * omega);
}



/*----------------------------------------------------------------------*/
 void LALvr3PN(REAL8 *vr, void *params ) 
{
  REAL8 A, dA, d2A, NA, DA, dDA, dNA, d2DA;
  REAL8 u, u2, u3, v, x1; 
  REAL8 eta,a4, a4p4eta, a4peta2, FDIS;
  
  pr3In *pr3in;
  pr3in = (pr3In *)params;

  eta = pr3in->eta;     
  u = 1./ pr3in->r;
  
  u2 = u*u;
  u3 = u2*u;

  
  a4 = (ninty4by3etc - 2. * pr3in->omegaS) * eta;
  a4p4eta = a4 + 4. * eta;
  a4peta2 = a4 + eta * eta;
  NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
  DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
  A = NA/DA;
  dNA = (a4 - 16. + 8. * eta);
  dDA = (a4p4eta + 4. * a4p4eta * u + 12. * a4peta2  * u2);
  d2DA = 4. * a4p4eta + 24. * a4peta2 * u;
  
  dA = (dNA * DA - NA * dDA)/ (DA*DA);
  d2A = (-NA * DA * d2DA - 2. * dNA * DA * dDA + 2. * NA * dDA * dDA)/pow(DA,3.);
  v = pow(pr3in->omega,oneby3);
  FDIS = -pr3in->in3copy.flux(v, pr3in->in3copy.coeffs)/(eta* pr3in->omega);
  x1 = -1./u2 * sqrt (-dA * pow(2.* u * A + u2 * dA, 3.) ) 
                / (2.* u * dA * dA + A*dA - u * A * d2A);
  *vr = FDIS * x1;
}

/*-------------------------------------------------------------------*/
/*                      pseudo-4PN functions                         */
/*-------------------------------------------------------------------*/

void
LALpphiInitP4PN(
            REAL8 *phase,
            REAL8 r,
            REAL8 eta,
            REAL8 omegaS
            )
{
  REAL8 u, u2, u3, u4, a4, a5, eta2, NA, DA, A, dA;

  eta2 = eta*eta;
  u = 1./r;
  u2 = u*u;
  u3 = u2*u;
  u4 = u2*u2;
  a4 = (ninty4by3etc - 2. * omegaS) * eta;
  a5 = 60.;
  NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
  DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
       - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
       + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
  A = NA/DA;
  dA = ( (32. - 24.*eta - 4.*a4 - a5*eta) * DA - NA *
         ( -(2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u
          - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2
          + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3))/(DA*DA);
  *phase = sqrt(-dA/(2.*u*A + u2 * dA));
/* why is it called phase? This is initial j!? */
}

/*-------------------------------------------------------------------*/
  void
LALprInitP4PN(
             LALStatus *status,
             REAL8 *pr,
             REAL8 p,
             void *params
             )
{
  REAL8   u, u2, u3, u4, p2, p3, p4, q2, A, DA, NA;
  REAL8  onebyD, AbyD, Heff, HReal, etahH;
  REAL8 eta, eta2, a4, a5, z3, r, vr, q;
  pr3In *ak;

  INITSTATUS(status, "LALprInit3PN", LALPRINIT3PN);
  ATTATCHSTATUSPTR(status);
  ak = (pr3In *) params;

  eta = ak->eta;
  vr = ak->vr;
  r = ak->r;
  q = ak->q;
  eta2 = eta*eta;


   p2 = p*p;
   p3 = p2*p;
   p4 = p2*p2;
   q2 = q*q;
   u = 1./ r;
   u2 = u*u;
   u3 = u2 * u;
   u4 = u2 * u2;
   z3 = 2. * (4. - 3. * eta) * eta;
   a4 = ninty4by3etc * eta;
   a5 = 60.;

   NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
   DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
        - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
        + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
   A = NA/DA;
   onebyD = 1. + 6.*eta*u2 + 2. * ( 26. - 3. * eta) * eta * u3 + 36.*eta2*u4;
   AbyD = A * onebyD;

   Heff = pow (A*(1. + AbyD * p2 + q*q * u2 + z3 * p4 * u2), 0.5);
   HReal = pow (1. + 2.*eta*(Heff - 1.), 0.5) / eta;
   etahH = eta*Heff*HReal;

   *pr = -vr +  A*(AbyD*p + 2. * z3 * u2 * p3)/etahH;
/* This sets pr = dH/dpr - vr, calls rootfinder, 
   gets value of pr s.t. dH/pr = vr */
   DETATCHSTATUSPTR(status);
   RETURN(status);
}


/*-------------------------------------------------------------------*/
static void
omegaofrP4PN (
             REAL8 *x,
             REAL8 r,
             void *params)
{
   REAL8 u, u2, u3, u4, a4, a5, eta, eta2, NA, DA, A, dA;

   /*include a status here ?*/
   pr3In *ak;
   ak = (pr3In *) params;
   eta = ak->eta;
   eta2 = eta*eta;

   u = 1./r;
   u2 = u*u;
   u3 = u2*u;
   u4 = u2*u2;
   a4 = ninty4by3etc * eta;
   a5 = 60.;
   NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
   DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
        - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
        + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
   A = NA/DA;
   dA = ( (32. - 24.*eta - 4.*a4 - a5*eta) * DA - NA *
          ( -(2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u
           - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2
           + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3))/(DA*DA); 

   *x = pow(u,1.5) * sqrt ( -0.5 * dA /(1. + 2.*eta * (A/sqrt(A+0.5 * u*dA)-1.)));

}


/*-------------------------------------------------------------------*/
void
LALrOfOmegaP4PN(
            LALStatus *status,
            REAL8 *x,
            REAL8 r,
            void *params)
{
  REAL8  omega1,omega2,eta ;
  pr3In *pr3in;

  status = NULL;
  pr3in = (pr3In *) params;
  eta = pr3in->eta;

  omega1 = pr3in->omega;
  omegaofrP4PN(&omega2,r, params);
  *x = -omega1 + omega2;

}


/*-------------------------------------------------------------------*/
static void
LALlightRingRadiusP4PN(
                      LALStatus *status,
                      REAL8 *x,
                      REAL8 r,
                      void *params
                      )
{
  REAL8 eta, eta2, u, u2, u3, u4, a4, a5, NA, DA, A, dA;
  rOfOmegaIn *rofomegain;
  status = NULL;
  rofomegain = (rOfOmegaIn *) params;
  eta = rofomegain->eta;
  eta2 = eta*eta;


  u = 1./r;
  u2 = u*u;
  u3 = u2*u;
  u4 = u2*u2;
  a4 = ninty4by3etc * eta;
  a5 = 60.;
  NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
  DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2 
       - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3 
       + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
  A = NA/DA;
  dA = ( (32. - 24.*eta - 4.*a4 - a5*eta) * DA - NA * 
         ( -(2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u 
          - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2 
          + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3))/(DA*DA);
  *x = 2 * A + dA * u;
}

/*-------------------------------------------------------------------*/
 void
LALHCapDerivativesP4PN(
                  REAL8Vector *values,
                  REAL8Vector *dvalues,
                  void *funcParams
                  )
{
   REAL8 r, s, p, q, u, u2, u3, u4, u5, p2, p3, p4, q2, Apot, DA, NA;
   REAL8  dA, onebyD, DonebyD, AbyD, Heff, HReal, etahH;
   REAL8 omega, v, eta, eta2, a4, z2, z30, z3, zeta2;
   REAL8 a5, c1;
   double dr, ds, dp, dq;

   InspiralDerivativesIn *ak;

   ak = (InspiralDerivativesIn *) funcParams;
   eta = ak->coeffs->eta;
   zeta2 = ak->coeffs->zeta2;

   r = values->data[0];
   s = values->data[1];
   p = values->data[2];
   q = values->data[3];

   p2 = p*p;
   p3 = p2*p;
   p4 = p2*p2;
   q2 = q*q;
   u = 1./r;
   u2 = u*u;
   u3 = u2 * u;
   u4 = u2 * u2;
   u5 = u*u4;
   z30 = 2.L * (4.L - 3.L * eta) * eta;
   z2 = 0.75L * z30 * zeta2,
   z3 = z30 * (1.L - zeta2);
   eta2 = eta*eta;

   a4 = ninty4by3etc * eta;
   a5 = 60.;

   NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
   DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
        - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
        + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
   Apot = NA/DA; /* This A(u) assume zeta2=0 (the default value) */

   onebyD = 1. + 6.*eta*u2 + (2.*eta * ( 26. - 3.*eta ) - z2)* u3 + 36.*eta2*u4;
   AbyD = Apot * onebyD;
   Heff = pow (Apot*(1. + AbyD * p2 + q*q * u2 + z3 * p4 * u2), 0.5);
   HReal = pow (1. + 2.*eta*(Heff - 1.), 0.5) / eta;
   dA = -u2 * ( (32. - 24.*eta - 4.*a4 - a5*eta) * DA - NA *
          ( -(2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u
          - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2
          + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3))/(DA*DA);

   DonebyD = -12.*eta*u3 - (6.*eta*(26. - 3.*eta) - z2)*u4 - 144.*eta2*u5;
   etahH = eta*Heff*HReal;

   dr = dvalues->data[0] = Apot*(AbyD*p +  z30 * u2 *(2.* p3
              + zeta2*(-0.5*p3 + 0.75*p*q2*u2)))/etahH;
   ds = dvalues->data[1] = omega = Apot * q * u2 * (1. + 0.75*z30*zeta2*p2*u2)/ etahH;
   v = pow(omega,oneby3);

   dp = dvalues->data[2] = -0.5 * Apot * (dA*Heff*Heff/(Apot*Apot) - 2.*q2*u3
              + (dA * onebyD + Apot * DonebyD) * p2
              + z30 * u3 *(-2.* p4+zeta2*(0.5*p4 - 3.0*p2*q2*u2))) / etahH;
   c1 = 1.+(u2 - 2.*u3*Apot/dA) * q2;/*below:dpphi/dt = F_RR*/
   dq = dvalues->data[3] = - ak->flux(v,ak->coeffs)/(eta * omega);
   /*
   fprintf(stdout, "%e %e %e %e %e %e %e %e %e %e %e\n", r, s, p, q, Heff, v, Apot, dr, ds, dp, dq);
   */
}


/*-------------------------------------------------------------------*/
 void LALvrP4PN(REAL8 *vr, void *params )
{
  REAL8 A, dA, d2A, NA, DA, dDA, dNA, d2DA;
  REAL8 u, u2, u3, u4, v, x1;
  REAL8 eta, eta2, a4, a5, FDIS;

  pr3In *pr3in;
  pr3in = (pr3In *)params;

  eta = pr3in->eta;
  u = 1./ pr3in->r;

  u2 = u*u;
  u3 = u2*u;
  u4 = u2*u2;
  eta2 = eta*eta;


  a4 = (ninty4by3etc - 2. * pr3in->omegaS) * eta;
  a5 = 60.;
  NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
  DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
       - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
       + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
  A = NA/DA;
  dNA = (32. - 24.*eta - 4.*a4 - a5*eta);
  dDA = - (2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u
       - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2
       + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3;
  d2DA = - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)
       - 6.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u
       + 12.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u2;

  dA = (dNA * DA - NA * dDA)/ (DA*DA);
  d2A = (-NA * DA * d2DA - 2. * dNA * DA * dDA + 2. * NA * dDA * dDA)/pow(DA,3.);
  v = pow(pr3in->omega,oneby3);
  FDIS = -pr3in->in3copy.flux(v, pr3in->in3copy.coeffs)/(eta* pr3in->omega);
  x1 = -1./u2 * sqrt (-dA * pow(2.* u * A + u2 * dA, 3.) )
                / (2.* u * dA * dA + A*dA - u * A * d2A);
  *vr = FDIS * x1;
}


/*-------------------------------------------------------------------*/

/*  <lalVerbatim file="LALEOBWaveformCP"> */
void 
LALEOBWaveform (
   LALStatus        *status,
   REAL4Vector      *signal,
   InspiralTemplate *params
   ) 
{ /* </lalVerbatim> */

   UINT4 count;
   InspiralInit paramsInit;   
   INITSTATUS(status, "LALEOBWaveform", LALEOBWAVEFORMC);
   ATTATCHSTATUSPTR(status);

   ASSERT(signal,  status, 
        LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signal->data,  status, 
        LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   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);
   ASSERT(params->totalMass > 0., status, 
        LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);

   LALInspiralSetup (status->statusPtr, &(paramsInit.ak), params);
   CHECKSTATUSPTR(status);
   LALInspiralChooseModel(status->statusPtr, &(paramsInit.func),
                                         &(paramsInit.ak), params);
   CHECKSTATUSPTR(status);

   memset(signal->data, 0, signal->length * sizeof( REAL4 ));

   /* Call the engine function */
   LALEOBWaveformEngine(status->statusPtr, signal, NULL, NULL, NULL, 
                        NULL, NULL, &count, params, &paramsInit);
   CHECKSTATUSPTR( status );

   DETATCHSTATUSPTR(status);
   RETURN(status);
}


NRCSID (LALEOBWAVEFORMTEMPLATESC,
"$Id: LALEOBWaveform.c,v 1.55 2008/05/29 17:08:13 spxcar Exp $");

/*  <lalVerbatim file="LALEOBWaveformTemplatesCP"> */

void 
LALEOBWaveformTemplates (
   LALStatus        *status,
   REAL4Vector      *signal1,
   REAL4Vector      *signal2,
   InspiralTemplate *params
   ) 
{ /* </lalVerbatim> */

   UINT4 count;
  
   InspiralInit paramsInit;
 
   INITSTATUS(status, "LALEOBWaveformTemplates", LALEOBWAVEFORMTEMPLATESC);
   ATTATCHSTATUSPTR(status);
   
   ASSERT(signal1,  status,
        LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signal2,  status,
        LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signal1->data,  status, 
        LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signal2->data,  status, 
        LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   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);
   ASSERT(params->totalMass > 0., status, 
        LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);

   LALInspiralSetup (status->statusPtr, &(paramsInit.ak), params);
   CHECKSTATUSPTR(status);
   LALInspiralChooseModel(status->statusPtr, &(paramsInit.func), 
                                        &(paramsInit.ak), params);
   CHECKSTATUSPTR(status);
   
   memset(signal1->data, 0, signal1->length * sizeof( REAL4 ));
   memset(signal2->data, 0, signal2->length * sizeof( REAL4 ));

   /* Call the engine function */
   LALEOBWaveformEngine(status->statusPtr, signal1, signal2, NULL, NULL, 
                           NULL, NULL, &count, params, &paramsInit);
   CHECKSTATUSPTR( status );

   DETATCHSTATUSPTR(status);
   RETURN(status);
}


/*=========================================================*/
/*======INJECTION =========================================*/
/*=========================================================*/

/*  <lalVerbatim file="LALEOBWaveformForInjectionCP"> */
void 
LALEOBWaveformForInjection (
                            LALStatus        *status,
                            CoherentGW       *waveform,
                            InspiralTemplate *params,
                            PPNParamStruc    *ppnParams
                            ) 
{ 
  /* </lalVerbatim> */
  UINT4 count, i;

  REAL4Vector *a=NULL;/* pointers to generated amplitude  data */
  REAL4Vector *h=NULL;/* pointers to generated polarization data */
  REAL4Vector *ff=NULL ;/* pointers to generated  frequency data */
  REAL8Vector *phi=NULL;/* pointer to generated phase data */

  REAL8 s;

  REAL8 phiC;/* phase at coalescence */
  CHAR message[256];
  InspiralInit paramsInit;