lal » packages » inject » src

GeneratePPNAmpCorInspiral.c

Go to the documentation of this file.
00001 /*
00002 
00003 *  Copyright (C) 2007 David McKechan, Thomas Cokelaer
00004 *
00005 *  This program is free software; you can redistribute it and/or modify
00006 *  it under the terms of the GNU General Public License as published by
00007 *  the Free Software Foundation; either version 2 of the License, or
00008 *  (at your option) any later version.
00009 *
00010 *  This program is distributed in the hope that it will be useful,
00011 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013 *  GNU General Public License for more details.
00014 *
00015 *  You should have received a copy of the GNU General Public License
00016 *  along with with program; see the file COPYING. If not, write to the
00017 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
00018 *  MA  02111-1307  USA
00019 */
00020 
00021 /************************** <lalVerbatim file="GeneratePPNAmpCorInspiralCV">
00022 Author: Creighton, T. D., McKechan David, Van Den Broeck Chris
00023 $Id: GeneratePPNAmpCorInspiral.c,v 1.33 2008/08/13 18:32:52 mckechan Exp $
00024 **************************************************** </lalVerbatim> */
00025 
00026 /********************************************************** <lalLaTeX>
00027 
00028 \providecommand{\lessim}{\stackrel{<}{\scriptstyle\sim}}
00029 
00030 \subsection{Module \texttt{GeneratePPNampCorInspiral.c}}
00031 \label{ss:GeneratePPNAmpCorInspiral.c}
00032 
00033 Computes a parametrized post-Newtonian inspiral waveform 
00034 with ampltidude corrections.
00035 
00036 \subsubsection*{Prototypes}
00037 \vspace{0.1in}
00038 \input{GeneratePPNAmpCorInspiralCP}
00039 \idx{LALGeneratePPNampCorInspiral()}
00040 
00041 \subsubsection*{Description}
00042 
00043 See GeneratePPNInspiral.c
00044 
00045 Phase computed to 3.5PN.
00046 Amplitude computed to 2.5PN.
00047 
00048 The ampitude corrrected gravitaional wave polarizations $h_+$ and $h_x$ 
00049 are stored in output.h.
00050 
00051 Warning! output.a is used to store the first three harmonics in 
00052 alternate values, (i.e. [a1(0),a2(0),a3(0),a1(dt),a2(dt),a3(dt)...]) as 
00053 this will be used for filtering with higher harmonic waveforms. 
00054 
00055 Although $h_{+,\times} are computed, $output.phi is also required for filtering.
00056 
00057 
00058 \subsubsection*{Algorithm}
00059 
00060 
00061 \subsubsection*{Uses}
00062 \begin{verbatim}
00063 LALMalloc()                   LALFree()
00064 LALSCreateVectorSequence()    LALSDestroyVectorSequence()
00065 LALSCreateVector()            LALSDestroyVector()
00066 LALDCreateVector()            LALDDestroyVector()
00067 LALSBisectionFindRoot()       LALSnprintf()
00068 \end{verbatim}
00069 
00070 \subsubsection*{Notes}
00071 
00072 \vfill{\footnotesize\input{GeneratePPNAmpCorInspiralCV}}
00073 
00074 ******************************************************* </lalLaTeX> */
00075 
00076 #include <math.h>
00077 #include <lal/LALStdio.h>
00078 #include <lal/LALStdlib.h>
00079 #include <lal/LALConstants.h>
00080 #include <lal/Units.h>
00081 #include <lal/FindRoot.h>
00082 #include <lal/AVFactories.h>
00083 #include <lal/SeqFactories.h>
00084 #include <lal/SimulateCoherentGW.h>
00085 #include <lal/GeneratePPNInspiral.h>
00086 
00087 NRCSID( GENERATEPPNAMPCORINSPIRALC, "$Id: GeneratePPNAmpCorInspiral.c,v 1.33 2008/08/13 18:32:52 mckechan Exp $" );
00088 
00089 /* Define some constants used in this module. */
00090 #define MAXORDER 8               /* Maximum number of N and PN terms */
00091 #define AMPMAXORDER 6            /* Maximum PN order in amplitude (plus one) */
00092 #define BUFFSIZE 1024            /* Number of timesteps buffered */
00093 #define ACCURACY (1.0e-8)        /* Accuracy of root finder */
00094 #define TWOTHIRDS (0.6666666667) /* 2/3 */
00095 #define ONEMINUSEPS (0.99999)    /* Something close to 1 */
00096 
00097 /* 
00098    A macro to computing the (normalized) frequency.  It appears in
00099    many places, including in the main loop, and I don't want the
00100    overhead of a function call.  The following variables are required
00101    to be defined and set outside of the macro:
00102 
00103    REAL4 c0, c1, c2, c3, c4, c5, c6, c7, p6; PN frequency coefficients
00104    BOOLEAN b0, b1, b2, b3, b4, b5, b6, b7;   whether to include each PN term
00105 
00106    The following variables must be defined outside the macro, but are
00107    set inside it:
00108 
00109    REAL4 x2, x3, x4, x5, x6, x7;  the input x raised to power 
00110    2, 3, 4, 5, 6 and 7 
00111 */
00112 
00113 #define FREQ( f, x )                                                 \
00114 do {                                                                 \
00115   x2 = (x)*(x);                                                      \
00116   x3 = x2*(x);                                                       \
00117   x4 = x3*(x);                                                       \
00118   x5 = x4*(x);                                                       \
00119   x6 = x5*(x);                                                       \
00120   x7 = x6*(x);                                                       \
00121   (f) = 0;                                                           \
00122   if ( b0 )                                                          \
00123     (f) += c0;                                                       \
00124   if ( b1 )                                                          \
00125     (f) += c1*(x);                                                   \
00126   if ( b2 )                                                          \
00127     (f) += c2*x2;                                                    \
00128   if ( b3 )                                                          \
00129     (f) += c3*x3;                                                    \
00130   if ( b4 )                                                          \
00131     (f) += c4*x4;                                                    \
00132   if ( b5 )                                                          \
00133     (f) += c5*x5;                                                    \
00134   if ( b6 )                                                          \
00135     (f) += (c6 + p6*( -107.0/2240.0*(-8.0)*(log(2.0*x))))*x6;        \
00136   if ( b7 )                                                          \
00137     (f) += c7*x7;                                                    \
00138   (f) *= x3;                                                         \
00139 } while (0)
00140 
00141 /* Definition of a data structure used by FreqDiff() below. */
00142 typedef struct tagFreqDiffParamStruc 
00143 {
00144   REAL4 *c;   /* PN coefficients of frequency series */
00145   BOOLEAN *b; /* whether to include each PN term */
00146   REAL4 p6;   /* synonym for p[6] */
00147   REAL4 y0;   /* normalized frequency being sought */
00148 } FreqDiffParamStruc;
00149 
00150 /* A function to compute the difference between the current and
00151    requested normalized frequency, used by the root bisector. */
00152 static void
00153 FreqDiff( LALStatus *stat, REAL4 *y, REAL4 x, void *p )
00154 {
00155   FreqDiffParamStruc *par;                  /* *p cast to its proper type */
00156   REAL4 c0, c1, c2, c3, c4, c5, c6, c7, p6; /* PN frequency coefficients */
00157   BOOLEAN b0, b1, b2, b3, b4, b5, b6, b7;   /* whether each order is nonzero */
00158   REAL4 x2, x3, x4, x5, x6, x7;             /* x^2, x^3, x^4, x^5, x^6, & x^7 */
00159 
00160   INITSTATUS( stat, "FreqDiff", GENERATEPPNAMPCORINSPIRALC );
00161   ASSERT( p, stat, 1, "Null pointer" );
00162 
00163   /* Set constants used by FREQ() macro. */
00164   par = (FreqDiffParamStruc *)( p );
00165   c0 = par->c[0];
00166   c1 = par->c[1];
00167   c2 = par->c[2];
00168   c3 = par->c[3];
00169   c4 = par->c[4];
00170   c5 = par->c[5];
00171   c6 = par->c[6];
00172   c7 = par->c[7];
00173   b0 = par->b[0];
00174   b1 = par->b[1];
00175   b2 = par->b[2];
00176   b3 = par->b[3];
00177   b4 = par->b[4];
00178   b5 = par->b[5];
00179   b6 = par->b[6];
00180   b7 = par->b[7];
00181   p6 = par->p6;
00182 
00183   /* Evaluate frequency and compare with reference. */
00184   FREQ( *y, x );
00185   *y -= par->y0;
00186   RETURN( stat );
00187 }
00188 
00189 /* Definition of a data buffer list for storing the waveform. */
00190 typedef struct tagPPNInspiralBuffer 
00191 {
00192   REAL4 h[2*BUFFSIZE];               /* polarisation data */
00193   REAL4 a[3*BUFFSIZE];               /* first three harmonics data */
00194   REAL8 phi[BUFFSIZE];               /* phase data */
00195   REAL4 f[BUFFSIZE];                 /* frequency data */
00196   struct tagPPNInspiralBuffer *next; /* next buffer in list */
00197 } PPNInspiralBuffer;
00198 
00199 /* Definition of a macro to free the tail of said list, from a given
00200    node onward. */
00201 #define FREELIST( node )                                             \
00202 do {                                                                 \
00203   PPNInspiralBuffer *herePtr = (node);                               \
00204   while ( herePtr ) {                                                \
00205     PPNInspiralBuffer *lastPtr = herePtr;                            \
00206     herePtr = herePtr->next;                                         \
00207     LALFree( lastPtr );                                              \
00208   }                                                                  \
00209 } while (0)
00210 
00211 
00212 /*********************************************************************
00213  * MAIN FUNCTION                                                     *
00214  *********************************************************************/
00215 
00216 /* <lalVerbatim file="GeneratePPNAmpCorInspiralCP"> */
00217 void
00218 LALGeneratePPNAmpCorInspiral( 
00219                               LALStatus     *stat,
00220                               CoherentGW    *output,
00221                               PPNParamStruc *params 
00222                             )
00223 { /* </lalVerbatim> */
00224 
00225   /* System-derived constants. */
00226   BOOLEAN b0, b1, b2, b3, b4, b5, b6, b7; /* whether each order is nonzero */
00227   BOOLEAN b[MAXORDER];                    /* vector of above coefficients */
00228   REAL4 c0, c1, c2, c3, c4, c5, c6, c7;   /* PN frequency coefficients */
00229   REAL4 c[MAXORDER];                      /* vector of above coefficients */
00230   REAL4 d0, d1, d2, d3, d4, d5, d6, d7;   /* PN phase coefficients */
00231   REAL4 e0, e1, e2, e3, e4, e5, e6, e7;   /* PN dy/dx coefficients */
00232   REAL4 p[MAXORDER];                      /* PN parameter values in phase */
00233   REAL4 q[AMPMAXORDER];                   /* PN parameter values in amplitude */ 
00234   REAL4 p6;                        /* synonym for p[6] */
00235   UINT4 ampOrder;                  /* Amplitude Order */
00236   INT4 Harmonics;                  /* Number of harmonics */
00237   REAL4 mTot, mu;                  /* total mass and reduced mass */
00238   REAL4 eta, etaInv;               /* mass ratio and its inverse */
00239   REAL4 phiC;                      /* phase at coalescence */
00240   REAL4 cosI, cos2I, cos4I, cos6I; /* cosine of system inclination */
00241   REAL4 sinI, sin2I, sin4I, sin5I; /* sine of system inclination */
00242 
00243   REAL4 fFac;          /* SI normalization for f and t */
00244   REAL4 f2aFac;        /* factor multiplying f in amplitude function */
00245   REAL4 fthree, ffour, ffive, fsix, fseven; /* powers in f2a to speed up 
00246                                                waveform construction */
00247   REAL4 preFac;        /* Overall prefactor in waveforms */
00248   REAL4 delta;         /* relative mass difference */
00249   REAL4 sd, scd;       /* sinI*delta, sd*cosI*/
00250 
00251   /* Integration parameters. */
00252   UINT4 i;               /* index over PN terms */
00253   UINT4 j;               /* index of leading nonzero PN term */
00254   UINT4 n, nMax;         /* index over timesteps, and its maximum + 1 */
00255   UINT4 nNext;           /* index where next buffer starts */
00256   REAL8 t, t0, dt;       /* dimensionless time, start time, and increment */
00257   REAL4 tStop = 0.0625;  /* time when orbit reaches minimum radius */
00258   REAL4 x, xStart, xMax; /* x = t^(-1/8), and its maximum range */
00259   REAL4 y, yStart, yMax; /* normalized frequency and its range and start time */
00260   REAL4 yOld, dyMax;     /* previous timestep y, and maximum y - yOld */
00261   REAL4 *f;              /* pointer to generated frequency data */
00262   REAL4 x2, x3, x4, x5, x6, x7;  /* x^2, x^3, x^4, x^5, x^6 and x^7 */
00263   
00264   /* Harmonic terms */
00265   REAL4 a1Pthree, a1Pfive, a1Psix, a1Pseven, a1Pmixsix;
00266   REAL4 a2Ptwo, a2Pfour, a2Pfive, a2Psix, a2Pseven, a2Pmixseven; 
00267   REAL4 a3Pthree, a3Pfive, a3Psix, a3Pseven, a3Pmixsix;
00268   REAL4 a4Pfour, a4Psix, a4Pseven, a4Pmixseven;
00269   REAL4 a5Pfive, a5Pseven;
00270   REAL4 a6Psix;
00271   REAL4 a7Pseven;
00272   REAL4 a1Cthree, a1Cfive, a1Csix, a1Cseven, a1Cmixsix;
00273   REAL4 a2Ctwo, a2Cfour, a2Cfive, a2Csix, a2Cseven, a2Cmixseven; 
00274   REAL4 a3Cthree, a3Cfive, a3Csix, a3Cseven, a3Cmixsix;
00275   REAL4 a4Cfour, a4Csix, a4Cseven, a4Cmixseven;
00276   REAL4 a5Cfive, a5Cseven;
00277   REAL4 a6Csix;
00278   REAL4 a7Cseven;
00279 
00280   REAL4 a1, a2, a3, a4, a5, a6, a7; /* generated amplitudes of harmonics */
00281   REAL4 a1mix, a2mix, a3mix, a4mix; /* generated amplitudes of harmonic 
00282                                        mixed terms */
00283   REAL4 *h;                /* pointer to generated hplus and hcross */
00284   REAL4 *a;                /* pointer to generated first three plus harmonics*/
00285   REAL8 *phi;              /* pointer to generated phase data */
00286 
00287   PPNInspiralBuffer *head, *here; /* pointers to buffered data */
00288 
00289   INITSTATUS( stat, "LALGeneratePPNAmpCorInspiral", GENERATEPPNAMPCORINSPIRALC);
00290   ATTATCHSTATUSPTR( stat );
00291 
00292   /*******************************************************************
00293    * CHECK INPUT PARAMETERS                                          *
00294    *******************************************************************/
00295 
00296   /* Dumb initialization to shut gcc up. */
00297   head = here = NULL;
00298   b0 = b1 = b2 = b3 = b4 = b5 = b6 = b7 = 0.0;
00299   c0 = c1 = c2 = c3 = c4 = c5 = c6 = c7 = 0.0;
00300   d0 = d1 = d2 = d3 = d4 = d5 = d6 = d7 = 0.0;
00301 
00302   /* Make sure parameter and output structures exist. */
00303   ASSERT( params, stat, GENERATEPPNINSPIRALH_ENUL,
00304     GENERATEPPNINSPIRALH_MSGENUL );
00305   ASSERT( output, stat, GENERATEPPNINSPIRALH_ENUL,
00306     GENERATEPPNINSPIRALH_MSGENUL );
00307 
00308   /* Make sure output fields don't exist. */
00309   ASSERT( !( output->h ), stat, GENERATEPPNINSPIRALH_EOUT,
00310     GENERATEPPNINSPIRALH_MSGEOUT );
00311   ASSERT( !( output->a ), stat, GENERATEPPNINSPIRALH_EOUT,
00312     GENERATEPPNINSPIRALH_MSGEOUT );
00313   ASSERT( !( output->f ), stat, GENERATEPPNINSPIRALH_EOUT,
00314     GENERATEPPNINSPIRALH_MSGEOUT );
00315   ASSERT( !( output->phi ), stat, GENERATEPPNINSPIRALH_EOUT,
00316     GENERATEPPNINSPIRALH_MSGEOUT );
00317   ASSERT( !( output->shift ), stat, GENERATEPPNINSPIRALH_EOUT,
00318     GENERATEPPNINSPIRALH_MSGEOUT );
00319 
00320   /* Get PN parameters, if they are specified; otherwise use
00321      3.5 post-Newtonian. */
00322   if ( params->ppn ) {
00323     ASSERT( params->ppn->data, stat, GENERATEPPNINSPIRALH_ENUL, 
00324         GENERATEPPNINSPIRALH_MSGENUL );
00325     j = params->ppn->length;
00326     if ( j > MAXORDER )
00327       j = MAXORDER;
00328     for ( i = 0; i < j; i++ )
00329       p[i] = params->ppn->data[i];
00330     for ( ; i < MAXORDER; i++ )
00331       p[i] = 0.0;
00332   } 
00333   else {
00334     p[0] = 1.0;
00335     p[1] = 0.0;
00336     p[2] = 1.0;
00337     p[3] = 1.0;
00338     p[4] = 1.0;
00339     p[5] = 1.0;
00340     for ( i = 6; i < MAXORDER; i++ )
00341       p[i] = 1.0;
00342   }
00343 
00344 
00345   
00346   /* Set PN parameters for amplitude */
00347   if ( (params->ampOrder < 0) || (params->ampOrder >= AMPMAXORDER) )
00348   {    
00349     if (params->ppn->length - 1 >= 5 )
00350       ampOrder = 5;
00351     else
00352       ampOrder = params->ppn->length - 1;      
00353   }
00354   else
00355     ampOrder = params->ampOrder;    
00356   
00357   q[0] = 1.0;
00358   for(i = 1; i < AMPMAXORDER; i++)
00359     q[i] = ( i <= ampOrder? 1.0 : 0.0 );
00360   
00361   /* Set number of harmonics in accordance with params->ampOrder*/
00362   /* Dominant harmonic is the 2nd */
00363   Harmonics = ampOrder + 2;
00364 
00365 
00366   /*******************************************************************
00367    * COMPUTE SYSTEM PARAMETERS                                       *
00368    *******************************************************************/
00369 
00370   /* Compute parameters of the system. */
00371   mTot = params->mTot;
00372   ASSERT( mTot != 0.0, stat, GENERATEPPNINSPIRALH_EMBAD,
00373     GENERATEPPNINSPIRALH_MSGEMBAD );
00374   eta = params->eta;
00375   ASSERT( eta != 0.0, stat, GENERATEPPNINSPIRALH_EMBAD,
00376     GENERATEPPNINSPIRALH_MSGEMBAD );
00377   etaInv = 2.0 / eta;
00378   mu = eta*mTot;
00379 
00380   sinI = sin( params->inc );
00381   sin2I = sinI*sinI;
00382   sin4I = sin2I*sin2I;
00383   sin5I = sin4I*sinI;
00384 
00385   cosI = cos( params->inc );
00386   cos2I = cosI*cosI;
00387   cos4I = cos2I*cos2I;
00388   cos6I = cos4I*cos2I;
00389   
00390   phiC = params->phi;
00391 
00392   preFac = -2.0*mu*LAL_MRSUN_SI/params->d; 
00393   delta = pow((1-4*eta), 0.5); 
00394   sd = sinI*delta;
00395   scd = sd*cosI;
00396 
00397 
00398 
00399   /* *************************************************************
00400    * COEFFICIENTS IN h_+ & h_x ***********************************
00401    *************************************************************** */
00402 
00403   /* The variables are annotated by the harmonic, then the polarisation and 
00404      then the frequency power by which they are multiplied. For instance 
00405      'a1Pthree' is the first harmonic, plus polarisation and is multiplied 
00406      by fthree */
00407 
00408   /* First harmonic plus */
00409   a1Pthree = sd*(5.0 + cos2I)/8.0;
00410 
00411   a1Pfive = - sd*(
00412                  19.0/64.0 + 5.0/16.0*cos2I - 1.0/192.0*cos4I 
00413                  + eta*(-49.0/96.0 + 1.0/8.0*cos2I + 1.0/96.0*cos4I)
00414                  );
00415 
00416   a1Psix = + sd*LAL_PI*(5.0 + cos2I)/8.0;
00417 
00418   a1Pseven = - sd*(
00419                   (1771.0 - 1667.0*cos2I)/5120 + (217*cos4I - cos6I)/9216.0
00420                   + eta*(
00421                         681.0/256.0 + (13.0*cos2I - 35*cos4I)/768.0 
00422                         + cos6I/2304.0
00423                         )
00424                   + eta*eta*(
00425                             -(3451.0 + 5.0*cos4I)/9216.0 
00426                             + (673.0*cos2I - cos6I)/3072.0
00427                             )
00428                   );  
00429 
00430   a1Pmixsix = - sd*(
00431                    11.0/40.0 + 5.0*log(2.0)/4.0 
00432                    + cos2I*(7.0/40.0 + log(2.0)/4.0)
00433                    );
00434 
00435 
00436   
00437   /* Second harmonic plus */
00438   a2Ptwo = (1.0 + cos2I);
00439 
00440   a2Pfour = - (
00441               19.0/6.0 + 3.0/2.0*cos2I - 1.0/3.0*cos4I 
00442               + eta*(-19.0/6.0 + 11.0/6.0*cos2I + cos4I)
00443               );
00444 
00445   a2Pfive = 2.0*LAL_PI*(1.0 + cos2I);
00446 
00447   a2Psix = - (
00448              11.0/60.0 + 33.0/10.0*cos2I + (29.0*cos4I - 1.0*cos6I)/24.0 
00449              + eta*(
00450                    353.0/36.0 - 3.0*cos2I - 251.0/72.0*cos4I 
00451                    + 5.0/24.0*cos6I
00452                    )
00453              + eta*eta*(-49.0/12.0 + 9.0/2.0*cos2I 
00454              - cos4I*(7.0 + 5.0*cos2I)/24.0)
00455              );
00456 
00457   a2Pseven = - LAL_PI*(
00458                       19.0/3.0 + 3.0*cos2I - 2.0/3.0*cos4I 
00459                       + eta*((-16.0 + 14.0*cos2I)/3.0 + 2*cos4I)
00460                       );     
00461 
00462   a2Pmixseven = - (
00463                   -9.0 + 14.0*cos2I + 7.0*cos4I 
00464                   + eta*(96.0 - 8.0*cos2I - 28.0*cos4I)
00465                   )/5.0;               
00466 
00467 
00468   
00469   /* Third harmonic plus */
00470   a3Pthree = -9.0/8.0*sd*(1.0 + cos2I);
00471 
00472   a3Pfive = - sd*(
00473                  - 657.0/128.0 -45.0/16.0*cos2I + 81.0/128.0*cos4I 
00474                  + eta*(225.0/64.0 - 9.0/8.0*cos2I - 81.0/64.0*cos4I)
00475                  );
00476 
00477   a3Psix = - sd*LAL_PI*27.0/8.0*(1.0 + cos2I);
00478 
00479   a3Pseven = - sd*(
00480                   3537.0/1024.0 
00481                   - (22977*cos2I + 15309.0*cos4I - 729.0*cos6I)/5120
00482                   + eta*(
00483                         -23829.0 + 5529.0*cos2I 
00484                         + 7749.0*cos4I -729.0*cos6I
00485                         )/1280.0
00486                   + eta*eta*(
00487                             29127.0 - 27267.0*cos2I - 1647.0*cos4I 
00488                             + 2187.0*cos6I
00489                             )/5120.0
00490                   );
00491 
00492   a3Pmixsix = - sd*(-189.0/40.0 + 27.0/4.0*log(1.5))*(1 + cos2I);         
00493 
00494 
00495 
00496   /* Fourth harmonic plus */      
00497   a4Pfour = 4.0/3.0*sin2I*(1.0 + cos2I)*(1.0 - 3.0*eta);
00498 
00499   a4Psix = - (
00500              118.0/15.0 - 16.0/5.0*cos2I - cos4I*(86.0 - 16.0*cos2I)/15.0
00501              + eta*(
00502                    -262.0/9.0 + 16.0*cos2I + 166.0/9.0*cos4I 
00503                    - 16.0/3.0*cos6I
00504                    )  
00505              + eta*eta*(14.0 - 16.0*cos2I + (-10.0*cos4I + 16.0*cos6I)/3.0)
00506              );
00507   
00508   a4Pseven = + 16.0*LAL_PI/3.0*(1.0 + cos2I)*sin2I*(1.0 - 3.0*eta);
00509 
00510   a4Pmixseven = - sin2I*(1.0 + cos2I)*(
00511                                       56.0/5.0 - 32.0*log(2.0)/3.0 
00512                                       - eta*(1193.0/30.0 -32.0*log(2.0))
00513                                       );  
00514 
00515 
00516 
00517   /* Fifth harmonic plus */
00518   a5Pfive = - sd*(625.0/384.0*sin2I*(1.0 + cos2I)*(1.0 - 2.0*eta));
00519 
00520   a5Pseven = - sd*(
00521                   (-108125.0 + 40625.0*cos2I + 83125.0*cos4I 
00522                   - 15625.0*cos6I
00523                   )/9216.0
00524                   + eta*(8125.0/256.0 - (
00525                                         40625.0*cos2I + 48125.0*cos4I 
00526                                         - 15625.0*cos6I
00527                                         )/2304.0
00528                         )
00529                   + eta*eta*(
00530                             (44375.0*cos4I - 119375.0)/9216.0 
00531                             + (40625.0*cos2I - 15625*cos6I)/3072.0)
00532                             );      
00533 
00534 
00535 
00536   /* Sixth harmonic plus */
00537   a6Psix = 81.0/40.0*sin4I*(1.0 + cos2I)*(1.0 + 5.0*eta*(eta - 1.0));
00538 
00539 
00540 
00541   /* Seventh harmonic plus */
00542   a7Pseven = delta*sin5I*117649.0/46080.0*(1 + cos2I)*(1 + eta*(3.0*eta - 4.0));
00543 
00544 
00545 
00546   /* First harmonic cross */
00547   a1Cthree = 3.0/4.0*scd;
00548 
00549   a1Cfive  = - scd*(21.0/32.0 - 5.0/96.0*cos2I + eta*(-23.0 + 5.0*cos2I)/48.0);
00550 
00551   a1Csix   = scd*3.0*LAL_PI/4.0;                
00552 
00553   a1Cseven =  - scd*(-913.0/7680.0 + 1891.0/11520.0*cos2I - 7.0/4608.0*cos4I
00554                  + eta*(1165.0/384.0 - 235.0/576.0*cos2I + 7.0/1152.0*cos4I)
00555            + eta*eta*(-1301.0/4608.0 + 301.0/2304.0*cos2I - 7.0/1536.0*cos4I)); 
00556 
00557   a1Cmixsix   = scd*(9.0/20.0 + 3.0*log(2.0)/2.0);
00558 
00559 
00560        
00561   /* Second harmonic cross */
00562 
00563   a2Ctwo = 2.0*cosI;
00564 
00565   a2Cfour = - cosI*(17.0/3.0 - 4.0/3.0*cos2I + eta*(-13.0/3.0 + 4.0*cos2I));         
00566   a2Cfive = 4*LAL_PI*cosI;
00567 
00568   a2Csix = - cosI*(
00569                   17.0/15.0 + 113.0/30.0*cos2I - 0.25*cos4I
00570                   + eta*(143.0/9.0 - 245.0/18.0*cos2I + 5.0/4.0*cos4I)
00571                   + eta*eta*(-14.0/3.0 + 35.0/6.0*cos2I - 5.0/4.0*cos4I)
00572                   );
00573 
00574   a2Cseven = - LAL_PI*cosI*(
00575                            (34.0 - 8.0*cos2I)/3.0 
00576                            - eta*(20.0/3.0 - 8.0*cos2I)
00577                            );
00578 
00579   a2Cmixseven = - cosI*(2.0 + (-22.0*cos2I + eta*(-154.0 + 94.0*cos2I))/5.0);
00580 
00581 
00582    
00583   /* Third harmonic cross */   
00584   a3Cthree = - 9.0/4.0*scd;    
00585 
00586   a3Cfive = - scd*(
00587                   -603.0/64.0 + 135.0/64.0*cos2I 
00588                   + eta*(171.0 - 135.0*cos2I)/32.0
00589                   );
00590 
00591   a3Csix = - scd*27.0/4.0*LAL_PI;
00592 
00593   a3Cseven = - scd*(
00594                    (12501.0 - 24138.0*cos2I + 1701.0*cos4I)/2560.0
00595                    + eta*(-19581.0 + 15642.0*cos2I - 1701.0*cos4I)/640.0
00596                    + eta*eta*(18903.0 - 22806.0*cos2I + 5103.0*cos4I)/2560.0
00597                    );
00598 
00599   a3Cmixsix = - scd*(189.0/20.0 - 27.0/2.0*log(1.5));
00600 
00601 
00602 
00603   /* Fourth harmonic cross */
00604   a4Cfour = cosI*sin2I*8.0/3.0*(1.0 - 3.0*eta);
00605 
00606   a4Csix = - cosI*(
00607                   44.0/3.0 - 268.0/15.0*cos2I + 16.0/5.0*cos4I
00608                   + eta*((-476.0 + 620.0*cos2I)/9.0 - 16.0*cos4I)
00609                   + eta*eta*((68.0 - 116.0*cos2I)/3.0 + 16.0*cos4I)
00610                   );
00611 
00612   a4Cseven = sin2I*cosI*32.0/3.0*LAL_PI*(1.0 - 3.0*eta);
00613 
00614   a4Cmixseven = - cosI*sin2I*(
00615                              -112.0/5.0 + 64.0*log(2.0)/3.0 
00616                              + eta*(1193.0/15.0 - 64.0*log(2.0))
00617                              );  
00618 
00619 
00620 
00621   /* Fifth harmonic cross */
00622   a5Cfive = - scd*(625.0/192.0*(1.0 - 2.0*eta)*sin2I); 
00623 
00624   a5Cseven = - scd*(
00625                    6875.0/256.0*cos2I 
00626                    - (101875.0 + 21875.0*cos4I)/4608.0
00627                    + eta*(
00628                          (66875.0 + 21875.0*cos4I)/1152.0 
00629                          - 44375.0/576.0*cos2I
00630                          )
00631                    + eta*eta*(
00632                              -100625.0/4608.0 + 83125.0/2304.0*cos2I 
00633                              - 21875.0/1536.0*cos4I
00634                              )
00635                    );
00636 
00637 
00638 
00639   /* Sixth harmonic cross */
00640   a6Csix = cosI*81.0/20.0*sin4I*(1.0 + 5.0*eta*(eta - 1.0));
00641 
00642 
00643 
00644   /* Seventh harmonic cross */
00645   a7Cseven = - scd*sin4I*117649.0/23040.0*(1.0 + eta*(3.0*eta - 4.0));
00646 
00647 
00648 
00649   /* *************************************************************
00650    * END OF COEFFICIENTS IN h_+ & h_x ****************************
00651    *************************************************************** */
00652 
00653 
00654 
00655   /* Compute frequency, phase, and amplitude factors. */
00656   fFac = 1.0 / ( 4.0*LAL_TWOPI*LAL_MTSUN_SI*mTot );
00657   dt   = -params->deltaT * eta / ( 5.0*LAL_MTSUN_SI*mTot );
00658   ASSERT( dt < 0.0, stat, GENERATEPPNINSPIRALH_ETBAD,
00659     GENERATEPPNINSPIRALH_MSGETBAD );
00660  
00661   f2aFac = LAL_PI*LAL_MTSUN_SI*mTot*fFac;
00662   ASSERT( params->d != 0.0, stat, GENERATEPPNINSPIRALH_EDBAD,
00663     GENERATEPPNINSPIRALH_MSGEDBAD );
00664 
00665   /* Compute PN expansion coefficients. 
00666      - Correction to the c5 term berlow in accordance with 
00667        erratum BFIJ, PRD 71 129902  
00668      - c6 does not include log at this stage but is added later 
00669      - p6 = p[6] used in FREQ() macro
00670   */
00671   c0 = c[0] =  p[0];
00672   c1 = c[1] =  p[1];
00673   c2 = c[2] =  p[2]*( 743.0/2688.0 + eta*11.0/32.0 );
00674   c3 = c[3] = -p[3]*( 3.0*LAL_PI/10.0 );
00675   c4 = c[4] =  p[4]*( 
00676                     1855099.0/14450688.0 + eta*56975.0/258048.0 +
00677                     eta*eta*371.0/2048.0 
00678                     );
00679   c5 = c[5] =  p[5]*( -7729.0/21504.0 + eta*13.0/256.0 )*LAL_PI;
00680   c6 = c[6] = -p[6]*(
00681                     720817631400877.0/288412611379200.0 
00682                     - 107.0*LAL_GAMMA/280.0 - LAL_PI*LAL_PI*53.0/200.0
00683                     + eta*(
00684                           - 25302017977.0/4161798144.0 
00685                           + LAL_PI*LAL_PI*451.0/2048.0
00686                           )
00687                     + eta*eta*30913.0/1835008.0
00688                     + eta*eta*eta*235925.0/1769472.0
00689                     );     
00690   c7 = c[7] = -p[7]*LAL_PI*(
00691                            377033378.0/867041280.0 + eta*977650.0/2580480.0 
00692                            - eta*eta*283538.0/2580480.0
00693                            );
00694   p6 = p[6];
00695  
00696   /* Compute expansion coefficients for series in phi and dy/dx. */
00697   d0 =  c0;
00698   d1 =  c1*5.0/4.0;
00699   d2 =  c2*5.0/3.0;
00700   d3 =  c3*5.0/2.0;
00701   d4 =  c4*5.0;
00702   d5 =  c5*5.0/8.0;
00703   d6 =  p6*(
00704            831032450749357.0/57682522275840.0 - LAL_PI*LAL_PI*53.0/40.0 
00705            - 107.0*LAL_GAMMA/56.0 
00706            + eta*(
00707                  -123292747421.0/4161798144.0 + LAL_PI*LAL_PI*2255.0/2048.0 
00708                  + 385.0/48.0*(-1987.0/3080) -55.0/16.0*(-11831.0/9240.0)
00709                  )
00710            + eta*eta*(154565.0/1835008.0 - eta*1179625.0/1769472.0)
00711            );
00712   d7 = -c7*5.0/2.0;
00713   e0 =  c0*3.0;
00714   e1 =  c1*4.0;
00715   e2 =  c2*5.0;
00716   e3 =  c3*6.0;
00717   e4 =  c4*7.0;
00718   e5 =  c5*8.0;
00719   e6 =  c6*9.0;
00720   e7 =  c7*10.0;
00721 
00722   /* Use Boolean variables to exclude terms that are zero. */
00723   b0 = b[0] = ( c0 == 0.0 ? 0 : 1 );
00724   b1 = b[1] = ( c1 == 0.0 ? 0 : 1 );
00725   b2 = b[2] = ( c2 == 0.0 ? 0 : 1 );
00726   b3 = b[3] = ( c3 == 0.0 ? 0 : 1 );
00727   b4 = b[4] = ( c4 == 0.0 ? 0 : 1 );
00728   b5 = b[5] = ( c5 == 0.0 ? 0 : 1 );
00729   b6 = b[6] = ( c6 == 0.0 ? 0 : 1 );
00730   b7 = b[7] = ( c7 == 0.0 ? 0 : 1 );
00731 
00732   /* Find the leading-order frequency term. */
00733   for ( j = 0; ( j < MAXORDER ) && ( b[j] == 0 ); j++ )
00734     ;
00735   if ( j == MAXORDER ) {
00736     ABORT( stat, GENERATEPPNINSPIRALH_EPBAD,
00737      GENERATEPPNINSPIRALH_MSGEPBAD );
00738   }
00739 
00740 
00741   /*******************************************************************
00742    * COMPUTE START TIME                                              *
00743    *******************************************************************/
00744 
00745   /* First, find the normalized start frequency, and the best guess as
00746      to the start times from each term.  We require the
00747      frequency to be increasing. */
00748   yStart =  2.0/(REAL4)(Harmonics)*params->fStartIn / fFac;
00749   
00750   if ( params->fStopIn == 0.0 )
00751     yMax = 1.0/(LAL_PI*pow(6.0, 1.5)*mTot*LAL_MTSUN_SI) / fFac;
00752   else if( params->fStopIn < 0.0 )
00753     yMax = -1.0*params->fStopIn / fFac;
00754   else {
00755     ASSERT( fabs( params->fStopIn ) > params->fStartIn, stat,
00756       GENERATEPPNINSPIRALH_EFBAD, GENERATEPPNINSPIRALH_MSGEFBAD );
00757     yMax = fabs( params->fStopIn ) / fFac;
00758   }
00759     
00760   if ( ( c[j]*fFac < 0.0 ) || ( yStart < 0.0 ) || ( yMax < 0.0 ) ) {
00761      ABORT( stat, GENERATEPPNINSPIRALH_EPBAD,
00762      GENERATEPPNINSPIRALH_MSGEPBAD );
00763   }
00764 
00765   xStart = pow( yStart/c[j], 1.0/( j + 3.0 ) );
00766   xMax = LAL_SQRT2;
00767     
00768   /* The above is exact if the leading-order term is the only one in
00769    the expansion.  Check to see if there are any other terms. */
00770   for ( i = j + 1; ( i < MAXORDER ) && ( b[i] == 0 ); i++ )
00771   ;
00772   if ( i < MAXORDER )  
00773   {
00774     /* There are other terms, so we have to use bisection to find the
00775        start time. */
00776     REAL4 xLow, xHigh; /* ultimately these will bracket xStart */
00777     REAL4 yLow, yHigh; /* the normalized frequency at these times */
00778 
00779     if ( xStart > 0.39*xMax )
00780       xStart = 0.39*xMax;
00781 
00782     /* If we are ignoring PN breakdown, adjust xMax (so that it won't
00783        interfere with the start time search) and tStop. */
00784     if ( params->fStopIn < 0.0 ) 
00785     {
00786       xMax = LAL_REAL4_MAX;
00787       tStop = 0.0;
00788     } 
00789 
00790     /* If our frequency is too high, step backwards and/or forwards
00791        until we have bracketed the correct frequency. */
00792     xLow = xHigh = xStart;
00793     FREQ( yHigh, xStart);
00794     yLow = yHigh;
00795     while ( yLow > yStart ) 
00796     {
00797       xHigh = xLow;
00798       yHigh = yLow;
00799       xLow *= 0.95;
00800       FREQ( yLow, xLow );
00801     }
00802   
00803     while ( yHigh < yStart ) 
00804     {
00805       xLow = xHigh;
00806       yLow = yHigh;
00807       xHigh *= 1.05;
00808       FREQ( yHigh, xHigh );
00809       /* Check for PN breakdown. */
00810       if ( ( yHigh < yLow ) || ( xHigh > xMax ) ) 
00811       {
00812         ABORT( stat, GENERATEPPNINSPIRALH_EFBAD,
00813          GENERATEPPNINSPIRALH_MSGEFBAD );
00814       }
00815     }
00816 
00817     /* We may have gotten lucky and nailed the frequency right on.
00818          Otherwise, find xStart by root bisection. */
00819     if ( yLow == yStart )
00820       xStart = xLow;
00821     else if ( yHigh == yStart )
00822       xStart = xHigh;
00823     else 
00824     {
00825       SFindRootIn in;
00826       FreqDiffParamStruc par;
00827       in.xmax = xHigh;
00828       in.xmin = xLow;
00829       in.xacc = ACCURACY;
00830       in.function = FreqDiff;
00831       par.c = c;
00832       par.b = b;
00833       par.p6 = p6;
00834       par.y0 = yStart;
00835 
00836       TRY( LALSBisectionFindRoot( stat->statusPtr, &(xStart), &in,
00837           (void *)( &par ) ), stat );
00838     }
00839  
00840   }
00841 
00842   /* If we are ignoring PN breakdown, adjust xMax and tStop, if they
00843      haven't been adjusted already. */
00844   else if ( params->fStopIn < 0.0 ) 
00845   {
00846     xMax = LAL_REAL4_MAX;
00847     tStop = 0.0;
00848   }
00849      
00850   /* Compute initial dimensionless time, record actual initial
00851      frequency (in case it is different), and record dimensional
00852      time-to-coalescence. */
00853   t0 = pow( xStart, -8.0 );
00854   FREQ( yStart, xStart );
00855   if ( yStart >= yMax ) 
00856   {
00857     ABORT( stat, GENERATEPPNINSPIRALH_EFBAD,
00858                    GENERATEPPNINSPIRALH_MSGEFBAD );
00859   }  
00860   params->fStart = yStart*fFac;
00861   params->tc = t0 * ( 5.0*LAL_MTSUN_SI*mTot ) / eta;
00862        
00863  
00864 
00865   /*******************************************************************
00866    * GENERATE WAVEFORM                                               *
00867    *******************************************************************/
00868 
00869   /* Set up data pointers and storage. */
00870   here = head = (PPNInspiralBuffer *)
00871   LALMalloc( sizeof(PPNInspiralBuffer) );
00872   if ( !here ) 
00873   {
00874     ABORT( stat, GENERATEPPNINSPIRALH_EMEM,
00875      GENERATEPPNINSPIRALH_MSGEMEM );
00876   }
00877   here->next = NULL;
00878   h = here->h;
00879   a = here->a;
00880   f = here->f;
00881   
00882   phi = here->phi;
00883   nMax = (UINT4)( -1 );
00884   if ( params->lengthIn > 0 )
00885     nMax = params->lengthIn;
00886   nNext = BUFFSIZE;
00887   if ( nNext > nMax )
00888     nNext = nMax;
00889 
00890 
00891   /* Start integrating!  Inner loop exits each time a new buffer is
00892      required.  Outer loop has no explicit test; when a termination
00893      condition is met, we jump directly from the inner loop using a
00894      goto statement.  All goto statements jump to the terminate: label
00895      at the end of the outer loop. */
00896   n = 0;
00897   t = t0;
00898   dyMax = 0.0;
00899   y = yOld = 0.0;
00900   x = xStart;
00901   
00902   while ( 1 ) 
00903   {
00904     while ( n < nNext ) 
00905     {
00906       REAL4 f2a; /* value inside 2/3 power in amplitude functions */
00907       REAL4 phase = 0.0; /* wave phase excluding overall constants */
00908       REAL4 dydx2 = 0.0; /* dy/dx divided by x^2 */
00909 
00910       /* Check if we're still in a valid PN regime. */
00911       if ( x > xMax ) 
00912       {
00913         params->termCode = GENERATEPPNINSPIRALH_EPNFAIL;
00914         params->termDescription = GENERATEPPNINSPIRALH_MSGEPNFAIL;
00915         goto terminate;
00916       }
00917 
00918      /* Compute the normalized frequency.  This also computes the
00919         variables x2, x3, x4, x5, x6 and x7 which are used later. */
00920      FREQ( y, x );
00921 
00922      if ( y > yMax ) 
00923      {
00924        params->termCode = GENERATEPPNINSPIRALH_EFSTOP;
00925        params->termDescription = GENERATEPPNINSPIRALH_MSGEFSTOP;
00926        goto terminate;
00927      }
00928 
00929 
00930      /* Check that frequency is still increasing. */
00931      if ( b0 )
00932        dydx2 += e0;
00933      if ( b1 )
00934        dydx2 += e1*x;
00935      if ( b2 )
00936        dydx2 += e2*x2;
00937      if ( b3 )
00938        dydx2 += e3*x3;
00939      if ( b4 )
00940        dydx2 += e4*x4;
00941      if ( b5 )
00942        dydx2 += e5*x5;
00943      if ( b6 )
00944        dydx2 += (e6 + (856.0/2240.0*(2.0 + 9.0*log(2.0*x))))*x6;
00945      if ( b7 )
00946        dydx2 += e7*x7;  
00947      
00948     if ( dydx2 < 0.0 ) 
00949     {
00950        params->termCode = GENERATEPPNINSPIRALH_EFNOTMON;
00951        params->termDescription = GENERATEPPNINSPIRALH_MSGEFNOTMON;
00952        goto terminate;
00953     }
00954 
00955      if ( y - yOld > dyMax )
00956        dyMax = y - yOld;
00957      *(f++) = fFac*y;
00958 
00959      /* Compute the phase. */
00960      if ( b0 )
00961        phase += d0;
00962      if ( b1 )
00963        phase += d1*x;
00964      if ( b2 )
00965        phase += d2*x2;
00966      if ( b3 )
00967        phase += d3*x3;
00968      if ( b4 )
00969        phase += d4*x4;
00970      if ( b5 )
00971        phase += d5*log(t)*x5;
00972      if ( b6 )
00973        phase += (d6 - 8.0*107.0*log(2.0*x)/448.0)*x6;      
00974      if ( b7 )
00975        phase += d7*x7;      
00976      /* etaInv absorbs the factor of 2! */
00977      phase *= t*x3*etaInv;
00978      *(phi++) = phiC - phase;
00979            
00980      /* Compute hplus and hcross */
00981      f2a = pow(f2aFac*y, TWOTHIRDS);
00982 
00983      /* powers of frequency */
00984      fthree = pow(f2a, 1.5);
00985      ffour  = pow(f2a, 2.0);
00986      ffive  = pow(f2a, 2.5);
00987      fsix   = pow(f2a, 3.0);
00988      fseven = pow(f2a, 3.5);
00989 
00990      /* PLUS */
00991      a1 =   q[1]*a1Pthree*fthree + q[3]*a1Pfive*ffive + q[4]*a1Psix*fsix 
00992           + q[5]*a1Pseven*fseven;  
00993      a1mix = q[4]*a1Pmixsix*fsix;
00994        
00995      a2 =   q[0]*a2Ptwo*f2a + q[2]*a2Pfour*ffour + q[3]*a2Pfive*ffive 
00996           + q[4]*a2Psix*fsix + q[5]*a2Pseven*fseven;     
00997      a2mix = q[5]*a2Pmixseven*fseven;               
00998        
00999      a3 =   q[1]*a3Pthree*fthree + q[3]*a3Pfive*ffive + q[4]*a3Psix*fsix 
01000           + q[5]*a3Pseven*fseven;    
01001      a3mix = q[4]*a3Pmixsix*fsix;         
01002      
01003      a4 = q[2]*a4Pfour*ffour + q[4]*a4Psix*fsix + q[5]*a4Pseven*fseven; 
01004      a4mix = q[5]*a4Pmixseven*fseven;      
01005      
01006      a5 = q[3]*a5Pfive*ffive + q[5]*a5Pseven*fseven;      
01007      
01008      a6 = q[4]*a6Psix*fsix;
01009      
01010      a7 = q[5]*a7Pseven*fseven;
01011      
01012      /* Store first three harmonics for filtering */
01013      *(a++) = a1Pthree*fthree;
01014      *(a++) = a2Ptwo*f2a;
01015      *(a++) = a3Pthree*fthree;
01016 
01017      *(h++) = preFac*(
01018                      a1*cos(1.0/2.0*(phiC - phase)) 
01019                      + a2*cos(2.0/2.0*(phiC - phase)) 
01020                      + a3*cos(3.0/2.0*(phiC - phase)) 
01021                      + a4*cos(4.0/2.0*(phiC - phase)) 
01022                      + a5*cos(5.0/2.0*(phiC - phase)) 
01023                      + a6*cos(6.0/2.0*(phiC - phase))
01024                      + a7*cos(7.0/2.0*(phiC - phase))
01025                      + a1mix*sin(1.0/2.0*(phiC - phase)) 
01026                      + a2mix*sin(2.0/2.0*(phiC - phase)) 
01027                      + a3mix*sin(3.0/2.0*(phiC - phase))
01028                      + a4mix*sin(4.0/2.0*(phiC - phase))
01029                      );      
01030 
01031 
01032      /* CROSS */
01033      a1 =   q[1]*a1Cthree*fthree + q[3]*a1Cfive*ffive + q[4]*a1Csix*fsix 
01034           + q[5]*a1Cseven*fseven;                
01035      a1mix = q[4]*a1Cmixsix*fsix;    
01036        
01037      a2 =   q[0]*a2Ctwo*f2a + q[2]*a2Cfour*ffour + q[3]*a2Cfive*ffive 
01038           + q[4]*a2Csix*fsix + q[5]*a2Cseven*fseven;    
01039      a2mix = q[5]*a2Cmixseven*fseven;
01040      
01041      a3 =   q[1]*a3Cthree*fthree + q[3]*a3Cfive*ffive + q[4]*a3Csix*fsix 
01042           + q[5]*a3Cseven*fseven;   
01043      a3mix = q[4]*a3Cmixsix*fsix;
01044      
01045      a4 = q[2]*a4Cfour*ffour + q[4]*a4Csix*fsix + q[5]*a4Cseven*fseven;
01046      a4mix = q[5]*a4Cmixseven*fseven;  
01047      
01048      a5 = q[3]*a5Cfive*ffive + q[5]*a5Cseven*fseven;
01049      
01050      a6 = q[4]*a6Csix*fsix;
01051      
01052      a7 = q[5]*a7Cseven*fseven;
01053      
01054      *(h++) = preFac*(
01055                      a1*sin(1.0/2.0*(phiC - phase)) 
01056                      + a2*sin(2.0/2.0*(phiC - phase)) 
01057                      + a3*sin(3.0/2.0*(phiC - phase)) 
01058                      + a4*sin(4.0/2.0*(phiC - phase)) 
01059                      + a5*sin(5.0/2.0*(phiC - phase)) 
01060                      + a6*sin(6.0/2.0*(phiC - phase))
01061                      + a7*sin(7.0/2.0*(phiC - phase))
01062                      + a1mix*cos(1.0/2.0*(phiC - phase)) 
01063                      + a2mix*cos(2.0/2.0*(phiC - phase)) 
01064                      + a3mix*cos(3.0/2.0*(phiC - phase))
01065                      + a4mix*cos(4.0/2.0*(phiC - phase))
01066                      );          
01067      
01068      n++;
01069      t = t0 + n*dt;
01070      yOld = y;
01071      if ( t <= tStop ) {
01072        params->termCode = GENERATEPPNINSPIRALH_ERTOOSMALL;
01073        params->termDescription = GENERATEPPNINSPIRALH_MSGERTOOSMALL;
01074        goto terminate;
01075      }
01076      x = pow( t, -0.125 );
01077 
01078    }
01079 
01080    /* We've either filled the buffer or we've exceeded the maximum
01081       length.  If the latter, we're done! */
01082    if ( n >= nMax ) {
01083      params->termCode = GENERATEPPNINSPIRALH_ELENGTH;
01084      params->termDescription = GENERATEPPNINSPIRALH_MSGELENGTH;
01085      
01086    }
01087 
01088 
01089    /* Otherwise, allocate the next buffer. */
01090    here->next =
01091       (PPNInspiralBuffer *)LALMalloc( sizeof(PPNInspiralBuffer) );
01092    here = here->next;
01093    if ( !here ) {
01094      FREELIST( head );
01095      ABORT( stat, GENERATEPPNINSPIRALH_EMEM,
01096        GENERATEPPNINSPIRALH_MSGEMEM );
01097    }
01098   
01099    here->next = NULL;
01100    h = here->h;
01101    a = here->a;
01102    f = here->f;
01103    phi = here->phi;
01104    nNext += BUFFSIZE;
01105    if ( nNext > nMax )
01106      nNext = nMax;
01107 
01108 
01109  }
01110 
01111 
01112   /*******************************************************************
01113    * CLEANUP                                                         *
01114    *******************************************************************/
01115 
01116   /* The above loop only exits by triggering one of the termination
01117      conditions, which jumps to the following point for cleanup and
01118      return. */
01119  terminate:
01120 
01121   /* First, set remaining output parameter fields. */
01122   params->dfdt = dyMax*fFac*params->deltaT;
01123   params->fStop = yOld*fFac;
01124   params->length = n;
01125 
01126   /* Allocate the output structures. */
01127   if ( ( output->h = (REAL4TimeVectorSeries *)
01128        LALMalloc( sizeof(REAL4TimeVectorSeries) ) ) == NULL ) 
01129   {
01130     FREELIST( head );
01131     ABORT( stat, GENERATEPPNINSPIRALH_EMEM,
01132      GENERATEPPNINSPIRALH_MSGEMEM );
01133   }
01134   memset( output->h, 0, sizeof(REAL4TimeVectorSeries) );
01135   
01136   if ( ( output->a = (REAL4TimeVectorSeries *)
01137        LALMalloc( sizeof(REAL4TimeVectorSeries) ) ) == NULL ) 
01138   {
01139     FREELIST( head );
01