21 #include <lal/LALStdio.h>
22 #include <lal/LALStdlib.h>
23 #include <lal/LALConstants.h>
24 #include <lal/Units.h>
25 #include <lal/AVFactories.h>
26 #include <lal/SeqFactories.h>
27 #include <lal/PulsarSimulateCoherentGW.h>
28 #include <lal/GenerateSpinOrbitCW.h>
158 REAL8 upsilon, argument;
167 REAL8 eccMinusOne, eccPlusOne;
179 GENERATESPINORBITCWH_MSGENUL );
181 GENERATESPINORBITCWH_MSGENUL );
185 GENERATESPINORBITCWH_MSGEOUT );
187 GENERATESPINORBITCWH_MSGEOUT );
189 GENERATESPINORBITCWH_MSGEOUT );
191 GENERATESPINORBITCWH_MSGEOUT );
196 GENERATESPINORBITCWH_MSGENUL );
197 nSpin =
params->f->length;
203 eccMinusOne = -
params->oneMinusEcc;
204 ecc = 1.0 + eccMinusOne;
205 eccPlusOne = 2.0 + eccMinusOne;
206 if ( eccMinusOne <= 0.0 ) {
208 GENERATESPINORBITCWH_MSGEECC );
211 vDotAvg =
params->angularSpeed
212 * sqrt( eccMinusOne * eccMinusOne * eccMinusOne / eccPlusOne );
218 GENERATESPINORBITCWH_MSGEFTL );
220 if ( vp <= 0.0 ||
dt <= 0.0 ||
f0 <= 0.0 || vDotAvg <= 0.0 ||
223 GENERATESPINORBITCWH_MSGESGN );
230 a = vp * eccMinusOne * cos( argument ) + ecc;
231 b = -vp * sqrt( eccMinusOne / eccPlusOne ) * sin( argument );
232 eCosOmega = ecc * cos( argument );
234 dxMax = 0.01 / (
f0 *
n *
dt );
238 if ( dxMax < 1.0e-15 ) {
241 " precision for this orbit" );
248 if (
f0 * tau * vp * vp * ecc / eccPlusOne > 0.25 )
250 " effects that are not included" );
255 tPeriObs = (
REAL8 )(
params->orbitEpoch.gpsSeconds -
257 tPeriObs += 1.0e-9 * (
REAL8 )(
params->orbitEpoch.gpsNanoSeconds -
260 spinOff = (
REAL8 )(
params->orbitEpoch.gpsSeconds -
261 params->spinEpoch.gpsSeconds );
262 spinOff += 1.0e-9 * (
REAL8 )(
params->orbitEpoch.gpsNanoSeconds -
263 params->spinEpoch.gpsNanoSeconds );
267 xMinus = 1.0 +
a * sinh( -1.0 ) + b *
cosh( -1.0 ) - b;
268 xPlus = -1.0 +
a * sinh( 1.0 ) + b *
cosh( 1.0 ) - b;
269 x = -vDotAvg * tPeriObs;
271 e = -log( -2.0 * (
x - xMinus ) / (
a - b ) - exp( 1.0 ) );
272 }
else if (
x <= 0 ) {
274 }
else if (
x <= xPlus ) {
277 e = log( 2.0 * (
x - xPlus ) / (
a + b ) - exp( 1.0 ) );
280 coshe =
cosh(
e ) - 1.0;
286 GENERATESPINORBITCWH_MSGEMEM );
294 GENERATESPINORBITCWH_MSGEMEM );
304 GENERATESPINORBITCWH_MSGEMEM );
370 fData =
output->f->data->data;
371 phiData =
output->phi->data->data;
372 for (
i = 0;
i <
n;
i++ ) {
374 x = vDotAvg * (
i *
dt - tPeriObs );
379 while ( fabs( dx = log(
e -
a * sinhe - b * coshe ) -
x ) > dxMax ) {
382 coshe =
cosh(
e ) - 1.0;
384 }
else if (
x > xPlus ) {
386 while ( fabs( dx = log( -
e +
a * sinhe + b * coshe ) -
x ) > dxMax ) {
389 coshe =
cosh(
e ) - 1.0;
395 while ( fabs( dx = -
e +
a * sinhe + b * coshe -
x ) > dxMax ) {
396 e -= dx / ( -1.0 +
a * coshe +
a + b * sinhe );
399 }
else if (
e > 1.0 ) {
403 coshe =
cosh(
e ) - 1.0;
409 ( ecc * sinhe -
e ) / vDotAvg + spinOff;
411 for (
j = 0;
j < nSpin;
j++ ) {
412 f += fSpin[
j] * tPow;
413 phi += fSpin[
j] * ( tPow *=
t ) / (
j + 2.0 );
417 upsilon = 2.0 * atan2( sqrt( eccPlusOne * coshe ),
418 sqrt( eccMinusOne * ( coshe + 2.0 ) ) );
419 f *=
f0 / ( 1.0 + vp * ( cos( argument + upsilon ) + eCosOmega )
422 if ( fabs(
f - fPrev ) >
df ) {
423 df = fabs(
f - fPrev );
425 *( fData++ ) = fPrev =
f;
426 *( phiData++ ) = phi +
phi0;
#define ABORT(statusptr, code, mesg)
#define ENDFAIL(statusptr)
#define TRY(func, statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define ASSERT(assertion, statusptr, code, mesg)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define BEGINFAIL(statusptr)
void LALSCreateVectorSequence(LALStatus *status, REAL4VectorSequence **vectorSequence, CreateVectorSequenceIn *vSeqParams)
#define GENERATESPINORBITCWH_EECC
Eccentricity out of range.
void LALGenerateHyperbolicSpinOrbitCW(LALStatus *stat, PulsarCoherentGW *output, SpinOrbitCWParamStruc *params)
Computes a continuous waveform with frequency drift and Doppler modulation from a hyperbolic orbital ...
#define GENERATESPINORBITCWH_ESGN
Sign error: positive parameter expected.
#define GENERATESPINORBITCWH_EMEM
Out of memory.
#define GENERATESPINORBITCWH_ENUL
Unexpected null pointer in arguments.
#define GENERATESPINORBITCWH_EOUT
Output field a, f, phi, or shift already exists.
#define GENERATESPINORBITCWH_EFTL
Periapsis motion is faster than light.
int LALWarning(LALStatus *status, const char *warning)
const LALUnit lalStrainUnit
const LALUnit lalHertzUnit
const LALUnit lalDimensionlessUnit
void LALDCreateVector(LALStatus *, REAL8Vector **, UINT4)
void LALDDestroyVector(LALStatus *, REAL8Vector **)
void LALSDestroyVector(LALStatus *, REAL4Vector **)
void LALSCreateVector(LALStatus *, REAL4Vector **, UINT4)
This structure stores a representation of a plane gravitational wave propagating from a particular po...
This structure stores the parameters for constructing a gravitational waveform with both a Taylor-pol...