21 #include <lal/LALStdio.h>
22 #include <lal/LALStdlib.h>
23 #include <lal/LALError.h>
24 #include <lal/DetectorSite.h>
25 #include <lal/DetResponse.h>
26 #include <lal/AVFactories.h>
28 #include <lal/Units.h>
29 #include <lal/TimeDelay.h>
30 #include <lal/LALBarycenter.h>
31 #include <lal/VectorOps.h>
32 #include <lal/PulsarSimulateCoherentGW.h>
33 #include <lal/SkyCoordinates.h>
34 #include <lal/SinCosLUT.h>
38 #define SIMULATECOHERENTGWH_ENUL 1
39 #define SIMULATECOHERENTGWH_EBAD 2
40 #define SIMULATECOHERENTGWH_ESIG 3
41 #define SIMULATECOHERENTGWH_EDIM 4
42 #define SIMULATECOHERENTGWH_EMEM 5
43 #define SIMULATECOHERENTGWH_EUNIT 6
47 #define SIMULATECOHERENTGWH_MSGENUL "Unexpected null pointer in arguments"
48 #define SIMULATECOHERENTGWH_MSGEBAD "A sampling interval is (effectively) zero"
49 #define SIMULATECOHERENTGWH_MSGESIG "Input signal must specify amplitude and phase functions"
50 #define SIMULATECOHERENTGWH_MSGEDIM "Amplitude must be a 2-dimensional vector"
51 #define SIMULATECOHERENTGWH_MSGEMEM "Memory allocation error"
52 #define SIMULATECOHERENTGWH_MSGEUNIT "Bad input units"
95 #define TCENTRE( time ) \
97 realIndex = delayOff + (time)*delayDt, \
98 intIndex = (INT4)floor( realIndex ), \
99 indexFrac = realIndex - intIndex, \
100 time + indexFrac*delayData[intIndex+1] \
101 + (1.0-indexFrac)*delayData[intIndex] \
241 INT4 shiftInit, shiftFinal;
245 REAL8 delayMin, delayMax;
265 REAL4 *aData, *fData, *shiftData, *plusData, *crossData;
266 REAL8 *phiData, *delayData;
267 REAL8 aOff, fOff, phiOff, shiftOff, polOff, delayOff;
268 REAL8 aDt, fDt, phiDt, shiftDt, polDt, delayDt;
276 REAL4 *aTransData = NULL, *phiTransData = NULL;
278 REAL8 phiFac = 1.0, fFac = 1.0;
295 SIMULATECOHERENTGWH_MSGENUL );
296 if ( !( CWsignal->
a ) ) {
298 SIMULATECOHERENTGWH_MSGESIG );
304 if ( !( CWsignal->
phi ) ) {
306 SIMULATECOHERENTGWH_MSGESIG );
318 if ( CWsignal->
shift ) {
326 if ( ( transfer = (
detector->transfer != NULL ) ) ) {
363 if ( CWsignal->
shift ) {
375 fFac = 1.0 /
detector->transfer->deltaF;
381 detector->heterodyneEpoch.gpsSeconds );
383 detector->heterodyneEpoch.gpsNanoSeconds );
386 LALWarning(
stat,
"REAL8 arithmetic is not sufficient to maintain"
387 " heterodyne phase to within a radian." );
396 if ( CWsignal->
shift ) {
415 outData =
output->data->data;
416 INT4 fLen = 0, shiftLen = 0;
424 if ( CWsignal->
shift ) {
442 location.
x =
detector->site->location[0];
443 location.
y =
detector->site->location[1];
444 location.
z =
detector->site->location[2];
456 delayDt =
output->deltaT / ( 2.0 * dtDelayBy2 );
457 nMax = (
UINT4 )(
output->data->length * delayDt ) + 3;
459 delayData = delay->
data;
473 for (
i = 0;
i < 3;
i++ ) {
483 for (
i = 0;
i < nMax;
i++ ) {
495 if ( tDelay < delayMin ) {
498 if ( tDelay > delayMax ) {
508 LALInfo(
stat,
"Detector site and ephemerides absent; simulating hplus with no"
509 " propagation delays" );
510 memset( delayData, 0, nMax *
sizeof(
REAL8 ) );
511 delayMin = delayMax = 0.0;
515 polDt =
output->deltaT / ( 2.0 * dtPolBy2 );
516 nMax = (
UINT4 )(
output->data->length * polDt ) + 3;
526 if ( polResponse.
pPlus ) {
529 if ( polResponse.
pCross ) {
537 SIMULATECOHERENTGWH_MSGEMEM );
579 if ( plusLen != crossLen ) {
580 XLALPrintError(
"plusLen = %d != crossLen = %d\n", plusLen, crossLen );
596 params.deltaT = 2.0 * dtPolBy2;
617 for (
i = 0;
i < nMax;
i++ ) {
629 INT4 phiTransLen = 0, aTransLen = 0;
631 nMax =
detector->transfer->data->length;
708 phiTransData = phiTransfer->
data;
709 phiTransLen = phiTransfer->
length;
710 aTransData = aTransfer->
data;
711 aTransLen = aTransfer->
length;
713 if ( aTransLen != phiTransLen ) {
714 XLALPrintError(
"aTransLen = %d != phiTransLen = %d\n", aTransLen, phiTransLen );
720 aOff = (
output->epoch.gpsSeconds -
723 aOff += (
output->epoch.gpsNanoSeconds -
726 phiOff = (
output->epoch.gpsSeconds -
729 phiOff += (
output->epoch.gpsNanoSeconds -
733 fOff = (
output->epoch.gpsSeconds -
736 fOff += (
output->epoch.gpsNanoSeconds -
742 if ( CWsignal->
shift ) {
743 shiftOff = (
output->epoch.gpsSeconds -
746 shiftOff += (
output->epoch.gpsNanoSeconds -
758 if ( aOff + (
i + delayMin )*aDt < 0.0 ) {
759 INT4 j = (
INT4 )floor( -aOff / aDt - delayMax );
763 while ( (
i < (
INT4 )(
output->data->length ) ) &&
764 ( aOff +
TCENTRE(
i )*aDt < 0.0 ) ) {
768 if ( phiOff + (
i + delayMin )*phiDt < 0.0 ) {
769 INT4 j = (
INT4 )( -phiOff / phiDt - delayMax );
773 while ( (
i < (
INT4 )(
output->data->length ) ) &&
774 ( phiOff +
TCENTRE(
i )*phiDt < 0.0 ) ) {
788 if ( aOff + (
n + delayMax )*aDt > nMax_a ) {
789 INT4 j = (
INT4 )( ( nMax_a - aOff ) / aDt - delayMin + 1.0 );
793 while ( (
n >= 0 ) &&
799 if ( phiOff + (
n + delayMax )*phiDt > nMax ) {
800 INT4 j = (
INT4 )( ( nMax - phiOff ) / phiDt - delayMin + 1.0 );
804 while ( (
n >= 0 ) &&
805 ( (
INT4 )floor( phiOff +
TCENTRE(
n )*phiDt ) > nMax ) ) {
818 if ( fOff + ( fInit + delayMin )*fDt < 0.0 ) {
819 INT4 j = (
INT4 )floor( -fOff / fDt - delayMax );
823 while ( ( fInit <=
n ) &&
824 ( fOff +
TCENTRE( fInit )*fDt < 0.0 ) ) {
830 if ( fOff + ( fFinal + delayMax )*fDt > nMax ) {
831 INT4 j = (
INT4 )( ( nMax - fOff ) / fDt - delayMin + 1.0 );
835 while ( ( fFinal >=
i ) &&
836 ( (
INT4 )floor( fOff +
TCENTRE( fFinal )*fDt ) > nMax ) ) {
846 if ( CWsignal->
shift ) {
848 if ( shiftOff + ( shiftInit + delayMin )*shiftDt < 0.0 ) {
849 INT4 j = (
INT4 )floor( -shiftOff / shiftDt - delayMax );
850 if ( shiftInit <
j ) {
853 while ( ( shiftInit <=
n ) &&
854 ( shiftOff +
TCENTRE( shiftInit )*shiftDt < 0.0 ) ) {
860 if ( shiftOff + ( shiftFinal + delayMax )*shiftDt > nMax ) {
861 INT4 j = (
INT4 )( ( nMax - shiftOff ) / shiftDt - delayMin + 1.0 );
862 if ( shiftFinal >
j ) {
865 while ( ( shiftFinal >=
i ) &&
866 ( (
INT4 )floor( shiftOff +
TCENTRE( shiftFinal )*shiftDt ) > nMax ) ) {
879 if ( ( nMax =
output->data->length -
n - 1 ) > 0 ) {
880 memset(
output->data->data +
n + 1, 0, nMax *
sizeof(
REAL4 ) );
886 nMax =
detector->transfer->data->length - 1;
890 for ( ;
i <=
n;
i++ ) {
899 REAL4 aTrans, phiTrans;
904 x = aOff + iCentre * aDt;
908 if (
j + 3 >= aLen ) {
909 XLALPrintError(
"Interpolation error: no point at or after last sample for {a1,a2}: i = %d, n = %d, j = %d, aLen = %d",
i,
n,
j, aLen );
913 a1 = frac * aData[
j + 2] + ( 1.0 - frac ) * aData[
j];
914 a2 = frac * aData[
j + 3] + ( 1.0 - frac ) * aData[
j + 1];
918 if ( (
i < shiftInit ) || (
i > shiftFinal ) ) {
921 x = shiftOff + iCentre * shiftDt;
925 if (
j + 1 >= shiftLen ) {
926 XLALPrintError(
"Interpolation error: no point at or after last sample for 'shift'" );
929 shift = frac * shiftData[
j + 1] + ( 1.0 - frac ) * shiftData[
j];
933 x = phiOff + iCentre * phiDt;
937 if (
j + 1 >= phiLen ) {
938 XLALPrintError(
"Interpolation error: no point at or after last sample for 'phi'" );
941 phi = frac * phiData[
j + 1] + ( 1.0 - frac ) * phiData[
j];
942 phi -= heteroFac *
i +
phi0;
946 if ( (
i < fInit ) || (
i > fFinal ) ) {
947 f = ( phiData[
j + 1] - phiData[
j] ) * phiFac;
950 x = fOff + iCentre * fDt;
953 if (
j + 1 >= fLen ) {
954 XLALPrintError(
"Interpolation error: no point at or after last sample for 'f'" );
957 f = frac * fData[
j + 1] + ( 1.0 - frac ) * fData[
j];
961 if ( (
x < 0.0 ) || (
x > nMax ) ) {
968 if (
j + 1 >= aTransLen ) {
969 XLALPrintError(
"Interpolation error: no point at or after last sample for '{aTrans,phiTrans}'" );
973 aTrans = frac * aTransData[
j + 1] + ( 1.0 - frac ) * aTransData[
j];
974 phiTrans = frac * phiTransData[
j + 1] + ( 1.0 - frac ) * phiTransData[
j];
989 oPlus = a1 * cs *
cp -
a2 *
ss * sp;
990 oCross = a1 *
ss *
cp +
a2 * cs * sp;
995 x = polOff +
i * polDt;
998 if (
j + 1 >= plusLen ) {
999 XLALPrintError(
"Interpolation error: no point at or after last sample for '{oPlus,oCross}'" );
1003 oPlus *= frac * plusData[
j + 1] + ( 1.0 - frac ) * plusData[
j];
1004 oCross *= frac * crossData[
j + 1] + ( 1.0 - frac ) * crossData[
j];
1007 outData[
i] = oPlus + oCross;
1013 LALWarning(
stat,
"Signal passed outside of the frequency domain"
1014 " of the transfer function (transfer function is"
1015 " treated as zero outside its specified domain)" );
1017 LALInfo(
stat,
"Signal frequency was estimated by differencing"
1018 " the signal phase" );
#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)
#define SIMULATECOHERENTGWH_EUNIT
Bad input units.
#define SIMULATECOHERENTGWH_ENUL
Unexpected null pointer in arguments.
#define SIMULATECOHERENTGWH_ESIG
Input signal must specify amplitude and phase functions.
#define SIMULATECOHERENTGWH_EMEM
Memory allocation error.
#define SIMULATECOHERENTGWH_EDIM
Amplitude must be a 2-dimensional vector.
#define SIMULATECOHERENTGWH_EBAD
A sampling interval is (effectively) zero.
void LALComputeDetAMResponseSeries(LALStatus *status, LALDetAMResponseSeries *pResponseSeries, const LALDetAndSource *pDetAndSource, const LALTimeIntervalAndNSample *pTimeInfo)
int XLALBarycenterEarth(EarthState *earth, const LIGOTimeGPS *tGPS, const EphemerisData *edat)
Computes the position and orientation of the Earth, at some arrival time , specified LIGOTimeGPS inp...
int XLALBarycenter(EmissionTime *emit, const BarycenterInput *baryinput, const EarthState *earth)
Transforms from detector arrival time in GPS (as specified in the LIGOTimeGPS structure) to pulse em...
#define XLAL_INIT_DECL(var,...)
int LALWarning(LALStatus *status, const char *warning)
int LALInfo(LALStatus *status, const char *info)
int XLALPulsarSimulateCoherentGW(REAL4TimeSeries *output, PulsarCoherentGW *CWsignal, PulsarDetectorResponse *detector)
FIXME: Temporary XLAL-wapper to LAL-function LALPulsarSimulateCoherentGW()
void LALPulsarSimulateCoherentGW(LALStatus *stat, REAL4TimeSeries *output, PulsarCoherentGW *CWsignal, PulsarDetectorResponse *detector)
Computes the response of a detector to a coherent gravitational wave.
int XLALSinCosLUT(REAL4 *sinx, REAL4 *cosx, REAL8 x)
Calculate sin(x) and cos(x) to roughly 1e-7 precision using a lookup-table and Taylor-expansion.
void LALConvertSkyCoordinates(LALStatus *stat, SkyPosition *output, SkyPosition *input, ConvertSkyParams *params)
COORDINATESYSTEM_EQUATORIAL
void LALGeocentricToGeodetic(LALStatus *stat, EarthPosition *location)
int XLALUnitCompare(const LALUnit *unit1, const LALUnit *unit2)
const LALUnit lalHertzUnit
const LALUnit lalDimensionlessUnit
LALUnit * XLALUnitMultiply(LALUnit *output, const LALUnit *unit1, const LALUnit *unit2)
void LALDCreateVector(LALStatus *, REAL8Vector **, UINT4)
void LALDDestroyVector(LALStatus *, REAL8Vector **)
void LALSDestroyVector(LALStatus *, REAL4Vector **)
void LALSCreateVector(LALStatus *, REAL4Vector **, UINT4)
void LALCVectorAngle(LALStatus *status, REAL4Vector *out, const COMPLEX8Vector *in)
void LALCVectorAbs(LALStatus *status, REAL4Vector *out, const COMPLEX8Vector *in)
void LALUnwrapREAL4Angle(LALStatus *status, REAL4Vector *out, const REAL4Vector *in)
#define XLAL_CHECK(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
Basic output structure of LALBarycenterEarth.c.
Basic output structure produced by LALBarycenter.c.
REAL8 deltaT
(TDB) - (GPS)
REAL4TimeSeries * pScalar
const LALDetector * pDetector
SkyPosition equatorialCoords
This structure stores a representation of a plane gravitational wave propagating from a particular po...
REAL4TimeSeries * f
A time-sampled sequence storing the instantaneous frequency , in Hz.
REAL4 psi
The polarization angle , in radians, as defined in Appendix B of .
SkyPosition position
The location of the source in the sky; this should be in equatorial celestial coordinates,...
REAL4TimeVectorSeries * a
A time-sampled two-dimensional vector storing the amplitudes and , in dimensionless strain.
REAL4TimeSeries * shift
A time-sampled sequence storing the polarization shift , in radians.
UINT4 dtPolBy2
A user defined half-interval time step for the polarisation response look-up table (will default to 3...
UINT4 dtDelayBy2
A user specified half-interval time step for the Doppler delay look-up table (will default to 400s if...
REAL8TimeSeries * phi
A time-sampled sequence storing the phase function , in radians.
This structure contains information required to determine the response of a detector to a gravitation...
REAL4VectorSequence * data