summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorLászló Veréb <lipiboy@gmail.com>2010-08-04 14:30:12 (GMT)
committerLászló Veréb <lipiboy@gmail.com>2010-08-04 15:34:06 (GMT)
commitd8c46ea014c6ea8dc22060be90c2dabef92ce02e (patch)
tree7166f2758ea1fb23d52163add89283b3f12c922a
parent9abf87262a96f19dd5563f271bd1a8050412d761 (diff)
The SpinQuadTaylor waveform generator code is integrated.
-rw-r--r--lalinspiral/src/GenerateInspiral.c43
-rw-r--r--lalinspiral/src/GenerateInspiral.h10
-rw-r--r--lalinspiral/src/LALInspiral.h17
-rw-r--r--lalinspiral/src/LALInspiralChooseModel.c11
-rw-r--r--lalinspiral/src/LALInspiralWave.c14
-rw-r--r--lalinspiral/src/LALSQTPNIntegrator.c72
-rw-r--r--lalinspiral/src/LALSQTPNIntegrator.h63
-rw-r--r--lalinspiral/src/LALSQTPNWaveform.c372
-rw-r--r--lalinspiral/src/LALSQTPNWaveform.h225
-rw-r--r--lalinspiral/src/LALSQTPNWaveformInterface.c296
-rw-r--r--lalinspiral/src/LALSQTPNWaveformInterface.h104
-rw-r--r--lalinspiral/src/Makefile.am6
-rw-r--r--lalinspiral/test/LALSQTPNWaveformTest.c139
-rw-r--r--lalmetaio/src/LIGOMetadataTables.h1
14 files changed, 1366 insertions, 7 deletions
diff --git a/lalinspiral/src/GenerateInspiral.c b/lalinspiral/src/GenerateInspiral.c
index f59118e..328ee98 100644
--- a/lalinspiral/src/GenerateInspiral.c
+++ b/lalinspiral/src/GenerateInspiral.c
@@ -1,5 +1,5 @@
/*
-* Copyright (C) 2007 Stas Babak, Drew Keppel, Duncan Brown, Eirini Messaritaki, Gareth Jones, Thomas Cokelaer
+* Copyright (C) 2007 Stas Babak, Drew Keppel, Duncan Brown, Eirini Messaritaki, Gareth Jones, Thomas Cokelaer, Laszlo Vereb
*
* 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
@@ -49,7 +49,7 @@ $Id$
\begin{description}
\item[\texttt{LALGenerateInspiral()}] create an inspiral binary
waveform generated either by the \texttt{inspiral} package (EOB,
-EOBNR, PadeT1, TaylorT1, TaylorT2, TaylorT3, SpinTaylor) or the
+EOBNR, PadeT1, TaylorT1, TaylorT2, TaylorT3, SpinTaylor, SpinQuadTaylor) or the
\texttt{inject} package (GeneratePPN). It is used in the module
\texttt{FindChirpSimulation} in \texttt{findchirp} package.
@@ -188,6 +188,12 @@ LALGenerateInspiral(
{
inspiralParams.approximant = approximant;
inspiralParams.order = order;
+ if (approximant == SpinQuadTaylor) {
+ xlalErrno = 0;
+ if (XLALGetSpinInteractionFromString(&inspiralParams.spinInteraction, thisEvent->waveform) == XLAL_FAILURE) {
+ ABORTXLAL(status);
+ }
+ }
/* We fill ppnParams */
LALGenerateInspiralPopulatePPN(status->statusPtr, ppnParams, thisEvent);
@@ -332,6 +338,33 @@ LALGetOrderFromString(
RETURN( status );
}
+int XLALGetSpinInteractionFromString(LALSpinInteraction *inter, CHAR *thisEvent) {
+ static const char *func = "XLALGetSpinInteractionFromString";
+
+ if (strstr(thisEvent, "ALL")) {
+ *inter = LAL_AllInter;
+ } else if (strstr(thisEvent, "NO")) {
+ *inter = LAL_NOInter;
+ } else {
+ *inter = LAL_SOInter;
+ if (strstr(thisEvent, "SO")) {
+ *inter |= LAL_SOInter;
+ }
+ if (strstr(thisEvent, "QM")) {
+ *inter |= LAL_QMInter;
+ }
+ if (strstr(thisEvent, "SELF")) {
+ *inter |= LAL_SSselfInter;
+ }
+ if (strstr(thisEvent, "SS")) {
+ *inter |= LAL_SSInter;
+ }
+ if (*inter == LAL_NOInter) {
+ XLAL_ERROR(func, XLAL_EDOM);
+ }
+ }
+ return XLAL_SUCCESS;
+}
/* <lalVerbatim file="LALGetApproxFromStringCP"> */
void
@@ -376,6 +409,10 @@ LALGetApproximantFromString(
{
*approximant = SpinTaylor;
}
+ else if ( strstr(thisEvent, "SpinQuadTaylor" ) )
+ {
+ *approximant = SpinQuadTaylor;
+ }
else if ( strstr(thisEvent, "PadeT1" ) )
{
*approximant = PadeT1;
@@ -537,6 +574,8 @@ LALGenerateInspiralPopulateInspiral(
inspiralParams->orbitTheta0 = thisEvent->theta0;
inspiralParams->orbitPhi0 = thisEvent->phi0;
+ inspiralParams->qmParameter[0] = thisEvent->qmParameter[0];
+ inspiralParams->qmParameter[1] = thisEvent->qmParameter[1];
DETATCHSTATUSPTR( status );
RETURN( status );
diff --git a/lalinspiral/src/GenerateInspiral.h b/lalinspiral/src/GenerateInspiral.h
index 9217941..cf5306d 100644
--- a/lalinspiral/src/GenerateInspiral.h
+++ b/lalinspiral/src/GenerateInspiral.h
@@ -1,5 +1,5 @@
/*
-* Copyright (C) 2007 Drew Keppel, Duncan Brown, Gareth Jones, Peter Shawhan, Thomas Cokelaer
+* Copyright (C) 2007 Drew Keppel, Duncan Brown, Gareth Jones, Peter Shawhan, Thomas Cokelaer, Laszlo Vereb
*
* 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
@@ -131,6 +131,14 @@ LALGetApproximantFromString(
Approximant *result
);
+/** Convert a string provided by the #CoherentGW structure in order to retrieve
+ * the approximant of the waveform to generate.
+ * @param[out] inter : the level of the spin interaction
+ * @param[in] thisEvent : string containing the spin interaction
+ * @return error code
+ */
+int XLALGetSpinInteractionFromString(LALSpinInteraction *inter, CHAR *thisEvent);
+
/* three function to populate the needed structures */
void
LALGenerateInspiralPopulatePPN(
diff --git a/lalinspiral/src/LALInspiral.h b/lalinspiral/src/LALInspiral.h
index bd7afd8..56db094 100644
--- a/lalinspiral/src/LALInspiral.h
+++ b/lalinspiral/src/LALInspiral.h
@@ -1,5 +1,5 @@
/*
-* Copyright (C) 2007 Stas Babak, David Churches, Drew Keppel, Duncan Brown, Jolien Creighton, David McKechan, Patrick Brady, Peter Shawhan, Reinhard Prix, B.S. Sathyaprakash, Anand Sengupta, Craig Robinson , Sean Seader, Thomas Cokelaer
+* Copyright (C) 2007 Stas Babak, David Churches, Drew Keppel, Duncan Brown, Jolien Creighton, David McKechan, Patrick Brady, Peter Shawhan, Reinhard Prix, B.S. Sathyaprakash, Anand Sengupta, Craig Robinson , Sean Seader, Thomas Cokelaer, Laszlo Vereb
*
* 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
@@ -185,6 +185,7 @@ by \texttt{TaylorT1} approximant (see Ref. \cite{dis2} for details). Outputs a f
Vallisneri including spin effects\cite{BCV03b}. Outputs a frequency-domain wave.
\item \texttt{SpinTaylorT3} Spinning case T3 models
\item \texttt{SpinTaylor} Spinning case PN models (should replace SpinTaylorT3 in the future)
+\item \texttt{SpinQuadTaylor} Spinning case PN models with quadrupole-monopole and self-spin interaction.
\item \texttt{FindChirpSP} The stationary phase templates implemented by FindChirpSPTemplate in the findchirp package (equivalent to TaylorF2 at twoPN order).
\item \texttt{GeneratePPN} The time domain templates generated by LALGeneratePPNInspiral() in the inject package (equivalent to TaylorT3 at twoPN order).
\item \texttt{FrameFile} The waveform contains arbitrary data read from a frame file.
@@ -455,6 +456,7 @@ typedef enum {
BCVSpin,
SpinTaylorT3,
SpinTaylor,
+ SpinQuadTaylor,
FindChirpSP,
FindChirpPTF,
GeneratePPN,
@@ -498,6 +500,17 @@ typedef enum {
/* </lalVerbatim> */
+/** Enumeration to specify which component will be used in the waveform
+ * generation. Their combination also can be used by the bitwise or.
+ **/
+typedef enum {
+ LAL_NOInter = 0,
+ LAL_SOInter = 1, ///< Spin-orbit interaction
+ LAL_SSInter = LAL_SOInter << 1, ///< Spin-spin interaction
+ LAL_SSselfInter = LAL_SSInter << 1, ///< Spin-spin-self interaction
+ LAL_QMInter = LAL_SSselfInter << 1, ///< quadrupole-monopole interaction
+ LAL_AllInter = LAL_QMInter << 1 ///< all interactions
+} LALSpinInteraction;
/* <lalLaTeX>
\idx[Type]{InputMasses}
@@ -610,6 +623,8 @@ tagInspiralTemplate
* For spinBCV searches, (in 4 dimensions) Gamma[0,...,9] would be required.
*/
REAL4 Gamma[10];
+ REAL4 qmParameter[2];
+ LALSpinInteraction spinInteraction;
struct tagInspiralTemplate *next;
struct tagInspiralTemplate *fine;
diff --git a/lalinspiral/src/LALInspiralChooseModel.c b/lalinspiral/src/LALInspiralChooseModel.c
index fadb104..b55274b 100644
--- a/lalinspiral/src/LALInspiralChooseModel.c
+++ b/lalinspiral/src/LALInspiralChooseModel.c
@@ -1,5 +1,5 @@
/*
-* Copyright (C) 2007 David Churches, Jolien Creighton, David McKechan, B.S. Sathyaprakash, Thomas Cokelaer, Duncan Brown
+* Copyright (C) 2007 David Churches, Jolien Creighton, David McKechan, B.S. Sathyaprakash, Thomas Cokelaer, Duncan Brown, Laszlo Vereb
*
* 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
@@ -508,6 +508,7 @@ LALInspiralChooseModel(
case TaylorF2:
case SpinTaylorT3:
case SpinTaylor:
+ case SpinQuadTaylor:
case IMRPhenomA:
case IMRPhenomB:
case IMRPhenomFA:
@@ -550,6 +551,7 @@ LALInspiralChooseModel(
case TaylorF2:
case SpinTaylorT3:
case SpinTaylor:
+ case SpinQuadTaylor:
case IMRPhenomA:
case IMRPhenomB:
case IMRPhenomFA:
@@ -590,6 +592,7 @@ LALInspiralChooseModel(
case TaylorF2:
case SpinTaylorT3:
case SpinTaylor:
+ case SpinQuadTaylor:
case IMRPhenomA:
case IMRPhenomB:
case IMRPhenomFA:
@@ -633,6 +636,7 @@ LALInspiralChooseModel(
case TaylorF2:
case SpinTaylorT3:
case SpinTaylor:
+ case SpinQuadTaylor:
/*
The value vlsoT4 is too large and doesn't work sometimes;
so we use vlsoT2.
@@ -680,6 +684,7 @@ LALInspiralChooseModel(
case TaylorF2:
case SpinTaylorT3:
case SpinTaylor:
+ case SpinQuadTaylor:
/*
The value vlsoT4 is too large and doesn't work with 2.5 PN
Taylor approximant; so we use vlsoT2.
@@ -728,6 +733,7 @@ LALInspiralChooseModel(
case TaylorF2:
case SpinTaylorT3:
case SpinTaylor:
+ case SpinQuadTaylor:
/*
vlsoT6 is as yet undetermined and vlsoT4 is too large in
certain cases (TaylorT2 crashes for (1.4,10)); using vlsoT2;
@@ -776,6 +782,7 @@ LALInspiralChooseModel(
case TaylorF2:
case SpinTaylorT3:
case SpinTaylor:
+ case SpinQuadTaylor:
ak->vn = ak->vlso = vlso = ak->vlsoT2;
f->dEnergy = dEt6;
f->flux = Ft7;
@@ -829,6 +836,7 @@ LALInspiralChooseModel(
case TaylorF2:
case SpinTaylorT3:
case SpinTaylor:
+ case SpinQuadTaylor:
case PadeT1:
case PadeF1:
case TaylorEt:
@@ -857,6 +865,7 @@ LALInspiralChooseModel(
case TaylorF2:
case SpinTaylorT3:
case SpinTaylor:
+ case SpinQuadTaylor:
case TaylorEt:
case TaylorT4:
case TaylorN:
diff --git a/lalinspiral/src/LALInspiralWave.c b/lalinspiral/src/LALInspiralWave.c
index a56d4f7..7687c62 100644
--- a/lalinspiral/src/LALInspiralWave.c
+++ b/lalinspiral/src/LALInspiralWave.c
@@ -1,5 +1,5 @@
/*
-* Copyright (C) 2007 David Churches, Duncan Brown, Jolien Creighton, David McKechan, B.S. Sathyaprakash, Thomas Cokelaer
+* Copyright (C) 2007 David Churches, Duncan Brown, Jolien Creighton, David McKechan, B.S. Sathyaprakash, Thomas Cokelaer, Laszlo Vereb
*
* 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
@@ -106,7 +106,7 @@ Depending on the user inputs one of the following functions is called:\\
\begin{itemize}
\item A time-domain waveform is returned when the {\tt approximant} is one of
- {\tt TaylorT1, TaylorT2, TaylorT3, PadeT1, EOB, SpinTaylorT3}
+ {\tt TaylorT1, TaylorT2, TaylorT3, PadeT1, EOB, SpinTaylorT3, SpinQuadTaylor}
\item A frequency-domain waveform is returned when the {\tt approximant} is one of
{\tt TaylorF1, TaylorF2, BCV}.
In these cases the code returns the real and imagninary parts of the
@@ -128,6 +128,7 @@ Depending on the user inputs one of the following functions is called:\\
#include <lal/LALNoiseModels.h>
#include <lal/LALStdlib.h>
#include <lal/GeneratePPNInspiral.h>
+#include <lal/LALSQTPNWaveformInterface.h>
NRCSID (LALINSPIRALWAVEC, "$Id$");
@@ -216,6 +217,9 @@ LALInspiralWave(
LALSTPNWaveform(status->statusPtr, signalvec, params);
CHECKSTATUSPTR(status);
break;
+ case SpinQuadTaylor:
+ TRY(LALSQTPNWaveform(status->statusPtr, signalvec, params), status);
+ break;
case AmpCorPPN:
LALInspiralAmplitudeCorrectedWave(status->statusPtr, signalvec, params);
CHECKSTATUSPTR(status);
@@ -314,6 +318,9 @@ LALInspiralWaveTemplates(
LALSTPNWaveformTemplates(status->statusPtr, signalvec1, signalvec2, params);
CHECKSTATUSPTR(status);
break;
+ case SpinQuadTaylor:
+ TRY(LALSTPNWaveformTemplates(status->statusPtr, signalvec1, signalvec2, params), status);
+ break;
case AmpCorPPN:
LALInspiralAmplitudeCorrectedWaveTemplates(status->statusPtr, signalvec1, signalvec2, params);
CHECKSTATUSPTR(status);
@@ -389,6 +396,9 @@ LALInspiralWaveForInjection(
LALSTPNWaveformForInjection(status->statusPtr, waveform, inspiralParams, ppnParams);
CHECKSTATUSPTR(status);
break;
+ case SpinQuadTaylor:
+ TRY(LALSQTPNWaveformForInjection(status->statusPtr, waveform, inspiralParams, ppnParams), status);
+ break;
case AmpCorPPN:
LALInspiralAmplitudeCorrectedWaveForInjection(status->statusPtr, waveform, inspiralParams, ppnParams);
CHECKSTATUSPTR(status);
diff --git a/lalinspiral/src/LALSQTPNIntegrator.c b/lalinspiral/src/LALSQTPNIntegrator.c
new file mode 100644
index 0000000..33c8272
--- /dev/null
+++ b/lalinspiral/src/LALSQTPNIntegrator.c
@@ -0,0 +1,72 @@
+/**
+ * @file LALSQTPNIntegrator.c
+ * Contains the function definitions needed by the integration method.
+ * @author László Veréb
+ * @date 2010.05.21.
+ */
+
+#include <lal/LALSQTPNIntegrator.h>
+#include <lal/LALSQTPNWaveform.h>
+
+NRCSID (LALSQTPNINTEGRATORC, "$Id LALSQTPNIntegrator.c$");
+
+int XLALSQTPNIntegratorInit(LALSQTPNIntegratorSystem *integrator, INT2 num, void *params,
+ int(*derivator)(REAL8, const REAL8[], REAL8[], void *)) {
+ static const char *func = "XLALSQTPNIntegratorSystem";
+
+ // Check for input errors
+ if (num <= 0) {
+ XLAL_ERROR(func, XLAL_EBADLEN);
+ }
+ if (!params) {
+ XLAL_ERROR(func, XLAL_EFAULT);
+ }
+ if (!derivator) {
+ XLAL_ERROR(func, XLAL_EFAULT);
+ }
+
+ // Initialise GSL integrator
+ integrator->type = gsl_odeiv_step_rkf45;
+ integrator->system.jacobian = NULL;
+ integrator->system.dimension = num;
+ integrator->system.params = params;
+ integrator->system.function = derivator;
+ XLAL_CALLGSL(integrator->step = gsl_odeiv_step_alloc(integrator->type, num));
+ XLAL_CALLGSL(integrator->control = gsl_odeiv_control_standard_new(1.0e-2, 1.0e-2, 1., 1.));
+ XLAL_CALLGSL(integrator->evolve = gsl_odeiv_evolve_alloc(num));
+
+ // Check if the integrator is correctly allocated
+ if (!(integrator->step) || !(integrator->control) || !(integrator->evolve)) {
+ XLALSQTPNIntegratorFree(integrator);
+ XLAL_ERROR(func, XLAL_ENOMEM);
+ }
+ return XLAL_SUCCESS;
+}
+
+void XLALSQTPNIntegratorFree(LALSQTPNIntegratorSystem *integrator) {
+ if (integrator->evolve) {
+ XLAL_CALLGSL(gsl_odeiv_evolve_free(integrator->evolve));
+ }
+ if (integrator->control) {
+ XLAL_CALLGSL(gsl_odeiv_control_free(integrator->control));
+ }
+ if (integrator->step) {
+ XLAL_CALLGSL(gsl_odeiv_step_free(integrator->step));
+ }
+}
+
+int XLALSQTPNIntegratorFunc(REAL8 values[], LALSQTPNIntegratorSystem *integrator, REAL8 step) {
+ static const char *func = "XLALSQTPNIntegratorFunc";
+ REAL8 time = 0., time_Old, step_X = step;
+ while (time < step) {
+ time_Old = time;
+ XLAL_CALLGSL(gsl_odeiv_evolve_apply(integrator->evolve,
+ integrator->control, integrator->step,
+ &(integrator->system), &time, step, &step_X, values));
+ if (time == time_Old) {
+ memset(values, 0, integrator->system.dimension * sizeof(REAL8));
+ XLAL_ERROR(func, XLAL_EFUNC);
+ }
+ }
+ return XLAL_SUCCESS;
+}
diff --git a/lalinspiral/src/LALSQTPNIntegrator.h b/lalinspiral/src/LALSQTPNIntegrator.h
new file mode 100644
index 0000000..14632a0
--- /dev/null
+++ b/lalinspiral/src/LALSQTPNIntegrator.h
@@ -0,0 +1,63 @@
+/**
+ * @file LALSQTPNIntegrator.h
+ * Contains the function declarations and structures needed by the
+ * integration method.
+ * @author László Veréb
+ * @date 2010.05.21.
+ */
+
+#ifndef LALSQTPNINTEGRATOR_H
+#define LALSQTPNINTEGRATOR_H
+
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_odeiv.h>
+
+#include <lal/LALInspiral.h>
+#include <lal/LALGSL.h>
+#include <lal/LALStatusMacros.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+NRCSID (LALSQTPNINTEGRATORH, "$Id LALSQTPNIntegrator.h$");
+
+/** The structure contains the integration method and its settings.
+ */
+typedef struct tagLALSQTPNIntegratorSystem{
+ const gsl_odeiv_step_type* type;
+ gsl_odeiv_step* step;
+ gsl_odeiv_control* control;
+ gsl_odeiv_evolve* evolve;
+ gsl_odeiv_system system;
+} LALSQTPNIntegratorSystem;
+
+/** The function initialize the integration method.
+ * @param[out] integrator : the structure containing the integration method
+ * @param[in] num : the number of the dynamic variables
+ * @param[in] params : the parameters used in the derivative function
+ * @param[in] derivator : pointer to the derivative function
+ */
+int XLALSQTPNIntegratorInit(LALSQTPNIntegratorSystem *integrator, INT2 num, void *params,
+ int(*derivator)(REAL8, const REAL8[], REAL8[], void *));
+
+/** The function deallocates the memory allocated for the integrator
+ * function.
+ * @param[in] integrator : the structure containing the integration method
+ */
+void XLALSQTPNIntegratorFree(LALSQTPNIntegratorSystem *integrator);
+
+/** The function evolves the system with the given time-step.
+ * @param[in,out] values : as input parameters the system's actual position,
+ * as ouput the system's next position.
+ * @param[in] integrator : the integration method
+ * @param[in] step : the step size
+ */
+int XLALSQTPNIntegratorFunc(REAL8 values[], LALSQTPNIntegratorSystem *integrator, REAL8 step);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* LALSQTPNINTEGRATOR_H */
diff --git a/lalinspiral/src/LALSQTPNWaveform.c b/lalinspiral/src/LALSQTPNWaveform.c
new file mode 100644
index 0000000..80d2d25
--- /dev/null
+++ b/lalinspiral/src/LALSQTPNWaveform.c
@@ -0,0 +1,372 @@
+/**
+ * @file LALSQTPNWaveform.c
+ * Contains the function definition to create GWforms.
+ * @author László Veréb
+ * @date 2010.05.21.
+ */
+
+#include <lal/LALSQTPNWaveform.h>
+#include <lal/LALSQTPNIntegrator.h>
+#include <lal/LALSQTPNWaveformInterface.h>
+
+NRCSID (LALSQTPNWAVEFORMC, "$Id LALSQTPN_Waveform.c$");
+
+#define UNUSED(expr) do { (void)(expr); } while (0)
+
+/** The macro function returns the square of the argument.
+ * Do not use with incrementing or decrementing operators!
+ * @param[in] a : the number
+ * @return the square of the number
+ */
+#define SQR(a) ((a)*(a))
+
+/** The macro function calculates the scalar product of two vectors.
+ * @param[in] a1 : the left vector
+ * @param[in] a2 : the right vector
+ * @return the product
+ */
+#define SCALAR_PRODUCT3(a1, a2) \
+ ((a1)[0] * (a2)[0] + (a1)[1] * (a2)[1] + (a1)[2] * (a2)[2]);
+
+/** The macro function calculates the vector product of two vectors.
+ * @param[in] left : the left vector
+ * @param[in] right : the right vector
+ * @param[out] product : the vector product
+ */
+#define VECTOR_PRODUCT3(left, right, product)\
+ (product)[0] = ((left)[1] * (right)[2] - (left)[2] * (right)[1]);\
+ (product)[1] = ((left)[2] * (right)[0] - (left)[0] * (right)[2]);\
+ (product)[2] = ((left)[0] * (right)[1] - (left)[1] * (right)[0]);
+
+void XLALSQTPNFillCoefficients(LALSQTPNWaveformParams * const params) {
+
+ // variable declaration and initialization
+ REAL8 thetahat = 1039. / 4620.;
+ REAL8 spin_MPow2[2];
+ REAL8 m_m[2] = { params->mass[1] / params->mass[0], params->mass[0]
+ / params->mass[1] };
+ REAL8 piPow2 = SQR(LAL_PI);
+ REAL8 etaPow2 = SQR(params->eta);
+ REAL8 etaPow3 = etaPow2 * params->eta;
+ INT2 i;
+ for (i = 0; i < 2; i++) {
+ spin_MPow2[i] = params->chiAmp[i] * SQR(params->mass[i])
+ /SQR(params->totalMass);
+ }
+
+ // calculating the coefficients
+ params->coeff.domegaGlobal = params->eta * 96. / 5.;
+ for (i = LAL_PNORDER_NEWTONIAN; i < LAL_PNORDER_PSEUDO_FOUR; i += 2) {
+ params->coeff.meco[i] = -0.5 * params->eta * (REAL8) (i + 2) / 3.;
+ }
+ switch (params->order) {
+ case LAL_PNORDER_THREE_POINT_FIVE:
+ params->coeff.domega[LAL_PNORDER_THREE_POINT_FIVE] = (-4415.
+ / 4032. + params->eta * 358675. / 6048. + etaPow2 * 91495.
+ / 1512.) * LAL_PI;
+ case LAL_PNORDER_THREE:
+ params->coeff.domega[LAL_PNORDER_THREE]
+ = (16447322263. / 139708800. - LAL_GAMMA * 1712. / 105.
+ + piPow2 * 16. / 3.) + (-273811877. / 1088640.
+ + piPow2 * 451. / 48. - thetahat * 88. / 3.)
+ * params->eta + etaPow2 * 541. / 896. - etaPow3
+ * 5605. / 2592.;
+ params->coeff.domegaLN = -856. / 105.;
+ params->coeff.meco[LAL_PNORDER_THREE] *= -675. / 64. + (209323.
+ / 4032. - 205. * piPow2 / 96. + (110. / 9.) * (1987.
+ / 3080.)) * params->eta - 155. * etaPow2 / 96. - 35.
+ * etaPow3 / 5184.;
+ case LAL_PNORDER_TWO_POINT_FIVE:
+ params->coeff.domega[LAL_PNORDER_TWO_POINT_FIVE] = -(4159. + 15876.
+ * params->eta) * LAL_PI / 672.;
+ case LAL_PNORDER_TWO:
+ params->coeff.domega[LAL_PNORDER_TWO] = 34103. / 18144.
+ + params->eta * 13661. / 2016. + etaPow2 * 59. / 18.;
+ params->coeff.domegaSSselfConst = 0.;
+ params->coeff.domegaQMConst = 0.;
+ if ((params->spinInteraction & LAL_SSInter) == LAL_SSInter) {
+ params->coeff.dchihSS[0] = spin_MPow2[1] / 2.;
+ params->coeff.dchihSS[1] = spin_MPow2[0] / 2.;
+ params->coeff.domegaSS[0] = 721. * params->eta
+ * params->chiAmp[0] * params->chiAmp[1] / 48.;
+ params->coeff.domegaSS[1] = -247. * params->coeff.domegaSS[0]
+ / 721.;
+ params->coeff.mecoSS = -spin_MPow2[0] * spin_MPow2[1];
+ }
+ if ((params->spinInteraction & LAL_SSselfInter) == LAL_SSselfInter) {
+ for (i = 0; i < 2; i++) {
+ params->coeff.domegaSSself[i] = -spin_MPow2[i] * params->chiAmp[i]
+ / 96.;
+ params->coeff.domegaSSselfConst -= 7. * params->coeff.domegaSSself[i];
+ }
+ }
+ if ((params->spinInteraction & LAL_QMInter) == LAL_QMInter) {
+ for (i = 0; i < 2; i++) {
+ params->coeff.domegaQM[i] = spin_MPow2[i] * params->chiAmp[i]
+ * params->qmParameter[i] * 7.5;
+ params->coeff.domegaQMConst -= params->coeff.domegaQM[i] / 3.;
+ params->coeff.dchihQM[i] = -params->qmParameter[i] * params->eta
+ * params->chiAmp[i] * 3. / 2.;
+ }
+ params->coeff.mecoQM = 2. * params->eta;
+ }
+ params->coeff.meco[LAL_PNORDER_TWO] *= (-81. + 57. * params->eta
+ - etaPow2) / 24.;
+ case LAL_PNORDER_ONE_POINT_FIVE:
+ params->coeff.domega[LAL_PNORDER_ONE_POINT_FIVE] = 4. * LAL_PI;
+ if (params->spinInteraction != 0) {
+ for (i = 0; i < 2; i++) {
+ params->coeff.dchihSO[i] = (4. + 3. * m_m[i]) * params->eta
+ / 2.;
+ params->coeff.dLNh[i] = -spin_MPow2[i] / params->eta;
+ params->coeff.domegaSO[i] = -spin_MPow2[i] * (113. + 75.
+ * m_m[i]) / 12.;
+ params->coeff.mecoSO[i] = -spin_MPow2[i] * 5. * params->eta
+ * (4. + 3. * m_m[i]) / 9.;
+ }
+ }
+ case LAL_PNORDER_ONE:
+ params->coeff.domega[LAL_PNORDER_ONE]
+ = -(743. + 924. * params->eta) / 336.;
+ params->coeff.meco[LAL_PNORDER_ONE] *= -(9. + params->eta) / 12.;
+ case LAL_PNORDER_HALF:
+ params->coeff.domega[LAL_PNORDER_HALF] = 0.;
+ case LAL_PNORDER_NEWTONIAN:
+ params->coeff.domega[LAL_PNORDER_NEWTONIAN] = 1.;
+ default:
+ break;
+ }
+}
+
+int LALSQTPNDerivator(REAL8 t, const REAL8 values[], REAL8 dvalues[], void * param) {
+
+ // variable declaration and initialization
+ LALSQTPNWaveformParams *params = param;
+ UNUSED(t);
+ const REAL8 *chi_p[2] = { values + LALSQTPN_CHIH1_1, values + LALSQTPN_CHIH2_1 };
+ UINT2 i, j, k; // indexes
+ memset(dvalues, 0, LALSQTPN_NUM_OF_VAR * sizeof(REAL8));
+ REAL8 omegaPowi_3[8];
+ omegaPowi_3[0] = 1.;
+ omegaPowi_3[1] = cbrt(values[LALSQTPN_OMEGA]);
+ for (i = 2; i < 8; i++) {
+ omegaPowi_3[i] = omegaPowi_3[i - 1] * omegaPowi_3[1];
+ }
+ REAL8 SS_Omega, SSself_Omega, QM_Omega;
+ SS_Omega = SSself_Omega = QM_Omega = 0.;
+ REAL8 chih1chih2, chih1xchih2[2][3], LNhchih[2], LNhxchih[2][3];
+ chih1chih2 = SCALAR_PRODUCT3(chi_p[0], chi_p[1]);
+ for (i = 0; i < 2; i++) {
+ LNhchih[i] = SCALAR_PRODUCT3(values + LALSQTPN_LNH_1, chi_p[i]);
+ VECTOR_PRODUCT3(values + LALSQTPN_LNH_1, chi_p[i], LNhxchih[i]);
+ }
+
+ // calculating domega and MECO without the spin components
+ for (i = LAL_PNORDER_NEWTONIAN; i <= params->order; i++) {
+ dvalues[LALSQTPN_OMEGA] += params->coeff.domega[i] * omegaPowi_3[i];
+ }
+ dvalues[LALSQTPN_MECO] += params->coeff.meco[0] / omegaPowi_3[1];
+ for (i = LAL_PNORDER_NEWTONIAN + 2; i <= params->order; i += 2) {
+ dvalues[LALSQTPN_MECO] += params->coeff.meco[i] * omegaPowi_3[i - 1];
+ }
+
+ // calculating the other derivatives and the domega and MECO with spin
+ // components
+ switch (params->order) {
+ case LAL_PNORDER_THREE_POINT_FIVE:
+ case LAL_PNORDER_THREE:
+ dvalues[LALSQTPN_OMEGA] += params->coeff.domegaLN
+ * log(16. * omegaPowi_3[2])
+ * omegaPowi_3[LAL_PNORDER_THREE];
+ case LAL_PNORDER_TWO_POINT_FIVE:
+ case LAL_PNORDER_TWO:
+ if ((params->spinInteraction & LAL_SSInter) == LAL_SSInter) {
+ // SS for domega
+ SS_Omega = params->coeff.domegaSS[0] * LNhchih[0] * LNhchih[1]
+ + params->coeff.domegaSS[1] * chih1chih2;
+ // SS for MECO
+ dvalues[LALSQTPN_MECO] += params->coeff.mecoSS * (chih1chih2 - 3
+ * LNhchih[0] * LNhchih[1]) * omegaPowi_3[3];
+ // SS for dchih
+ for (i = 0; i < 2; i++) {
+ k = (i + 1) % 2; // the opposite index
+ VECTOR_PRODUCT3(chi_p[k], chi_p[i], chih1xchih2[i]);
+ for (j = 0; j < 3; j++) {
+ // the 3*index is used, to acces the first spin, if index=0,
+ // otherwise the second spin
+ dvalues[LALSQTPN_CHIH1_1 + 3 * i + j] += params->coeff.dchihSS[i] *
+ (chih1xchih2[i][j] - 3. * LNhchih[k] * LNhxchih[i][j]) * omegaPowi_3[6];
+ }
+ }
+ }
+ if ((params->spinInteraction & LAL_SSselfInter) == LAL_SSselfInter) {
+ // SSself for domega
+ SSself_Omega = params->coeff.domegaSSselfConst;
+ for (i = 0; i < 2; i++) {
+ SSself_Omega += params->coeff.domegaSSself[i] * SQR(LNhchih[i]);
+ }
+ }
+ if ((params->spinInteraction & LAL_QMInter) == LAL_QMInter) {
+ QM_Omega = params->coeff.domegaQMConst;
+ for (i = 0; i < 2; i++) {
+ QM_Omega += params->coeff.domegaQM[i] * LNhchih[i];
+ // QM for dchih
+ for (j = 0; j < 3; j++) {
+ // the 3*index is used, to acces the first spin, if index=0,
+ // otherwise the second spin
+ dvalues[LALSQTPN_CHIH1_1 + 3 * i + j] += params->coeff.dchihQM[i]
+ * LNhchih[i] * LNhxchih[i][j] * omegaPowi_3[6];
+ }
+ }
+ // QM for MECO
+ dvalues[LALSQTPN_MECO] += params->coeff.mecoQM * QM_Omega * omegaPowi_3[3];
+ }
+ dvalues[LALSQTPN_OMEGA] += (QM_Omega + SSself_Omega + SS_Omega)
+ * omegaPowi_3[LAL_PNORDER_TWO];
+ case LAL_PNORDER_ONE_POINT_FIVE:
+ if (params->spinInteraction != 0) {
+ // SO for domega and MECO
+ for (i = 0; i < 2; i++) {
+ dvalues[LALSQTPN_OMEGA] += params->coeff.domegaSO[i] * LNhchih[i]
+ * omegaPowi_3[LAL_PNORDER_ONE_POINT_FIVE];
+ dvalues[LALSQTPN_MECO] += params->coeff.mecoSO[i] * LNhchih[i] * omegaPowi_3[2];
+ }
+ // dLNh and SO for dchih
+ for (i = 0; i < 3; i++) {
+ dvalues[LALSQTPN_CHIH1_1 + i] += params->coeff.dchihSO[0]
+ * LNhxchih[0][i] * omegaPowi_3[5];
+ dvalues[LALSQTPN_CHIH2_1 + i] += params->coeff.dchihSO[1]
+ * LNhxchih[1][i] * omegaPowi_3[5];
+ dvalues[LALSQTPN_LNH_1 + i] += (params->coeff.dLNh[0] * dvalues[LALSQTPN_CHIH1_1
+ + i] + params->coeff.dLNh[1] * dvalues[LALSQTPN_CHIH2_1 + i])
+ * omegaPowi_3[1];
+ }
+ }
+ case LAL_PNORDER_ONE:
+ case LAL_PNORDER_HALF:
+ case LAL_PNORDER_NEWTONIAN:
+ default:
+ break;
+ }
+ dvalues[LALSQTPN_OMEGA] *= params->coeff.domegaGlobal * omegaPowi_3[7]
+ * omegaPowi_3[4];
+ dvalues[LALSQTPN_PHASE] = values[LALSQTPN_OMEGA] + values[LALSQTPN_LNH_3] * (values[LALSQTPN_LNH_2]
+ * dvalues[LALSQTPN_LNH_1] - values[LALSQTPN_LNH_1] * dvalues[LALSQTPN_LNH_2])
+ / (SQR(values[LALSQTPN_LNH_1]) + SQR(values[LALSQTPN_LNH_2]));
+ return GSL_SUCCESS;
+}
+
+void LALSQTPNGenerator(LALStatus *status, LALSQTPNWave *waveform, LALSQTPNWaveformParams *params) {
+ INITSTATUS(status, "LALSQTPNGenerator", LALSQTPNWAVEFORMC);
+ ATTATCHSTATUSPTR(status);
+ ASSERT(params, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(waveform, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+
+
+ // variable declaration and initialization
+ UINT4 i = 0; // index
+ REAL8 time = 0.;
+ REAL8 LNhztol = 1.0e-8;
+ REAL8 alpha, amp, temp1, temp2;
+ const REAL8 geometrized_m_total = params->totalMass * LAL_MTSUN_SI;
+ const REAL8 freq_Step = geometrized_m_total * LAL_PI;
+ const REAL8 step = params->samplingTime / geometrized_m_total;
+ REAL8 values[LALSQTPN_NUM_OF_VAR], dvalues[LALSQTPN_NUM_OF_VAR];
+ LALSQTPNIntegratorSystem integrator;
+ xlalErrno = 0;
+ if (XLALSQTPNIntegratorInit(&integrator, LALSQTPN_NUM_OF_VAR, params, LALSQTPNDerivator)) {
+ if(XLAL_ENOMEM == XLALClearErrno()) {
+ ABORT(status, LALINSPIRALH_EMEM, LALINSPIRALH_MSGEMEM);
+ } else {
+ ABORTXLAL(status);
+ }
+ }
+
+ // initializing the dynamic variables
+ values[LALSQTPN_PHASE] = params->phi;
+ values[LALSQTPN_OMEGA] = params->lowerFreq * freq_Step;
+ values[LALSQTPN_LNH_1] = sin(params->inclination);
+ values[LALSQTPN_LNH_2] = 0.;
+ values[LALSQTPN_LNH_3] = cos(params->inclination);
+ values[LALSQTPN_MECO] = 0.;
+ for (i = 0; i < 3; i++) {
+ values[LALSQTPN_CHIH1_1 + i] = params->chih[0][i];
+ values[LALSQTPN_CHIH2_1 + i] = params->chih[1][i];
+ }
+
+ // filling the LALSQTPNCoefficients
+ xlalErrno = 0;
+ XLALSQTPNFillCoefficients(params);
+ if (xlalErrno) {
+ ABORTXLAL(status);
+ }
+ LALSQTPNDerivator(time, values, dvalues, params);
+ dvalues[LALSQTPN_MECO] = -1.; // to be able to start the loop
+ i = 0;
+ do {
+ alpha = atan2(values[LALSQTPN_LNH_2], values[LALSQTPN_LNH_1]);
+ amp = params->signalAmp * pow(values[LALSQTPN_OMEGA], 2. / 3.);
+
+ // writing the output
+ if (waveform->hp || waveform->hc || waveform->h) {
+ temp1 = -0.5*amp*cos(2.*values[LALSQTPN_PHASE])
+ * (values[LALSQTPN_LNH_3] * values[LALSQTPN_LNH_3] + 1.);
+ temp2 = amp * sin(2.*values[LALSQTPN_PHASE]) * values[LALSQTPN_LNH_3];
+ if (waveform->h) {
+ waveform->hp->data[2*i] = temp1 * cos(2.*alpha) + temp2 * sin(2.*alpha);
+ waveform->hc->data[2*i+1] = temp1 * sin(2.*alpha) - temp2 * cos(2.*alpha);
+ }
+ if (waveform->hp) {
+ waveform->hp->data[i] = temp1 * cos(2.*alpha) + temp2 * sin(2.*alpha);
+ }
+ if (waveform->hc) {
+ waveform->hc->data[i] = temp1 * sin(2.*alpha) - temp2 * cos(2.*alpha);
+ }
+ }
+ if (waveform->waveform) {
+ waveform->waveform->a->data->data[2*i] = -amp*0.5*(1. + values[LALSQTPN_LNH_3]
+ * values[LALSQTPN_LNH_3]);
+ waveform->waveform->a->data->data[2*i + 1] = -amp * values[LALSQTPN_LNH_3];
+ waveform->waveform->phi->data->data[i] = 2. * (values[LALSQTPN_PHASE] - params->phi);
+ waveform->waveform->shift->data->data[i] = 2. * alpha;
+ waveform->waveform->f->data->data[i] = values[LALSQTPN_OMEGA] / freq_Step;
+ }
+
+ // evolving
+ time = i++ * params->samplingTime;
+ xlalErrno = 0;
+ if(XLALSQTPNIntegratorFunc(values, &integrator, step)) {
+ ABORTXLAL(status);
+ }
+ // if one of the variables is nan, the PN approximation braked down
+ if (isnan(values[LALSQTPN_PHASE]) || isnan(values[LALSQTPN_OMEGA])
+ || isnan(values[LALSQTPN_LNH_1]) || isnan(values[LALSQTPN_LNH_2])
+ || isnan(values[LALSQTPN_LNH_3]) || isnan(values[LALSQTPN_CHIH1_1])
+ || isnan(values[LALSQTPN_CHIH1_2]) || isnan(values[LALSQTPN_CHIH1_3])
+ || isnan(values[LALSQTPN_CHIH2_1]) || isnan(values[LALSQTPN_CHIH2_2])
+ || isnan(values[LALSQTPN_CHIH2_3])) {
+ break;
+ }
+ LALSQTPNDerivator(time, values, dvalues, params);
+ if ((waveform->waveform && i == waveform->waveform->f->data->length) ||
+ (waveform->hp && i == waveform->hp->length) ||
+ (waveform->hc && i == waveform->hc->length)) {
+ XLALSQTPNIntegratorFree(&integrator);
+ ABORT(status, LALINSPIRALH_ESIZE, LALINSPIRALH_MSGESIZE);
+ }
+ } while (dvalues[LALSQTPN_MECO] < 0. && dvalues[LALSQTPN_OMEGA] > 0.0 && SQR(values[LALSQTPN_LNH_3])
+ < 1. - LNhztol && values[LALSQTPN_OMEGA] / freq_Step < params->samplingFreq
+ / 2.);
+ if (waveform->hp || waveform->hc){
+ params->finalFreq = values[LALSQTPN_OMEGA] / (LAL_PI * geometrized_m_total);
+ params->coalescenceTime = time;
+ }
+ if (waveform->waveform->a){
+ params->finalFreq = waveform->waveform->f->data->data[i-1];
+ }
+
+ waveform->length = i;
+ XLALSQTPNIntegratorFree(&integrator);
+ DETATCHSTATUSPTR(status);
+ RETURN(status);
+}
diff --git a/lalinspiral/src/LALSQTPNWaveform.h b/lalinspiral/src/LALSQTPNWaveform.h
new file mode 100644
index 0000000..1cd78a9
--- /dev/null
+++ b/lalinspiral/src/LALSQTPNWaveform.h
@@ -0,0 +1,225 @@
+/**
+ * @file LALSQTPNWaveform.h
+ * Contains the enumerations, structures and functions declarations to create GWforms.
+ * <b>Notations</b><br />
+ * \f$M=M_{in}M_\odot G/c^3\f$<br />
+ * \f$\hat{L_N}=\left(\sin\iota,0,\cos\iota\right)\f$ is the direction of the
+ * orbital angular momentum in radiation frame see [1].<br />
+ * \f$\omega=\pi f_L\f$<br />
+ * \f$\hat{\chi}_i\f$, \f$M_{in}\f$, \f$\iota\f$, \f$f_L\f$ are defined in
+ * LALSQTPNFillParams() function.<br />
+ * <b>References</b><br />
+ * [1] L. E. Kidder, Phys.Rev. D52, 821 (1995)<br />
+ * [2] Alessandra Buonanno, Yanbei Chen, and Michele Vallisneri, Phys.Rev. D67 (2003) 104025; Erratum-ibid. D74 (2006) 029904<br />
+ * [3] Balázs Mikóczi, Mátyás Vasúth, László Á. Gergely, Phys.Rev. D71 (2005) 124043
+ * @author László Veréb
+ * @date 2010.05.21.
+ */
+
+#ifndef LALSQTPNWAVEFORM_H
+#define LALSQTPNWAVEFORM_H
+
+#include <math.h>
+
+#include <lal/LALStatusMacros.h>
+#include <lal/LALInspiral.h>
+#include <lal/Units.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+NRCSID (LALSQTPNWAVEFORMH, "$Id$ LALSQTPNWaveform.h");
+
+typedef struct tagLALSQTPNWave{
+ CoherentGW *waveform;
+ REAL4Vector *h;
+ REAL4Vector *hp;
+ REAL4Vector *hc;
+ UINT4 length;
+} LALSQTPNWave;
+
+/** The structure contains the coefficients for calculating the derivatives
+ * of the evolving quantities.
+ */
+typedef struct tagLALSQTPNCoefficients{
+ ///@name coefficients for domega
+ //@{
+ REAL8 domegaGlobal; ///< global coefficient for domega
+ REAL8 domega[LAL_PNORDER_PSEUDO_FOUR]; ///< coefficients for domega for every PN order
+ REAL8 domegaSO[3]; ///< the spin-orbit coefficients for domega
+ REAL8 domegaSS[3]; ///< the spin1-spin2 coefficients for domega
+ REAL8 domegaSSself[3]; ///< the spin-selft coefficients for domega
+ REAL8 domegaSSselfConst; ///< the constant spin-selft coefficients for domega
+ REAL8 domegaQM[3]; ///< the quadropole-monopole coefficients for domega
+ REAL8 domegaQMConst; ///< the constant quadropole-monopole coefficients for domega
+ REAL8 domegaLN; ///< coefficient for the ln component in domega
+ //@}
+ ///@name coefficients for dchih and dLNh
+ REAL8 dchihSO[3]; ///< the spin-orbit coefficients for dchih
+ REAL8 dchihSS[3]; ///< the spin1-spin2 coefficientd for dchih
+ REAL8 dchihQM[3]; ///< the quadropole-monopole coefficients for dchih
+ REAL8 dLNh[3]; ///< coefficients for dLNh
+ //@}
+ ///@name coefficients for MECO
+ //@{
+ REAL8 meco[8]; ///< coefficients for MECO-test
+ REAL8 mecoSO[3]; ///< spin-orbit coefficients for MECO
+ REAL8 mecoSS; ///< spin1-spin2 coefficients for MECO
+ REAL8 mecoQM; ///< quadropole-monopole coefficients for MECO
+ //@}
+} LALSQTPNCoefficients;
+
+/** The structure contains the system's and the generator's parameters.
+ */
+typedef struct tagLALSQTPNWaveformParams{
+ ///@name mass-Parameters
+ //@{
+ REAL8 mass[2]; ///< masses of the BHs in \f$M_\odot\f$
+ REAL8 totalMass; ///< total mass in \f$M_\odot\f$
+ REAL8 chirpMass; ///< chirp mass in \f$M_\odot\f$
+ REAL8 mu; ///< reduced mass in \f$M_\odot\f$
+ REAL8 eta; ///< symmetric mass ratio
+ //@}
+ ///@name spin-parameters
+ //@{
+ REAL8 chi[2][3]; ///< components of the normalized spin
+ REAL8 chih[2][3]; ///< components of the unity-vectors of the normalized spin
+ REAL8 chiAmp[2]; ///< amplitude of the normalized spin
+ //@}
+ ///@name other system-parameters
+ //@{
+ REAL8 qmParameter[2]; ///< flatness of the BHs or NSs
+ REAL8 distance; ///< distance to the source in \f$Mps\f$
+ REAL8 inclination; ///< inclination of the system \f$rad\f$
+ REAL8 phi; ///< the initial phase (currently not in use)
+ //@}
+ ///@name other parameters
+ //@{
+ REAL8 signalAmp; ///< the amplitude of the signal
+ REAL8 lowerFreq; ///< the detectors sensitivityband's lower border in \f$Hz\f$
+ REAL8 samplingFreq; ///< sampling frequency in \f$Hz\f$
+ REAL8 samplingTime; ///< sampling time in \f$s\f$
+ REAL8 finalFreq; ///< the final frequency
+ REAL8 coalescenceTime; ///< the time at the coalescence
+ LALPNOrder order; ///< the Post_Newtonian order of the GW generation
+ LALSpinInteraction spinInteraction; ///< which spin interaction will be included in the generation
+ LALSQTPNCoefficients coeff; ///< coefficients for the deriving the parameters
+ //@}
+} LALSQTPNWaveformParams;
+
+/** The function fills the #LALSQTPNCoefficients structure with the needed
+ * coefficients for calculating the derived dynamic variables with the LALSQTPNDerivator() function.
+ *
+ * The orders above 2PN are incomplete, so use them if you want to try their
+ * effects.
+ * @param[in,out] params : the LALSQTPN_Generator's parameters
+ */
+void XLALSQTPNFillCoefficients(LALSQTPNWaveformParams * const params);
+
+/** The function calculates the derived values.
+ * The formulae are:
+ * \f{center}
+ * \newcommand{\OM}[1]{\left(M\omega\right)^{#1/3}}
+ * \newcommand{\derTpM}[1]{\dfrac{d#1}{d\left(t/M\right)}}
+ * \newcommand{\SPU}[2]{\left(\hat{#1}\hat{#2}\right)}
+ * \newcommand{\VPU}[2]{\left(\hat{#1}\times\hat{#2}\right)}
+ * \begin{gather}
+ * \begin{gathered}\label{eq:LALSQTPNchih}
+ * \derTpM{{\hat{\chi}}_i}={SO}_{\chi i}\OM{5}+\left({SS}_{\chi i}+{QM}_{\chi i}\right)\OM{6},\\
+ * {SO}_{\chi i}=\frac{\eta}{2}\left(4+3\frac{m_j}{m_i}\right)\left(\hat{L_N}\times\hat{\chi_i}\right),\\
+ * {SS}_{\chi i}=\frac{1}{2}\frac{\chi_jm_j^2}{M^2}\left[\hat{\chi_j}-3\left(\hat{L_N}\hat{\chi_j}\right)\hat{L_N}\right]\times\hat{\chi_i},\\
+ * {QM}_{\chi i}=-\frac{3}{2}\eta\chi_iw_i\left(\hat{L_N}\hat{\chi_i}\right)\left(\hat{L_N}\times\hat{\chi_i}\right),
+ * \end{gathered}\\[15pt]
+ * \derTpM{\hat{L_N}}=\sum_i-\frac{1}{\eta}\frac{\chi_im_i^2}{M^2}\derTpM{\hat{\chi_i}},\\[15pt]
+ * \begin{gathered}\label{eq:LALSQTPNomega}
+ * \begin{split}
+ * \derTpM{\left(M\omega\right)}&=\frac{96\eta}{5}\OM{11}\bigg[1-\frac{743+924\eta}{336}\OM{2}\bigg.\\&
+ * \quad+\left(4\pi+SO_{\omega}\right)\OM{3}+\left(\frac{34103}{18144}+\frac{13661}{2016}\eta+\frac{59}{18}\eta^2\bigg.\\&
+ * \quad+\bigg.\bigg.SS_{\omega}+SSself_{\omega}+QM_{\omega}\bigg)\OM{4}\bigg],
+ * \end{split}\\
+ * SO_{\omega}=\sum_{i\ne j}-\frac{1}{12}\frac{\chi_im_i^2}{M^2}\left(113+75\frac{m_j}{m_i}\right)\SPU{L_N}{\chi_i},\\
+ * SS_{\omega}=\frac{\eta\chi_1\chi_2}{48}\left[721\SPU{L_N}{\chi_1}\SPU{L_N}{\chi_2}-247\SPU{\chi_1}{\chi_2}\right],\\
+ * SSself_{\omega}=\sum_{i}\frac{1}{96}\frac{\chi_im_i}{M^2}\chi_i\left[7-\SPU{L_N}{\chi_i}^2\right],\\
+ * QM_{\omega}=\sum_{i}\frac{5}{2}\frac{\chi_im_i^2}{M^2}\chi_iw_i\left[3\SPU{L_N}{\chi_i}-1\right],
+ * \end{gathered}\\[15pt]
+ * \begin{gathered}
+ * \begin{split}
+ * MECO&=-0.5\eta\frac{2}{3}\OM{-1}+0.5\eta\frac{4}{3}\frac{9+\eta}{12}\OM{1}\\&\quad+
+ * SO_{MECO}\OM{2}+\Big(0.5\eta\frac{6}{3}\frac{81-57\eta+\eta^2}{24}\Big.\\&\quad+
+ * \Big.SS_{MECO}+QM_{MECO}\Big)\OM{3},
+ * \end{split}\\
+ * SO_{MECO}=\sum_{i}-\frac{5}{9}\eta\frac{\chi_im_i^2}{M^2}\left(4+3\frac{m_j}{m_i}\right)\SP{\hat{L}_N}{\hat{\chi}_i};\quad
+ * QM_{MECO}=2\eta QM_{\omega}\\
+ * SS_{MECO}=-\displaystyle\frac{\chi_1m_1^2}{M^2}\frac{\chi_2m_2^2}{M^2}\left[\SP{\hat{\chi}_1}{\hat{\chi}_2}-3\SP{\hat{L}_N}{\hat{\chi}_1}\SP{\hat{L}_N}{\hat{\chi}_2}\right]
+ * \end{gathered}\\[15pt]
+ * \begin{gathered}
+ * \derTpM{\Phi}=M\omega-\derTpM{\alpha}\cos\iota\\\alpha=\arctan\dfrac{\hat{L_N}_y}{\hat{L_N}_x}
+ * \end{gathered}
+ * \end{gather}
+ * \f}
+ * Equation (1) is given by (Eq. (2)-(3)) of [2] except of the quadropole-monopole contribution
+ * (\f$QM_{\chi_i}\f$) that was added here, equation (2) is given by (Eq. (9)) of [2],
+ * equation (3) is given by (Eq. (7)-(9)) of [3],
+ * equation (4) is given by (Eq. (11),(13)) of [2] except of the quadropole-monopole
+ * contribution (\f$QM_{MECO}\f$) that was added here, equation
+ * (5) is given by (Eq. (4.31)) of [1].<br />
+ */
+/** The equation's constant parts are calculated by the LALSQTPNFillCoefficients() function.
+ * This function is used by the LALSQTPNIntegratorFunc() function to evolve the system with a given time.
+ * @param[in] t : evolution time, not in used
+ * @param[in] values : the values to be derived
+ * @param[out] dvalues : the derived values and the last element is the MECO
+ * @param[in] params : the LALSQTPN_Generator's parameters
+ */
+int LALSQTPNDerivator(REAL8 t, const REAL8 values[], REAL8 dvalues[],
+ void * params);
+
+/** Enumeration to index the dynamic variables in the LALSQTPNGenerator function.
+ */
+typedef enum {
+ LALSQTPN_PHASE, ///< index of the phase
+ LALSQTPN_OMEGA, ///< index of the \f$M\omega\f$
+ LALSQTPN_LNH_1, ///< index of the \f$\hat{L}_N\f$'s x component
+ LALSQTPN_LNH_2, ///< index of the \f$\hat{L}_N\f$'s y component
+ LALSQTPN_LNH_3, ///< index of the \f$\hat{L}_N\f$'s z component
+ LALSQTPN_CHIH1_1, ///< index of the \f$\hat{\chi}_1\f$'s x component
+ LALSQTPN_CHIH1_2, ///< index of the \f$\hat{\chi}_1\f$'s y component
+ LALSQTPN_CHIH1_3, ///< index of the \f$\hat{\chi}_1\f$'s z component
+ LALSQTPN_CHIH2_1, ///< index of the \f$\hat{\chi}_2\f$'s x component
+ LALSQTPN_CHIH2_2, ///< index of the \f$\hat{\chi}_2\f$'s y component
+ LALSQTPN_CHIH2_3, ///< index of the \f$\hat{\chi}_2\f$'s z component
+ LALSQTPN_MECO, ///< index of the MECO
+ LALSQTPN_NUM_OF_VAR ///< number of the dynamic variables
+} LALSQTPNGeneratorVariables;
+
+/** The function generates the parameters of the waveform.
+ * The formulae are:
+ * \f{center}{
+ * \begin{gather*}
+ * \newcommand{\derTpM}[1]{\dfrac{d#1}{d\left(t/M\right)}}
+ * a_1=-0.5A\left(M\omega\right)^{2/3}\left(1+\cos^2\iota\right),\quad
+ * a_2=-A\left(M\omega\right)^{2/3}\cos\iota
+ * \end{gather}
+ * \f}
+ *
+ * The waveform-parameters are generated by evolving the system with the LALSQTPNIntegratorFunc() function.
+ * The derived values of the dynamic variables [\f$\left(M\omega\right)\f$, \f$\hat{L_N}_x\f$,
+ * \f$\hat{L_N}_y\f$, \f$\hat{L_N}_z\f$, \f$\Phi\f$] are calculated by the LALSQTPNDerivator() function.
+ * The function returns the amplitude \f$a_1\f$, \f$a_2\f$, polarization
+ * shift \f$\alpha\f$, and phase \f$\Phi\f$.
+ * The \f$\alpha\f$ is defined with equation (5) in the documentation of the
+ * LALSQTPNDerivator() function.
+ * @param[in,out] status : LAL universal status structure
+ * @param[out] waveform : the generated waveform
+ * @param[in] params : the input parameters
+ */
+void
+LALSQTPNGenerator(LALStatus *status, LALSQTPNWave *waveform,
+ LALSQTPNWaveformParams *params);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* LALSQTPN_WAVEFORM_H */
diff --git a/lalinspiral/src/LALSQTPNWaveformInterface.c b/lalinspiral/src/LALSQTPNWaveformInterface.c
new file mode 100644
index 0000000..272f545
--- /dev/null
+++ b/lalinspiral/src/LALSQTPNWaveformInterface.c
@@ -0,0 +1,296 @@
+/**
+ * @file LALSQTPNWaveformInterface.c
+ * Contains function definitions to integrate the SpinQuadTaylor code into the other parts of the LALSuit.
+ * If you want to run the program use the LALSQTPNWaveformTest.c file int the
+ * test directory.
+ * @author László Veréb
+ * @date 2010.06.27.
+ */
+
+#include <lal/LALSQTPNWaveformInterface.h>
+#include <lal/LALSQTPNWaveform.h>
+
+NRCSID (LALSQTPNWAVEFORMINTERFACEC, "$Id LALSQTPNWaveformInterface.c$");
+
+/** The macro function returns the square of the argument.
+ * Do not use with incrementing or decrementing operators!
+ * @param[in] a : the number
+ * @return the square of the number
+ */
+#define SQR(a) ((a)*(a))
+
+void LALSQTPNWaveformTemplates (LALStatus *status, REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *params) {
+
+ InspiralInit paramsInit;
+
+ INITSTATUS(status, "LALSTPNWaveform", LALSQTPNWAVEFORMINTERFACEC);
+ ATTATCHSTATUSPTR(status);
+
+ ASSERT(signalvec1, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(signalvec1->data, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(signalvec2, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(signalvec2->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);
+ LALSQTPNWaveformParams wave_Params;
+ LALSQTPNWave wave;
+
+ TRY(LALInspiralInit(status->statusPtr, params, &paramsInit), status);
+
+ memset(signalvec1->data, 0, signalvec1->length * sizeof(REAL4));
+ memset(signalvec2->data, 0, signalvec2->length * sizeof(REAL4));
+ XLALSQTPNFillParams(&wave_Params, params);
+
+ wave.waveform = NULL;
+ wave.h = NULL;
+ wave.hp = signalvec1;
+ wave.hc = signalvec2;
+
+ /* Call the engine function */
+ LALSQTPNGenerator(status->statusPtr, &wave, &wave_Params);
+ CHECKSTATUSPTR(status);
+
+ DETATCHSTATUSPTR(status);
+ RETURN(status);
+}
+
+void LALSQTPNWaveform (LALStatus *status, REAL4Vector *signalvec, InspiralTemplate *params){
+ INITSTATUS(status, "LALSQTPNWaveform", LALSQTPNWAVEFORMINTERFACEC);
+ ATTATCHSTATUSPTR(status);
+ InspiralInit paramsInit;
+ LALSQTPNWaveformParams wave_Params;
+ LALSQTPNWave wave;
+ memset(&wave, 0, sizeof(LALSQTPNWave));
+ wave.h = NULL;
+ wave.hp = signalvec;
+ wave.hc = NULL;
+ wave.waveform = NULL;
+
+ ASSERT(signalvec, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(signalvec->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);
+ XLALSQTPNFillParams(&wave_Params, params);
+ wave_Params.distance *= LAL_PC_SI * 1.e6;
+ wave_Params.signalAmp /= LAL_PC_SI * 1.e6;
+
+ /* Call the engine function */
+ LALSQTPNGenerator(status->statusPtr, &wave, &wave_Params);
+ params->tC = wave_Params.coalescenceTime;
+ CHECKSTATUSPTR(status);
+ DETATCHSTATUSPTR(status);
+ RETURN(status);
+}
+
+void LALSQTPNWaveformForInjection(LALStatus *status, CoherentGW *waveform,
+ InspiralTemplate *params, PPNParamStruc *ppnParams) {
+ INITSTATUS(status, "LALSQTPNWaveformInterface", LALSQTPNWAVEFORMINTERFACEC);
+ ATTATCHSTATUSPTR(status);
+ // variable declaration and initialization
+ UINT4 count, i;
+ InspiralInit paramsInit;
+
+ ASSERT(params, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(waveform, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(!(waveform->a), status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(!(waveform->f), status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(!(waveform->phi), status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ ASSERT(!(waveform->shift), status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
+ // Compute some parameters
+ LALInspiralInit(status->statusPtr, params, &paramsInit);
+ CHECKSTATUSPTR(status);
+ if (paramsInit.nbins == 0) {
+ DETATCHSTATUSPTR(status);
+ RETURN(status);
+ }
+
+ // Allocate the waveform structures.
+ xlalErrno = 0;
+ XLALSQTPNAllocateCoherentGW(waveform, paramsInit.nbins);
+
+ LALSQTPNWave wave;
+ wave.h = NULL;
+ wave.hp = NULL;
+ wave.hc = NULL;
+ wave.waveform = waveform;
+ LALSQTPNWaveformParams wave_Params;
+
+ // filling the parameters
+ XLALSQTPNFillParams(&wave_Params, params);
+
+ // calling the engine function
+ LALSQTPNGenerator(status->statusPtr, &wave, &wave_Params);
+ BEGINFAIL(status) {
+ XLALSQTPNDestroyCoherentGW(waveform);
+ } ENDFAIL(status);
+ params->fFinal = wave_Params.finalFreq;
+ count = wave.length;
+ for (i = 0; i < count; i++) {
+ if (waveform->phi->data->data[i] != 0.) {
+ break;
+ }
+ if (i == count - 1) {
+ XLALSQTPNDestroyCoherentGW(waveform);
+ DETATCHSTATUSPTR( status );
+ RETURN( status );
+ }
+ }
+
+ {
+ if (waveform->a != NULL) {
+ waveform->f->data->length = waveform->phi->data->length = waveform->shift->data->length = count;
+ waveform->a->data->length = 2 * count;
+ for (i = 0; i < wave.length; i++) {
+ // (PPNParamStruct)ppnParams->phi === (InspiralTemplate)params->startPhase === (SimInspiralTable)injparams/this_event->coa_phase it is set to 0 in LALSQTPNWaveformTest.c at line 83.
+ waveform->phi->data->data[i] = waveform->phi->data->data[i] - waveform->phi->data->data[wave.length-1] + ppnParams->phi;
+ }
+ waveform->a->deltaT = waveform->f->deltaT = waveform->phi->deltaT
+ = waveform->shift->deltaT = 1. / params->tSampling;
+
+ waveform->a->sampleUnits = lalStrainUnit;
+ waveform->f->sampleUnits = lalHertzUnit;
+ waveform->phi->sampleUnits = lalDimensionlessUnit;
+ waveform->shift->sampleUnits = lalDimensionlessUnit;
+
+ waveform->position = ppnParams->position;
+ waveform->psi = ppnParams->psi;
+
+ snprintf(waveform->a->name, LALNameLength,
+ "STPN inspiral amplitudes");
+ snprintf(waveform->f->name, LALNameLength,
+ "STPN inspiral frequency");
+ snprintf(waveform->phi->name, LALNameLength, "STPN inspiral phase");
+ snprintf(waveform->shift->name, LALNameLength,
+ "STPN inspiral polshift");
+ }
+ // --- fill some output ---
+ ppnParams->tc = (REAL8) (count - 1) / params->tSampling;
+ ppnParams->length = count;
+ ppnParams->dfdt = ((REAL4) (waveform->f->data->data[count - 1]
+ - waveform->f->data->data[count - 2])) * ppnParams->deltaT;
+ ppnParams->fStop = params->fFinal;
+ ppnParams->termCode = GENERATEPPNINSPIRALH_EFSTOP;
+ ppnParams->termDescription = GENERATEPPNINSPIRALH_MSGEFSTOP;
+
+ ppnParams->fStart = ppnParams->fStartIn;
+ } // end phase condition
+ DETATCHSTATUSPTR(status);
+ RETURN(status);
+}
+
+int XLALSQTPNAllocateCoherentGW(CoherentGW *wave, UINT4 length) {
+ static const char *func = "LALSQTPNAllocateCoherentGW";
+ if (!wave) {
+ XLAL_ERROR(func, XLAL_EFAULT);
+ }
+ if (length <= 0) {
+ XLAL_ERROR(func, XLAL_EBADLEN);
+ }
+ if (wave->a || wave->f || wave->phi || wave->shift) {
+ XLAL_ERROR(func, XLAL_EFAULT);
+ }
+ wave->a = (REAL4TimeVectorSeries *)LALMalloc(sizeof(REAL4TimeVectorSeries));
+ wave->f = (REAL4TimeSeries *)LALMalloc(sizeof(REAL4TimeSeries));
+ wave->phi = (REAL8TimeSeries *)LALMalloc(sizeof(REAL8TimeSeries));
+ wave->shift = (REAL4TimeSeries *)LALMalloc(sizeof(REAL4TimeSeries));
+ if (!(wave->a && wave->f && wave->phi && wave->shift)) {
+ XLALSQTPNDestroyCoherentGW(wave);
+ XLAL_ERROR(func, XLAL_ENOMEM);
+ }
+ xlalErrno = 0;
+ wave->a->data = XLALCreateREAL4VectorSequence(length, 2);
+ wave->f->data = XLALCreateREAL4Vector(length);
+ wave->phi->data = XLALCreateREAL8Vector(length);
+ wave->shift->data = XLALCreateREAL4Vector(length);
+ if (!(wave->a->data && wave->f->data && wave->phi->data && wave->shift->data)) {
+ XLALSQTPNDestroyCoherentGW(wave);
+ XLAL_ERROR(func, XLAL_ENOMEM);
+ }
+ return XLAL_SUCCESS;
+}
+
+void XLALSQTPNDestroyCoherentGW(CoherentGW *wave) {
+ //static const char *func = "LALSQTPNDestroyCoherentGW";
+ if (wave->a) {
+ if (wave->a->data) {
+ XLALDestroyREAL4VectorSequence(wave->a->data);
+ }
+ XLALFree(wave->a);
+
+ }
+ if (wave->f) {
+ if (wave->f->data) {
+ XLALDestroyREAL4Vector(wave->f->data);
+ }
+ XLALFree(wave->f);
+ }
+ if (wave->phi) {
+ if (wave->phi->data) {
+ XLALDestroyREAL8Vector(wave->phi->data);
+ }
+ XLALFree(wave->phi);
+ }
+ if (wave->shift) {
+ if (wave->shift->data) {
+ XLALDestroyREAL4Vector(wave->shift->data);
+ }
+ XLALFree(wave->shift);
+ }
+}
+
+void XLALSQTPNFillParams(LALSQTPNWaveformParams *wave, InspiralTemplate *params) {
+ wave->mass[0] = params->mass1;
+ wave->mass[1] = params->mass2;
+ wave->totalMass = wave->mass[0] + wave->mass[1];
+ wave->mu = wave->mass[0] * wave->mass[1] / wave->totalMass;
+ wave->eta = wave->mu / wave->totalMass;
+ wave->chirpMass = wave->totalMass * pow(wave->eta, 3. / 5.);
+ wave->chiAmp[0] = wave->chiAmp[1] = 0.;
+ INT2 i;
+ for (i = 0; i < 3; i++) {
+ wave->chi[0][i] = params->spin1[i];
+ wave->chi[1][i] = params->spin2[i];
+ wave->chiAmp[0] += SQR(wave->chi[0][i]);
+ wave->chiAmp[1] += SQR(wave->chi[1][i]);
+ }
+ wave->chiAmp[0] = sqrt(wave->chiAmp[0]);
+ wave->chiAmp[1] = sqrt(wave->chiAmp[1]);
+ for (i = 0; i < 3; i++) {
+ if (wave->chiAmp[0] != 0.) {
+ wave->chih[0][i] = wave->chi[0][i] / wave->chiAmp[0];
+ } else {
+ wave->chih[0][i] = 0.;
+ }
+ if (wave->chiAmp[1] != 0.) {
+ wave->chih[1][i] = wave->chi[1][i] / wave->chiAmp[1];
+ } else {
+ wave->chih[1][i] = 0.;
+ }
+ }
+ wave->qmParameter[0] = params->qmParameter[0];
+ wave->qmParameter[1] = params->qmParameter[1];
+ wave->distance = params->distance;
+ wave->inclination = params->inclination;
+ wave->lowerFreq = params->fLower;
+ wave->samplingFreq = params->tSampling;
+ wave->samplingTime = 1. / wave->samplingFreq;
+ wave->phi = 0.;
+ wave->signalAmp = 4. * wave->totalMass * wave->eta * LAL_MRSUN_SI / wave->distance;
+ wave->order = params->order;
+ wave->spinInteraction = params->spinInteraction;
+ if (wave->spinInteraction != 0)
+ wave->spinInteraction |= LAL_SOInter;
+}
+
diff --git a/lalinspiral/src/LALSQTPNWaveformInterface.h b/lalinspiral/src/LALSQTPNWaveformInterface.h
new file mode 100644
index 0000000..15302a8
--- /dev/null
+++ b/lalinspiral/src/LALSQTPNWaveformInterface.h
@@ -0,0 +1,104 @@
+/**
+ * @file LALSQTPNWaveformInterface.h
+ * Contains function declarations to integrate the SpinQuadTaylor code into the other parts of the LALSuit.
+ * @author László Veréb
+ * @date 2010.06.27.
+ */
+
+#ifndef LALSQTPNWAVEFORMINTERFACE_H
+#define LALSQTPNWAVEFORMINTERFACE_H
+
+#include <lal/Units.h>
+#include <lal/LALInspiral.h>
+#include <lal/SeqFactories.h>
+#include <lal/LALStatusMacros.h>
+#include <lal/GenerateInspiral.h>
+#include <lal/LALSQTPNWaveform.h>
+#include <stdio.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+NRCSID (LALSQTPNWAVEFORMINTERFACEH, "$Id LALSQTPNWaveformInterface.h$");
+
+#define LALSQTPN_MSGPPNPARAMS "the PPNParamsStruct structure is null"
+#define LALSQTPN_MSGINSPIRALTEMPLATE "the InspiralTemplate structure is null"
+#define LALSQTPN_ZEROLENGTH 3
+#define LALSQTPN_MSGZEROLENGTH "the given length is not positive"
+
+void LALSQTPNWaveformTemplates (LALStatus *status, REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *params);
+
+/** The function returns the generated waveform.
+ * @param[in,out] status : LAL universal status structure
+ * @param[out] signalvec : array containing the waveform \f$(h_+, h_\times)\f$
+ * @param[in] params : structure containing the inspiral parameters
+ */
+void LALSQTPNWaveform (LALStatus *status, REAL4Vector *signalvec, InspiralTemplate *params);
+
+/** The function returns the generated waveform for injection.
+ * @param[in,out] status : LAL universal status structure
+ * @param[out] wave_out : structure containing the waveform \f$(a_1,
+ * a_2, \Phi, \alpha)\f$
+ * @param[in] params : structure containing the inspiral parameters
+ * @param[in] ppnParams : parameters for restricted post-Newtonian waveform
+ */
+void LALSQTPNWaveformForInjection(LALStatus *status, CoherentGW *wave_out,
+ InspiralTemplate *params, PPNParamStruc *ppnParams);
+
+/** The function allocates memory for the waveform's \f$a_1\f$, \f$a_2\f$,
+ * \f$\Phi\f$ and \f$\alpha\f$.
+ * @param[out] waveform : pointer to the allocated waveform
+ * @param[in] length : the length of the waveform
+ */
+int XLALSQTPNAllocateCoherentGW(CoherentGW *wave, UINT4 length);
+
+/** The function deallocates memory of the waveform.
+ * @param[out] waveform : pointer to the allocated waveform
+ */
+void XLALSQTPNDestroyCoherentGW(CoherentGW *wave);
+
+/** The function calculates the parameters from the InspiralTemplate
+ * structure. <em>The used parameters are:</em>
+ * <ul>
+ * <li>masses of the BHs (or NSs) \f$m_i\f$ in \f$M_\odot\f$</li>
+ * <li>the spin components \f$\chi_{ij}\f$, the values of \f$\sqrt{\sum_j\chi_{ij}}\f$, are between 0 and 1</li>
+ * <li>the quadrupole parameters \f$w_i\in(4,8)\f$ for NSs [1] and \f$w_i=1\f$ for BHs[2] are 1 (default 1)</li>
+ * <li>the inclination (angle between the line of sight and Newtonian orbital angular momentum) \f$\iota\f$ in \f$rad\f$
+ * <li>the initial frequency \f$f_L\f$ in \f$Hz\f$</li>
+ * <li>the distance \f$d\f$ in \f$Mpc\f$</li>
+ * <li>the sampling time \f$t_s\f$ in \f$s\f$</li>
+ * <li>the PN order, see #LALPNOrder</li>
+ * <li>level of accuracy in including spin and quadrupole contributions, see
+ * #LALSQTPNSpinInteraction</li>
+ * <ul><br />
+ * <em>The calculated parameters:</em>
+ * \f{center}{
+ * \begin{gather}
+ * \displaystyle M_{in}=m_1+m_2,\quad
+ * \mu=\frac{m_1m_2}{M_{in}},\quad
+ * \eta=\frac{\mu}{M_{in}},\\
+ * \chi_i=\sqrt{\sum_{j}\chi_{ij}^2},\quad
+ * \hat{\chi}_{ij}=\dfrac{\chi_{ij}}{\chi_i},\\
+ * f_s=t_s^{-1}\\
+ * A=\frac{4\cdot\eta M_{in}M_\odot\displaystyle\frac{G}{c^2}}{d\cdot3.0856775807\cdot10^{16}\cdot10^6}
+ * \end{gather}
+ * \f}
+ * and the initial phase \f$\phi=0\f$
+ * Assuming that:
+ * <ul>
+ * <li>masses are positive</li>
+ * <li>eta is positive</li>
+ * <li>sampling frequency is positive</li>
+ * <li>distance is positive</li>
+ * </ul>
+ * @param[out] wave : the used parameters
+ * @param[in] params : the inspiral parameters
+ */
+void XLALSQTPNFillParams(LALSQTPNWaveformParams *wave, InspiralTemplate *params);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* LALSQTPNWAVEFOMRINTERFACE_H */
diff --git a/lalinspiral/src/Makefile.am b/lalinspiral/src/Makefile.am
index 265cc3c..c8522e6 100644
--- a/lalinspiral/src/Makefile.am
+++ b/lalinspiral/src/Makefile.am
@@ -30,6 +30,9 @@ lalinspiralinclude_HEADERS = \
LALNoiseModelsInspiral.h \
LALSBBH-Metric.h \
LALSimInspiral.h \
+ LALSQTPNIntegrator.h \
+ LALSQTPNWaveform.h \
+ LALSQTPNWaveformInterface.h \
LIGOLwXMLInspiralRead.h \
LIGOLwXMLRingdownRead.h \
LIGOMetadataInspiralUtils.h \
@@ -173,6 +176,9 @@ liblalinspiral_la_SOURCES = \
LALRectangleVertices.c \
LALRungeKutta4.c \
LALSTPNWaveform.c \
+ LALSQTPNIntegrator.c \
+ LALSQTPNWaveform.c \
+ LALSQTPNWaveformInterface.c \
LALSimInspiral.c \
LALSimInspiralPNMode.c \
LALTrigScanCluster.c \
diff --git a/lalinspiral/test/LALSQTPNWaveformTest.c b/lalinspiral/test/LALSQTPNWaveformTest.c
new file mode 100644
index 0000000..e3532b9
--- /dev/null
+++ b/lalinspiral/test/LALSQTPNWaveformTest.c
@@ -0,0 +1,139 @@
+/**
+ * @file LALSQTPNWaveformTest.c
+ * The user interface for the SpinQuadTaylor program.
+ * This file is an example howto use the SpinQuadTaylor program.
+ * The input parameters are:<br/>
+ * <em>Waveform parameters:</em>(note:\f$i=1,2\f$ and \f$j=x,y,z\f$)
+ * <ul>
+ * <li>masses of the BHs (or NSs) \f$m_i\f$ in \f$M_\odot\f$</li>
+ * <li>the spin components \f$\chi_{ij}\f$, the values of \f$\sqrt{\sum_j\chi_{ij}}\f$, are between 0 and 1</li>
+ * <li>the quadrupole parameters \f$w_i\in(4,8)\f$ for NSs [1] and \f$w_i=1\f$ for BHs[2] are 1 (default 1)</li>
+ * <li>the inclination (angle between the line of sight and Newtonian orbital angular momentum) \f$\iota\f$ in \f$rad\f$</li>
+ * <li>the initial frequency \f$f_L\f$ in \f$Hz\f$</li>
+ * <li>the distance \f$d\f$ in \f$Mpc\f$</li>
+ * <li>the sampling time \f$t_s\f$ in \f$s\f$</li>
+ * </ul>
+ *<em>Program parameters:</em>
+ * <ul>
+ * <li>the name of the output file (default out.dat)</li>
+ * <li>the PN order (newtonian, oneHalfPN, onePN, onePointFivePN, twoPN(default), twoPointFivePN, threePN, threePointFivePN)</li>
+ * <li>level of accuracy in including spin and quadrupole contributions
+ * (NO, SO, SS, SELF, QM, ALL(default))</li>
+ * </ul>
+ * The output file contains three coloums: elapsed time, \f$h_+\f$, \f$h_\times\f$.
+ * \f{center}
+ * \begin{gather*}
+ * h_+=a_1\cos\left(2\alpha\right)\cos\left(2\Phi\right)-a_2\sin\left(2\alpha\right)\sin\left(2\Phi\right),\\
+ * h_\times=a_1\sin\left(2\alpha\right)\cos\left(2\Phi\right)+a_2\cos\left(2\alpha\right)\sin\left(2\Phi\right)
+ * \end{gather*}
+ * \f}
+ * with \f$a_i\f$ amplitudes, \f$\alpha\f$ polarization shift, \f$\Phi\f$ phase (Eq. (4.28)-(4.30) of [3] up to leading order (The \f$\Phi\f$ is shifted by \f$\pi\f$ with respect to [3]). We note that \f$\Theta=0\f$ in leading order because we use radiation frame).<br />
+ * \f$a_1\f$, \f$a_2\f$ are defined in LALSQTPNGenerator() function, \f$\alpha\f$ and \f$\Phi\f$ phase is defined in LALSQTPNDerivator() function.<br />
+ * <b>References</b><br />
+ * [1] E. Poisson, Phys.Rev. D57 5287 (1998)<br />
+ * [2] K. S. Thorne, Rev.Mod.Phys. 52 299 (1980)<br />
+ * [3] L. E. Kidder, Phys.Rev.D 52, 821 (1995)<br />
+ * @author László Veréb
+ * @date 2010.06.27.
+ */
+
+#include <lal/GenerateInspiral.h>
+#include <lal/LALStdlib.h>
+#include <lal/LALInspiral.h>
+#include <lal/GeneratePPNInspiral.h>
+#include <lal/LALSQTPNWaveformInterface.h>
+
+NRCSID(LALSQTPNWAVEFORMTESTC, "$Id: LALSQTPNWaveformTest.c,v 0.1 2010/05/21");
+
+int lalDebugLevel = 0; ///< the debug level
+
+/** The main program.
+ */
+int main(int argc, char *argv[]) {
+
+ // variable declaration and initialization
+ static LALStatus mystatus;
+ CoherentGW thewaveform;
+ SimInspiralTable injParams;
+ PPNParamStruc ppnParams;
+ const char *filename;
+ FILE *outputfile;
+ INT4 i, length;
+ REAL8 dt;
+ char PNString[50];
+
+ if (argc == 17) {
+ sprintf(PNString, "SpinQuadTaylor%s%s", argv[15], argv[16]);
+ filename = argv[14];
+ } else if (argc == 16) {
+ sprintf(PNString, "SpinQuadTaylor%sALL", argv[15]);
+ filename = argv[14];
+ } else if (argc == 15) {
+ sprintf(PNString, "SpinQuadTaylorTwoPNALL");
+ filename = argv[14];
+ } else if (argc == 14) {
+ sprintf(PNString, "SpinQuadTaylorTwoPNALL");
+ filename = "out.dat";
+ } else if (argc != 14) {
+ printf(
+ " 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16\n");
+ printf(
+ "Correct parameter order: m1 m2 S1x S1y S1z S2x S2y S2z incl f_lower f_final distance dt output PNorder SpinInter\n");
+ return (1);
+ }
+
+ memset(&mystatus, 0, sizeof(LALStatus));
+ memset(&thewaveform, 0, sizeof(CoherentGW));
+ memset(&injParams, 0, sizeof(SimInspiralTable));
+ memset(&ppnParams, 0, sizeof(PPNParamStruc));
+
+ // setting the parameters
+ injParams.mass1 = atof(argv[1]);
+ injParams.mass2 = atof(argv[2]);
+ injParams.spin1x = atof(argv[3]);
+ injParams.spin1y = atof(argv[4]);
+ injParams.spin1z = atof(argv[5]);
+ injParams.spin2x = atof(argv[6]);
+ injParams.spin2y = atof(argv[7]);
+ injParams.spin2z = atof(argv[8]);
+ injParams.qmParameter[0] = 1.;//atof(argv[9]);
+ injParams.qmParameter[1] = 1.;//atof(argv[10]);
+ injParams.inclination = atof(argv[9]);
+ injParams.f_lower = atof(argv[10]);
+ injParams.f_final = atof(argv[11]);
+ injParams.distance = atof(argv[12]);
+ ppnParams.deltaT = atof(argv[13]);
+ injParams.polarization = 0;
+ LALSnprintf(injParams.waveform, LIGOMETA_WAVEFORM_MAX * sizeof(CHAR),
+ PNString);
+
+ //interface(&mystatus, &thewaveform, &injParams, &ppnParams);
+ LALGenerateInspiral(&mystatus, &thewaveform, &injParams, &ppnParams);
+ if (mystatus.statusCode) {
+ fprintf( stderr, "LALSQTPNWaveformTest: error generating waveform\n" );
+ return mystatus.statusCode;
+ }
+ // and finally save in a file
+ outputfile = fopen(filename, "w");
+
+ length = thewaveform.f->data->length;
+ dt = thewaveform.phi->deltaT;
+ REAL8 a1, a2, phi, shift;
+ for(i = 0; i < length; i++) {
+ a1 = thewaveform.a->data->data[2*i];
+ a2 = thewaveform.a->data->data[2*i+1];
+ phi = thewaveform.phi->data->data[i] - thewaveform.phi->data->data[0];
+ shift = thewaveform.shift->data->data[i];
+
+ fprintf(outputfile,"% 10.7e\t% 10.7e\t% 10.7e\n",
+ i*dt,
+ a1*cos(shift)*cos(phi) - a2*sin(shift)*sin(phi),
+ a1*sin(shift)*cos(phi) + a2*cos(shift)*sin(phi));
+ }
+ fclose(outputfile);
+ XLALSQTPNDestroyCoherentGW(&thewaveform);
+ puts("Done.");
+ LALCheckMemoryLeaks();
+ return mystatus.statusCode;
+}
+
diff --git a/lalmetaio/src/LIGOMetadataTables.h b/lalmetaio/src/LIGOMetadataTables.h
index 3187253..4343631 100644
--- a/lalmetaio/src/LIGOMetadataTables.h
+++ b/lalmetaio/src/LIGOMetadataTables.h
@@ -658,6 +658,7 @@ tagSimInspiralTable
INT4 amp_order;
CHAR taper[LIGOMETA_INSPIRALTAPER_MAX];
INT4 bandpass;
+ REAL4 qmParameter[2];
}
SimInspiralTable;
/* </lalVerbatim> */