30 #include <gsl/gsl_sf_expint.h>
31 #include <lal/LALSimulation.h>
32 #include <lal/LALDetectors.h>
33 #include <lal/DetResponse.h>
35 #include <lal/Units.h>
36 #include <lal/TimeDelay.h>
37 #include <lal/SkyCoordinates.h>
38 #include <lal/TimeSeries.h>
39 #include <lal/TimeSeriesInterp.h>
40 #include <lal/FrequencySeries.h>
41 #include <lal/TimeFreqFFT.h>
42 #include <lal/Window.h>
65 for(n = 1; n && (
x & (
x + 1)); n *= 2)
95 XLALPrintError(
"%s(): error: can't identify instrument from string \"%s\"\n", __func__,
string);
110 REAL8TimeSeries * XLALSimQuasiPeriodicInjectionREAL8TimeSeries(
REAL8TimeSeries *aplus,
REAL8TimeSeries *across,
REAL8TimeSeries *frequency,
REAL8TimeSeries *phi,
REAL8TimeSeries *psi,
SkyPosition *position,
LALDetector *
detector, EphemerisData *ephemerides,
LIGOTimeGPS *start,
REAL8 deltaT,
UINT4 length,
COMPLEX16FrequencySeries *response )
112 REAL8 slowDeltaT = 600.0;
123 slowLength = floor(0.5 + (length*
deltaT + 16.0)/slowDeltaT);
136 for ( j = 0; j < injection->length; ++j ) {
155 injection->data->data[j] = fpl*apl*cos(phase) + fcr*acr*sin(phase);
182 double T = kernel_data->
T;
184 double *stop = cached_kernel + kernel_length;
187 for(
i = -(kernel_length - 1) / 2; cached_kernel < stop;
i++) {
188 double x =
i + residual;
192 double Si2 = gsl_sf_Si(
LAL_PI * (
x +
T));
193 double Si3 = gsl_sf_Si(
LAL_PI * (
x -
T));
196 *cached_kernel++ = 0.;
257 REAL8 right_ascension,
267 const int kernel_length = 67 + 48 * lround(2.0 * arm_length_samples);
269 const unsigned det_resp_interval = round(0.25 / hplus->
deltaT) < 1 ? 1 : round(0.25 / hplus->
deltaT);
272 LALREAL8TimeSeriesInterp *xinterp = NULL;
273 LALREAL8TimeSeriesInterp *yinterp = NULL;
305 sprintf(
name,
"%s injection",
detector->frDetector.prefix);
372 if(!(
i % det_resp_interval)) {
376 XLALComputeDetAMResponseParts(&armlen, &xcos, &ycos, &fxplus, &fyplus, &fxcross, &fycross,
detector, right_ascension, declination, psi,
XLALGreenwichMeanSiderealTime(&t));
391 if(!xinterp || !yinterp)
407 if(!(
i % det_resp_interval)) {
411 XLALComputeDetAMResponseParts(&armlen, &xdata.
armcos, &ydata.
armcos, &fxplus, &fyplus, &fxcross, &fycross,
detector, right_ascension, declination, psi,
XLALGreenwichMeanSiderealTime(&t));
486 const double noop_threshold = 1
e-4;
491 const unsigned aperiodicity_suppression_buffer = 32768;
492 double start_sample_int;
493 double start_sample_frac;
503 XLALPrintError(
"%s(): error: input sample rates or heterodyne frequencies do not match\n", __func__);
517 if(start_sample_frac < -0.5) {
518 start_sample_frac += 1.0;
519 start_sample_int -= 1.0;
520 }
else if(start_sample_frac > +0.5) {
521 start_sample_frac -= 1.0;
522 start_sample_int += 1.0;
527 if(fabs(start_sample_frac) > noop_threshold || response) {
540 if(i < h->
data->length) {
542 XLALPrintError(
"%s(): error: source time series too long\n", __func__);
555 if(!tilde_h || !plan) {
571 const double f = tilde_h->
f0 +
i * tilde_h->
deltaF;
596 int j = floor((f - response->
f0) / response->
deltaF + 0.5);
599 else if((
unsigned) j > response->
data->
length - 1)
621 if(tilde_h->
f0 == 0.0)
628 if(tilde_h->
f0 == 0.0)
656 XLALPrintError(
"%s(): error: oops, internal sample rate mismatch\n", __func__);
736 const double noop_threshold = 1
e-4;
741 const unsigned aperiodicity_suppression_buffer = 32768;
742 double start_sample_int;
743 double start_sample_frac;
753 XLALPrintError(
"%s(): error: input sample rates or heterodyne frequencies do not match\n", __func__);
767 if(start_sample_frac < -0.5) {
768 start_sample_frac += 1.0;
769 start_sample_int -= 1.0;
770 }
else if(start_sample_frac > +0.5) {
771 start_sample_frac -= 1.0;
772 start_sample_int += 1.0;
777 if(fabs(start_sample_frac) > noop_threshold || response) {
789 if(i < h->
data->length) {
791 XLALPrintError(
"%s(): error: source time series too long\n", __func__);
804 if(!tilde_h || !plan) {
820 const double f = tilde_h->
f0 +
i * tilde_h->
deltaF;
843 int j = floor((f - response->
f0) / response->
deltaF + 0.5);
846 else if((
unsigned) j > response->
data->
length - 1)
867 if(tilde_h->
f0 == 0.0)
874 if(tilde_h->
f0 == 0.0)
901 XLALPrintError(
"%s(): error: oops, internal sample rate mismatch\n", __func__);
955 REAL8FFTPlan *fwdplan,
956 REAL8FFTPlan *revplan,
993 &fxcross, &fycross,
detector, ra, dec, psi, gmst);
1010 if (offrac < -0.5) {
1013 }
else if (offrac > 0.5) {
1026 for (j = 0; j < (int)segment->
data->
length; ++j)
1027 if (j >= offset && j < (
int)hplus->
data->
length + offset) {
1036 for (j = 0; j < (int)segment->
data->
length; ++j)
1037 if (j >= offset && j < (
int)hcross->
data->
length + offset) {
1049 double f = work1->
f0 + k * work1->
deltaF;
1062 gplus = Tx * fxplus + Ty * fyplus;
1063 gcross = Tx * fxcross + Ty * fycross;
1099 REAL4FFTPlan *fwdplan,
1100 REAL4FFTPlan *revplan,
1137 &fxcross, &fycross,
detector, ra, dec, psi, gmst);
1154 if (offrac < -0.5) {
1157 }
else if (offrac > 0.5) {
1170 for (j = 0; j < (int)segment->
data->
length; ++j)
1171 if (j >= offset && j < (
int)hplus->
data->
length + offset) {
1180 for (j = 0; j < (int)segment->
data->
length; ++j)
1181 if (j >= offset && j < (
int)hcross->
data->
length + offset) {
1193 double f = work1->
f0 + k * work1->
deltaF;
1206 gplus = Tx * fxplus + Ty * fyplus;
1207 gcross = Tx * fxcross + Ty * fycross;
1261 const double nominal_segdur = 2.0;
1262 const double max_time_delay = 0.1;
1263 const size_t strides_per_segment = 2;
1276 REAL8FFTPlan *fwdplan = NULL;
1277 REAL8FFTPlan *revplan = NULL;
1289 if (response == NULL) {
1306 stride = seglen / strides_per_segment;
1307 padlen = max_time_delay / target->
deltaT;
1309 ovrlap -= 2 * padlen;
1357 nsteps = ((length%stride) ? (1 + length/stride) : (length/stride));
1368 if (!work1 || !work2) {
1377 if (!fwdplan || !revplan) {
1389 for (step = 0; step < nsteps; ++step) {
1407 for (j = padlen; j < seglen - padlen; ++j)
1409 if (step && j - padlen < ovrlap) {
1411 double x = (double)(j - padlen) / ovrlap;
1413 + (1.0 -
x) * h->
data->
data[j + offset];
1425 for (j = 0; j < stride - padlen; ++j)
1427 for ( ; j < stride; ++j) {
1428 double fac = window->
data->
data[j - (stride - padlen)];
1484 const double nominal_segdur = 2.0;
1485 const double max_time_delay = 0.1;
1486 const size_t strides_per_segment = 2;
1499 REAL4FFTPlan *fwdplan = NULL;
1500 REAL4FFTPlan *revplan = NULL;
1512 if (response == NULL) {
1529 stride = seglen / strides_per_segment;
1530 padlen = max_time_delay / target->
deltaT;
1532 ovrlap -= 2 * padlen;
1580 nsteps = ((length%stride) ? (1 + length/stride) : (length/stride));
1591 if (!work1 || !work2) {
1600 if (!fwdplan || !revplan) {
1612 for (step = 0; step < nsteps; ++step) {
1630 for (j = padlen; j < seglen - padlen; ++j)
1632 if (step && j - padlen < ovrlap) {
1634 double x = (double)(j - padlen) / ovrlap;
1636 + (1.0 -
x) * h->
data->
data[j + offset];
1648 for (j = 0; j < stride - padlen; ++j)
1650 for ( ; j < stride; ++j) {
1651 double fac = window->
data->
data[j - (stride - padlen)];
1695 REAL8FFTPlan *fwdplan,
1696 REAL8FFTPlan *revplan,
1728 (
const REAL4(*)[3])(uintptr_t)
detector->response, ra, dec, psi, gmst);
1745 if (offrac < -0.5) {
1748 }
else if (offrac > 0.5) {
1762 for (j = 0; j < (int)segment->
data->
length; ++j)
1763 if (j >= offset && j < (
int)hplus->
data->
length + offset)
1765 * ( fplus * hplus->
data->
data[j - offset]
1766 + fcross * hcross->
data->
data[j - offset] );
1779 double f = work->
f0 + k * work->
deltaF;
1819 REAL4FFTPlan *fwdplan,
1820 REAL4FFTPlan *revplan,
1852 (
const REAL4(*)[3])(uintptr_t)
detector->response, ra, dec, psi, gmst);
1869 if (offrac < -0.5) {
1872 }
else if (offrac > 0.5) {
1886 for (j = 0; j < (int)segment->
data->
length; ++j)
1887 if (j >= offset && j < (
int)hplus->
data->
length + offset)
1889 * ( fplus * hplus->
data->
data[j - offset]
1890 + fcross * hcross->
data->
data[j - offset] );
1903 double f = work->
f0 + k * work->
deltaF;
1966 const double nominal_segdur = 2.0;
1967 const double max_time_delay = 0.1;
1968 const size_t strides_per_segment = 2;
1980 REAL8FFTPlan *fwdplan = NULL;
1981 REAL8FFTPlan *revplan = NULL;
1993 if (response == NULL) {
2010 stride = seglen / strides_per_segment;
2011 padlen = max_time_delay / target->
deltaT;
2013 ovrlap -= 2 * padlen;
2061 nsteps = ((length%stride) ? (1 + length/stride) : (length/stride));
2079 if (!fwdplan || !revplan) {
2091 for (step = 0; step < nsteps; ++step) {
2100 hplus, hcross, work, fwdplan, revplan, window, ra, dec,
2110 for (j = padlen; j < seglen - padlen; ++j)
2112 if (step && j - padlen < ovrlap) {
2114 double x = (double)(j - padlen) / ovrlap;
2116 + (1.0 -
x) * h->
data->
data[j + offset];
2128 for (j = 0; j < stride - padlen; ++j)
2130 for ( ; j < stride; ++j) {
2131 double fac = window->
data->
data[j - (stride - padlen)];
2191 const double nominal_segdur = 2.0;
2192 const double max_time_delay = 0.1;
2193 const size_t strides_per_segment = 2;
2205 REAL4FFTPlan *fwdplan = NULL;
2206 REAL4FFTPlan *revplan = NULL;
2218 if (response == NULL) {
2235 stride = seglen / strides_per_segment;
2236 padlen = max_time_delay / target->
deltaT;
2238 ovrlap -= 2 * padlen;
2286 nsteps = ((length%stride) ? (1 + length/stride) : (length/stride));
2304 if (!fwdplan || !revplan) {
2316 for (step = 0; step < nsteps; ++step) {
2325 hplus, hcross, work, fwdplan, revplan, window, ra, dec,
2335 for (j = padlen; j < seglen - padlen; ++j)
2337 if (step && j - padlen < ovrlap) {
2339 double x = (double)(j - padlen) / ovrlap;
2341 + (1.0 -
x) * h->
data->
data[j + offset];
2353 for (j = 0; j < stride - padlen; ++j)
2355 for ( ; j < stride; ++j) {
2356 double fac = window->
data->
data[j - (stride - padlen)];
static double beta(const double a, const double b, const sysq *system)
Internal function that computes the spin-orbit couplings.
static int XLALSimComputeStrainSegmentREAL8TimeSeries(REAL8TimeSeries *segment, const REAL8TimeSeries *hplus, const REAL8TimeSeries *hcross, COMPLEX16FrequencySeries *work1, COMPLEX16FrequencySeries *work2, REAL8FFTPlan *fwdplan, REAL8FFTPlan *revplan, REAL8Window *window, double ra, double dec, double psi, LALDetector *detector, const COMPLEX16FrequencySeries *response)
static unsigned long round_up_to_power_of_two(unsigned long x)
static int XLALSimComputeLWLStrainSegmentREAL8TimeSeries(REAL8TimeSeries *segment, const REAL8TimeSeries *hplus, const REAL8TimeSeries *hcross, COMPLEX16FrequencySeries *work, REAL8FFTPlan *fwdplan, REAL8FFTPlan *revplan, REAL8Window *window, double ra, double dec, double psi, LALDetector *detector, const COMPLEX16FrequencySeries *response)
static int XLALSimComputeStrainSegmentREAL4TimeSeries(REAL4TimeSeries *segment, const REAL4TimeSeries *hplus, const REAL4TimeSeries *hcross, COMPLEX8FrequencySeries *work1, COMPLEX8FrequencySeries *work2, REAL4FFTPlan *fwdplan, REAL4FFTPlan *revplan, REAL4Window *window, double ra, double dec, double psi, LALDetector *detector, const COMPLEX8FrequencySeries *response)
static int XLALSimComputeLWLStrainSegmentREAL4TimeSeries(REAL4TimeSeries *segment, const REAL4TimeSeries *hplus, const REAL4TimeSeries *hcross, COMPLEX8FrequencySeries *work, REAL4FFTPlan *fwdplan, REAL4FFTPlan *revplan, REAL4Window *window, double ra, double dec, double psi, LALDetector *detector, const COMPLEX8FrequencySeries *response)
static void highfreq_kernel(double *cached_kernel, int kernel_length, double residual, void *data)
This kernel function incorporates beyond-long-wavelength effects.
static REAL8TimeSeries * error(const REAL8TimeSeries *s1, const REAL8TimeSeries *s0)
LALREAL8TimeSeriesInterp * XLALREAL8TimeSeriesInterpCreate(const REAL8TimeSeries *series, int kernel_length, void(*kernel)(double *, int, double, void *), void *kernel_data)
void XLALREAL8TimeSeriesInterpDestroy(LALREAL8TimeSeriesInterp *interp)
REAL8 XLALREAL8TimeSeriesInterpEval(LALREAL8TimeSeriesInterp *interp, const LIGOTimeGPS *t, int bounds_check)
#define LAL_CHECK_VALID_SERIES(s, val)
#define LAL_CHECK_CONSISTENT_TIME_SERIES(s1, s2, val)
#define LAL_CHECK_COMPATIBLE_TIME_SERIES(s1, s2, val)
const LALDetector lalCachedDetectors[LAL_NUM_DETECTORS]
void XLALComputeDetAMResponseParts(double *armlen, double *xcos, double *ycos, double *fxplus, double *fyplus, double *fxcross, double *fycross, const LALDetector *detector, double ra, double dec, double psi, double gmst)
void XLALComputeDetAMResponse(double *fplus, double *fcross, const REAL4 D[3][3], const double ra, const double dec, const double psi, const double gmst)
COMPLEX16 XLALComputeDetArmTransferFunction(double beta, double mu)
COMPLEX8FrequencySeries * XLALCreateCOMPLEX8FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
COMPLEX16FrequencySeries * XLALCreateCOMPLEX16FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
void XLALDestroyCOMPLEX8FrequencySeries(COMPLEX8FrequencySeries *series)
void XLALDestroyCOMPLEX16FrequencySeries(COMPLEX16FrequencySeries *series)
void * XLALMalloc(size_t n)
REAL8TimeSeries * XLALSimDetectorStrainREAL8TimeSeries(const REAL8TimeSeries *hplus, const REAL8TimeSeries *hcross, REAL8 right_ascension, REAL8 declination, REAL8 psi, const LALDetector *detector)
Transforms the waveform polarizations into a detector strain.
int XLALSimInjectLWLDetectorStrainREAL4TimeSeries(REAL4TimeSeries *target, const REAL4TimeSeries *hplus, const REAL4TimeSeries *hcross, double ra, double dec, double psi, LALDetector *detector, const COMPLEX8FrequencySeries *response)
Computes strain for a detector and injects into target time series.
int XLALSimInjectDetectorStrainREAL8TimeSeries(REAL8TimeSeries *target, const REAL8TimeSeries *hplus, const REAL8TimeSeries *hcross, double ra, double dec, double psi, LALDetector *detector, const COMPLEX16FrequencySeries *response)
Computes strain for a detector and injects into target time series.
int XLALSimInjectLWLDetectorStrainREAL8TimeSeries(REAL8TimeSeries *target, const REAL8TimeSeries *hplus, const REAL8TimeSeries *hcross, double ra, double dec, double psi, LALDetector *detector, const COMPLEX16FrequencySeries *response)
Computes strain for a detector and injects into target time series.
int XLALSimInjectDetectorStrainREAL4TimeSeries(REAL4TimeSeries *target, const REAL4TimeSeries *hplus, const REAL4TimeSeries *hcross, double ra, double dec, double psi, LALDetector *detector, const COMPLEX8FrequencySeries *response)
Computes strain for a detector and injects into target time series.
int XLALSimAddInjectionREAL4TimeSeries(REAL4TimeSeries *target, REAL4TimeSeries *h, const COMPLEX8FrequencySeries *response)
Adds a detector strain time series to detector data.
const LALDetector * XLALDetectorPrefixToLALDetector(const char *string)
Turn a detector prefix string into a LALDetector structure.
int XLALSimAddInjectionREAL8TimeSeries(REAL8TimeSeries *target, REAL8TimeSeries *h, const COMPLEX16FrequencySeries *response)
Adds a detector strain time series to detector data.
REAL4FFTPlan * XLALCreateReverseREAL4FFTPlan(UINT4 size, int measurelvl)
void XLALDestroyREAL8FFTPlan(REAL8FFTPlan *plan)
REAL4FFTPlan * XLALCreateForwardREAL4FFTPlan(UINT4 size, int measurelvl)
REAL8FFTPlan * XLALCreateReverseREAL8FFTPlan(UINT4 size, int measurelvl)
void XLALDestroyREAL4FFTPlan(REAL4FFTPlan *plan)
REAL8FFTPlan * XLALCreateForwardREAL8FFTPlan(UINT4 size, int measurelvl)
REAL8 XLALTimeDelayFromEarthCenter(const double detector_earthfixed_xyz_metres[3], double source_right_ascension_radians, double source_declination_radians, const LIGOTimeGPS *gpstime)
int XLALREAL8TimeFreqFFT(COMPLEX16FrequencySeries *freq, const REAL8TimeSeries *tser, const REAL8FFTPlan *plan)
int XLALREAL4FreqTimeFFT(REAL4TimeSeries *tser, const COMPLEX8FrequencySeries *freq, const REAL4FFTPlan *plan)
int XLALREAL8FreqTimeFFT(REAL8TimeSeries *tser, const COMPLEX16FrequencySeries *freq, const REAL8FFTPlan *plan)
int XLALREAL4TimeFreqFFT(COMPLEX8FrequencySeries *freq, const REAL4TimeSeries *tser, const REAL4FFTPlan *plan)
REAL8TimeSeries * XLALResizeREAL8TimeSeries(REAL8TimeSeries *series, int first, size_t length)
REAL4TimeSeries * XLALResizeREAL4TimeSeries(REAL4TimeSeries *series, int first, size_t length)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
REAL8TimeSeries * XLALAddREAL8TimeSeries(REAL8TimeSeries *arg1, const REAL8TimeSeries *arg2)
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
void XLALDestroyREAL4TimeSeries(REAL4TimeSeries *series)
REAL4TimeSeries * XLALAddREAL4TimeSeries(REAL4TimeSeries *arg1, const REAL4TimeSeries *arg2)
REAL4TimeSeries * XLALCreateREAL4TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
const LALUnit lalDimensionlessUnit
REAL4Window * XLALCreateTukeyREAL4Window(UINT4 length, REAL4 beta)
void XLALDestroyREAL8Window(REAL8Window *window)
void XLALDestroyREAL4Window(REAL4Window *window)
REAL8Window * XLALCreateTukeyREAL8Window(UINT4 length, REAL8 beta)
#define XLAL_ERROR_NULL(...)
#define XLAL_REAL8_FAIL_NAN
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_IS_REAL8_FAIL_NAN(val)
REAL8 XLALGreenwichMeanSiderealTime(const LIGOTimeGPS *gpstime)
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
REAL8 XLALGPSModf(REAL8 *iptr, const LIGOTimeGPS *epoch)
int XLALGPSCmp(const LIGOTimeGPS *t0, const LIGOTimeGPS *t1)
REAL8 XLALGPSDiff(const LIGOTimeGPS *t1, const LIGOTimeGPS *t0)