23 #include <gsl/gsl_linalg.h>
24 #include <gsl/gsl_rng.h>
25 #include <gsl/gsl_randist.h>
27 #include <lal/LALConstants.h>
28 #include <lal/LALDetectors.h>
30 #include <lal/FrequencySeries.h>
31 #include <lal/Sequence.h>
32 #include <lal/TimeFreqFFT.h>
33 #include <lal/Units.h>
34 #include <lal/LALSimSGWB.h>
36 #include <lal/LALSimReadData.h>
45 # define CLEANUP_AND_RETURN(errnum) do { \
46 if (htilde) for (i = 0; i < numDetectors; ++i) XLALDestroyCOMPLEX16FrequencySeries(htilde[i]); \
47 XLALFree(htilde); XLALDestroyREAL8FFTPlan(plan); gsl_matrix_free(R); \
48 if (errnum) XLAL_ERROR(errnum); else return 0; \
50 REAL8FFTPlan *plan = NULL;
61 deltaF = 1.0 / (length * h[0]->
deltaT);
62 psdfac = 0.3 * pow(H0 /
LAL_PI, 2.0);
89 for (k = 1; k < length/2; ++k) {
90 double f = k * deltaF;
91 double sigma = 0.5 * sqrt(psdfac * OmegaGW->
data->
data[k] * pow(f, -3.0) / deltaF);
95 gsl_matrix_set_identity(R);
107 gsl_matrix_set(R,
i, j, Rij);
108 gsl_matrix_set(R, j,
i, Rij);
112 gsl_linalg_cholesky_decomp(R);
117 double re = gsl_ran_gaussian_ziggurat(rng,
sigma);
118 double im = gsl_ran_gaussian_ziggurat(rng,
sigma);
120 htilde[
i]->
data->
data[k] += gsl_matrix_get(R,
i, j) * re;
121 htilde[
i]->
data->
data[k] += I * gsl_matrix_get(R,
i, j) * im;
132 # undef CLEANUP_AND_RETURN
156 size_t klow =
flow / deltaF;
162 for (k = 1; k < klow; ++k)
165 for (; k < length - 1; ++k)
168 OmegaGW->
data->
data[length - 1] = 0.0;
189 size_t klow =
flow / deltaF;
195 for (k = 1; k < klow; ++k)
198 for (; k < length - 1; ++k)
201 OmegaGW->
data->
data[length - 1] = 0.0;
286 # define CLEANUP_AND_RETURN(errnum) do { \
287 if (overlap) for (i = 0; i < numDetectors; ++i) XLALDestroyREAL8Sequence(overlap[i]); \
289 if (errnum) XLAL_ERROR(errnum); else return 0; \
303 if (h[
i]->
data->length != length
311 || (
size_t)floor(0.5 + 1.0/(
deltaT * OmegaGW->
deltaF)) != length)
321 }
else if (stride == length) {
343 for (j = 0; j < length - stride; ++j) {
344 double x = cos(
LAL_PI*j/(2.0 * (length - stride)));
345 double y = sin(
LAL_PI*j/(2.0 * (length - stride)));
356 # undef CLEANUP_AND_RETURN
444 deltaF = 1.0/(length * h[0]->
deltaT);
546 deltaF = 1.0/(length * h[0]->
deltaT);
585 if (
N == (
size_t)(-1))
589 deltaF = f[
N-1]/(length - 2);
591 for (
i = 0;
i <
N; ++
i)
592 Omega[
i] = log(Omega[
i]);
594 klow =
flow / deltaF;
601 for (k = 1; k < klow; ++k)
605 for (; k < length - 1; ++k) {
606 fk =
flow + k * deltaF;
607 while (f[
i] < fk &&
i <
N - 1)
609 x = (f[
i] - fk) / (f[
i] - f[
i-1]);
610 Omegak =
x * Omega[
i-1] + (1.0 -
x) * Omega[
i];
611 OmegaGW->
data->
data[k] = exp(Omegak);
614 OmegaGW->
data->
data[length - 1] = 0.0;
631 #include <lal/TimeSeries.h>
636 const double flow = 40.0;
638 const double srate = 16384.0;
646 size_t stride = length / 2;
651 rng = gsl_rng_alloc(gsl_rng_default);
658 for (j = 0; j < stride; ++j)
659 printf(
"%.9f\t%e\t%e\t%e\n",
XLALGPSGetREAL8(&seg[0]->
epoch) + j /
srate, seg[0]->
data->data[j], seg[1]->data->data[j], seg[2]->data->data[j]);
672 const double recdur = 1024.0;
673 const double segdur = 8.0;
674 const double srate = 4096.0;
675 const double flow = 1.0;
695 double deltaF = 1.0 /
segdur;
696 size_t reclen = recdur *
srate;
698 size_t stride = seglen / 2;
699 size_t numseg = 1 + (reclen - seglen)/stride;
700 size_t klow =
flow / deltaF;
704 if (reclen != (numseg - 1)*stride + seglen) {
705 fprintf(stderr,
"warning: numseg, seglen, reclen, and stride not commensurate\n");
706 reclen = (numseg - 1)*stride + seglen;
709 rng = gsl_rng_alloc(gsl_rng_default);
715 snprintf(seg[
i]->
name,
sizeof(seg[
i]->
name),
"%2s:%s", prefix, seg[
i]->
name);
720 for (
l = 0;
l < numseg; ++
l) {
725 memcpy(rec[
i]->
data->data +
l*stride, seg[
i]->data->data, (
l == numseg - 1 ? seglen : stride) *
sizeof(*rec[
i]->data->data));
728 fp = fopen(
"sgwb.dat",
"w");
729 for (j = 0; j < reclen; ++j) {
744 fp = fopen(
"sgwb-psd.dat",
"w");
745 for (k = klow; k < seglen/2; ++k) {
746 double f = k * deltaF;
748 fprintf(fp,
"%.9f\t%e", f, actual);
758 memset(psd[0]->
data->data, 0, psd[0]->data->length *
sizeof(*psd[0]->data->data));
759 for (
l = 0;
l < numseg; ++
l) {
762 for (j = 0; j < seglen; ++j) {
769 for (k = klow; k < seglen/2; ++k) {
770 double re = creal(htilde1->
data->
data[k]);
771 double im = cimag(htilde1->
data->
data[k]);
772 psd[0]->
data->
data[k] += fac * (re * re + im * im);
775 fp = fopen(
"sgwb-orf.dat",
"w");
776 for (k = klow; k < seglen/2; ++k) {
777 double f = k * deltaF;
778 double csdfac = 0.15 * pow(H0 / LAL_PI, 2.0) * pow(f, -3.0) *
Omega0 *
segdur;
static double sigma(const double a, const double b, const sysq *system)
Internal function that computes the spin-spin couplings.
size_t XLALSimReadDataFile2Col(double **xdat, double **ydat, LALFILE *fp)
Read a two-column data file.
LALFILE * XLALSimReadDataFileOpen(const char *fname)
Opens a specified data file, searching default path if necessary.
static int XLALSimSGWBSegment(REAL8TimeSeries **h, const LALDetector *detectors, size_t numDetectors, const REAL8FrequencySeries *OmegaGW, double H0, gsl_rng *rng)
#define CLEANUP_AND_RETURN(errnum)
int main(int argc, char *argv[])
const LALDetector lalCachedDetectors[LAL_NUM_DETECTORS]
int XLALFileClose(LALFILE *file)
REAL8FrequencySeries * XLALCreateREAL8FrequencySeries(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 XLALDestroyREAL8FrequencySeries(REAL8FrequencySeries *series)
void XLALDestroyCOMPLEX16FrequencySeries(COMPLEX16FrequencySeries *series)
int XLALSimSGWBPowerLawSpectrum(REAL8TimeSeries **h, const LALDetector *detectors, size_t numDetectors, size_t stride, double Omegaref, double alpha, double fref, double flow, double H0, gsl_rng *rng)
Routine that may be used to generate sequential segments of stochastic background gravitational wave ...
int XLALSimSGWBFlatSpectrum(REAL8TimeSeries **h, const LALDetector *detectors, size_t numDetectors, size_t stride, double Omega0, double flow, double H0, gsl_rng *rng)
Routine that may be used to generate sequential segments of stochastic background gravitational wave ...
REAL8FrequencySeries * XLALSimSGWBOmegaGWPowerLawSpectrum(double Omegaref, double alpha, double fref, double flow, double deltaF, size_t length)
Creates a frequency series that contains a power law SGWB spectrum with the specified power Omegaref ...
int XLALSimSGWB(REAL8TimeSeries **h, const LALDetector *detectors, size_t numDetectors, size_t stride, const REAL8FrequencySeries *OmegaGW, double H0, gsl_rng *rng)
Routine that may be used to generate sequential segments of stochastic background gravitational wave ...
REAL8FrequencySeries * XLALSimSGWBOmegaGWFlatSpectrum(double Omega0, double flow, double deltaF, size_t length)
Creates a frequency series that contains a flat SGWB spectrum with the specified power Omega0 above s...
REAL8FrequencySeries * XLALSimSGWBOmegaGWNumericalSpectrumFromFile(const char *fname, size_t length)
double XLALSimSGWBOverlapReductionFunction(double f, const LALDetector *detector1, const LALDetector *detector2)
Computes the overlap reduction function between two detectors at a specified frequency.
void XLALDestroyREAL8FFTPlan(REAL8FFTPlan *plan)
REAL8FFTPlan * XLALCreateReverseREAL8FFTPlan(UINT4 size, int measurelvl)
REAL8FFTPlan * XLALCreateForwardREAL8FFTPlan(UINT4 size, int measurelvl)
REAL8Sequence * XLALCreateREAL8Sequence(size_t length)
int XLALREAL8AverageSpectrumWelch(REAL8FrequencySeries *spectrum, const REAL8TimeSeries *tseries, UINT4 seglen, UINT4 stride, const REAL8Window *window, const REAL8FFTPlan *plan)
int XLALREAL8TimeFreqFFT(COMPLEX16FrequencySeries *freq, const REAL8TimeSeries *tser, const REAL8FFTPlan *plan)
int XLALREAL8FreqTimeFFT(REAL8TimeSeries *tser, const COMPLEX16FrequencySeries *freq, const REAL8FFTPlan *plan)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
const LALUnit lalStrainUnit
const LALUnit lalSecondUnit
const LALUnit lalDimensionlessUnit
LALUnit * XLALUnitMultiply(LALUnit *output, const LALUnit *unit1, const LALUnit *unit2)
void XLALDestroyREAL8Window(REAL8Window *window)
REAL8Window * XLALCreateHannREAL8Window(UINT4 length)
#define XLAL_ERROR_NULL(...)
void XLALAbortErrorHandler(const char *func, const char *file, int line, int errnum)
XLALErrorHandlerType * XLALSetErrorHandler(XLALErrorHandlerType *newHandler)
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
int XLALGPSCmp(const LIGOTimeGPS *t0, const LIGOTimeGPS *t1)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
LALDetector detectors[LAL_NUM_DETECTORS]