20#include <lal/XLALError.h>
21#include <lal/LALBarycenter.h>
22#include <lal/LALInitBarycenter.h>
23#include <lal/UniversalDopplerMetric.h>
24#include <lal/LogPrintf.h>
25#include <lal/CWMakeFakeData.h>
26#include <lal/LALConstants.h>
27#include <lal/ExtrapolatePulsarSpins.h>
28#include <lal/ComputeFstat.h>
29#include <lal/DetectorStates.h>
30#include <lal/LFTandTSutils.h>
31#include <lal/LALString.h>
40main(
int argc,
char *argv[] )
60 injectSqrtSX.sqrtSn[X] = 0;
61 assumeSqrtSX.sqrtSn[X] = 1.0 + 2.0 * X;
111 Doppler.Delta = -0.5;
112 Doppler.fkdot[0] = Freq;
113 Doppler.fkdot[1] = -1
e-9;
114 Doppler.refTime = refTime;
116 Doppler.asini = asini;
119 Doppler.period = Period;
129 REAL8 dPeriod = 3600;
130 UINT4 numFreqBins = 1000;
132 UINT4 numf1dotPoints = 2;
133 UINT4 numSkyPoints = 2;
134 UINT4 numPeriodPoints = 2;
137 spinRange.refTime = refTime;
138 memcpy( &spinRange.fkdot, &injectSources->
data[0].
Doppler.
fkdot,
sizeof( spinRange.fkdot ) );
139 spinRange.fkdotBand[0] = ( numFreqBins - 1 ) * dFreq - 10 *
LAL_REAL8_EPS;
140 spinRange.fkdotBand[1] = ( numf1dotPoints - 1 ) * df1dot - 10 *
LAL_REAL8_EPS;
142 Doppler.fkdot[0] -= 0.4 * spinRange.fkdotBand[0];
144 REAL8 minCoverFreq, maxCoverFreq;
154 FstatInput *input_seg1[FMETHOD_END], *input_seg2[FMETHOD_END];
155 FstatResults *results_seg1[FMETHOD_END], *results_seg2[FMETHOD_END];
156 for (
UINT4 iMethod = FMETHOD_START; iMethod < FMETHOD_END; iMethod ++ ) {
157 results_seg1[iMethod] = results_seg2[iMethod] = NULL;
165 optionalArgs.
prevInput = input_seg1[iMethod];
174 for (
UINT4 iSky = 0; iSky < numSkyPoints; iSky ++ ) {
175 for (
UINT4 if1dot = 0; if1dot < numf1dotPoints; if1dot ++ ) {
176 for (
UINT4 iPeriod = 0; iPeriod < numPeriodPoints; iPeriod ++ ) {
182 for (
UINT4 iMethod = FMETHOD_START; iMethod < FMETHOD_END; iMethod ++ ) {
186 if ( firstMethod == FMETHOD_START ) {
187 firstMethod = iMethod;
201 if ( iMethod != firstMethod ) {
214 if ( first_SRC_a == NULL ) {
219 for (
UINT4 X = 0; X < first_SRC_a->
length; ++X ) {
239 for (
UINT4 X = 0; X < first_SRC_a->
length; ++X ) {
250 VectorComparison tol = { .relErr_L1 = 5
e-5, .relErr_L2 = 5
e-5, .angleV = 5
e-5, .relErr_atMaxAbsx = 5
e-5, .relErr_atMaxAbsy = 5
e-5 };
261 Doppler.period += dPeriod;
265 Doppler.fkdot[1] += df1dot;
269 Doppler.Alpha += dSky;
290 FstatResults *results_input_slice = NULL, *results_sft_slice = NULL;
291 FstatInput *input_sft_slice = NULL, *input_slice = NULL;
305 printFstatResults2File( results_input_slice,
"FstatInputTimeslice", numSkyPoints, numf1dotPoints, numPeriodPoints, whatToCompute );
306 printFstatResults2File( results_sft_slice,
"SFTCatalogTimeslice", numSkyPoints, numf1dotPoints, numPeriodPoints, whatToCompute );
309 XLALPrintInfo(
"Comparing results between SFTCatalogTimeslice and FstatInputTimeslice\n" );
311 XLALPrintError(
"Comparison between SFTCatalogTimeslice and FstatInputTimeslice failed\n" );
316 for (
UINT4 iMethod = FMETHOD_START; iMethod < FMETHOD_END; iMethod ++ ) {
356 tol.relErr_L1 = 2.5e-2;
357 tol.relErr_L2 = 2.2e-2;
359 tol.relErr_atMaxAbsx = 2.1e-2;
360 tol.relErr_atMaxAbsy = 2.1e-2;
369 v1.length = numFreqBins;
370 v2.length = numFreqBins;
371 v1.data = result1->
twoF;
372 v2.data = result2->
twoF;
389 c1.length = numFreqBins;
390 c2.length = numFreqBins;
392 c1.data = result1->
Fa;
393 c2.data = result2->
Fa;
397 c1.data = result1->
Fb;
398 c2.data = result2->
Fb;
413 snprintf( fname,
sizeof( fname ) - 1,
"twoF%s-iSky%02d-if1dot%02d-iPeriod%02d.dat", MethodName, iSky, if1dot, iPeriod );
419 fprintf(
fp,
"%20.16g %10.4g %10.4g %10.4g %10.4g %10.4g\n",
420 Freq_k, results->
twoF[
k],
421 crealf( results->
Fa[
k] ), cimagf( results->
Fa[
k] ),
422 crealf( results->
Fb[
k] ), cimagf( results->
Fb[
k] )
426 Freq_k, results->
twoF[
k] );
int main(int argc, char *argv[])
static int compareFstatResults(const FstatResults *result1, const FstatResults *result2)
static int printFstatResults2File(const FstatResults *results, const char *MethodName, const INT4 iSky, const INT4 if1dot, const INT4 iPeriod, FstatQuantities whatToCompute)
void LALCheckMemoryLeaks(void)
PulsarParamsVector * XLALCreatePulsarParamsVector(UINT4 numPulsars)
Create zero-initialized PulsarParamsVector for numPulsars.
void XLALDestroyPulsarParamsVector(PulsarParamsVector *ppvect)
Destructor for PulsarParamsVector type.
FstatMethodType
Different methods available to compute the -statistic, falling into two broad classes:
const CHAR * XLALGetFstatInputMethodName(const FstatInput *input)
Returns the human-readable name of the -statistic method being used by a FstatInput structure.
const FstatOptionalArgs FstatOptionalArgsDefaults
Global initializer for setting FstatOptionalArgs to default values.
int XLALFstatInputTimeslice(FstatInput **slice, const FstatInput *input, const LIGOTimeGPS *minStartGPS, const LIGOTimeGPS *maxStartGPS)
Create and return an FstatInput 'timeslice' for given input FstatInput object [must be using LALDemod...
int XLALExtractResampledTimeseries(MultiCOMPLEX8TimeSeries **multiTimeSeries_SRC_a, MultiCOMPLEX8TimeSeries **multiTimeSeries_SRC_b, FstatInput *input)
Extracts the resampled timeseries from a given FstatInput structure (must have been initialized previ...
void XLALDestroyFstatInput(FstatInput *input)
Free all memory associated with a FstatInput structure.
FstatQuantities
Bit-field of -statistic quantities which can be computed by XLALComputeFstat().
FstatInput * XLALCreateFstatInput(const SFTCatalog *SFTcatalog, const REAL8 minCoverFreq, const REAL8 maxCoverFreq, const REAL8 dFreq, const EphemerisData *ephemerides, const FstatOptionalArgs *optionalArgs)
Create a fully-setup FstatInput structure for computing the -statistic using XLALComputeFstat().
int XLALComputeFstat(FstatResults **Fstats, FstatInput *input, const PulsarDopplerParams *doppler, const UINT4 numFreqBins, const FstatQuantities whatToCompute)
Compute the -statistic over a band of frequencies.
int XLALFstatMethodIsAvailable(FstatMethodType method)
Return true if given FstatMethodType corresponds to a valid and available Fstat method,...
void XLALDestroyFstatResults(FstatResults *Fstats)
Free all memory associated with a FstatResults structure.
@ FMETHOD_RESAMP_GENERIC
Resamp: generic implementation
@ FMETHOD_DEMOD_BEST
Demod: best guess of the fastest available hotloop
@ FMETHOD_RESAMP_BEST
Resamp: best guess of the fastest available implementation
@ FSTATQ_2F
Compute multi-detector .
@ FSTATQ_FAFB
Compute multi-detector and .
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
#define XLAL_INIT_DECL(var,...)
void * XLALCalloc(size_t m, size_t n)
int XLALCompareREAL4Vectors(VectorComparison *result, const REAL4Vector *x, const REAL4Vector *y, const VectorComparison *tol)
Compare two REAL4 vectors using various different comparison metrics.
int XLALCompareCOMPLEX8Vectors(VectorComparison *result, const COMPLEX8Vector *x, const COMPLEX8Vector *y, const VectorComparison *tol)
Compare two COMPLEX8 vectors using various different comparison metrics.
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
LIGOTimeGPSVector * XLALMakeTimestamps(LIGOTimeGPS tStart, REAL8 Tspan, REAL8 Tsft, REAL8 Toverlap)
Given a start-time, Tspan, Tsft and Toverlap, returns a list of timestamps covering this time-stretch...
SFTCatalog * XLALMultiAddToFakeSFTCatalog(SFTCatalog *catalog, const LALStringVector *detectors, const MultiLIGOTimeGPSVector *timestamps)
Multi-detector and multi-timestamp wrapper of XLALAddToFakeSFTCatalog().
void XLALDestroyMultiTimestamps(MultiLIGOTimeGPSVector *multiTS)
Destroy a MultiLIGOTimeGPSVector timestamps vector.
int XLALSFTCatalogTimeslice(SFTCatalog *slice, const SFTCatalog *catalog, const LIGOTimeGPS *minStartGPS, const LIGOTimeGPS *maxStartGPS)
Set a SFT catalog 'slice' to a timeslice of a larger SFT catalog 'catalog', with entries restricted t...
void XLALDestroyStringVector(LALStringVector *vect)
LALStringVector * XLALCreateStringVector(const CHAR *str1,...)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
int XLALGPSCmp(const LIGOTimeGPS *t0, const LIGOTimeGPS *t1)
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
Struct of optional 'advanced level' and (potentially method-specific) arguments to be passed to the ...
MultiNoiseFloor * injectSqrtSX
Single-sided PSD values for fake Gaussian noise to be added to SFT data.
PulsarParamsVector * injectSources
Vector of parameters of CW signals to simulate and inject.
MultiNoiseFloor * assumeSqrtSX
Single-sided PSD values to be used for computing SFT noise weights instead of from a running median o...
BOOLEAN resampFFTPowerOf2
Resamp: round up FFT lengths to next power of 2; see FstatMethodType.
FstatInput * prevInput
An FstatInput structure from a previous call to XLALCreateFstatInput(); may contain common workspace ...
FstatMethodType FstatMethod
Method to use for computing the -statistic.
XLALComputeFstat() computed results structure.
UINT4 numDetectors
Number of detectors over which the were computed.
FstatQuantities whatWasComputed
Bit-field of which -statistic quantities were computed.
REAL8 dFreq
Spacing in frequency between each computed -statistic.
UINT4 numFreqBins
Number of frequencies at which the were computed.
REAL4 * twoF
If whatWasComputed & FSTATQ_2F is true, the multi-detector values computed at numFreqBins frequencie...
COMPLEX8 * Fa
If whatWasComputed & FSTATQ_PARTS is true, the multi-detector and computed at numFreqBins frequenci...
PulsarDopplerParams doppler
Doppler parameters, including the starting frequency, at which the were computed.
Multi-IFO container for COMPLEX8 resampled timeseries.
COMPLEX8TimeSeries ** data
array of COMPLEX8TimeSeries (pointers)
UINT4 length
number of IFOs
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
array of detector-specific 'noise floors' (ie PSD values), assumed constant over the frequency-band o...
REAL8 aCross
Signal amplitude (cross)
REAL8 psi
polarization angle psi
REAL8 aPlus
Signal amplitude (plus)
REAL8 phi0
initial signal-phase (at some reference time)
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
PulsarSpins fkdot
Intrinsic spins: [Freq, f1dot, f2dot, ... ] where fkdot = d^kFreq/dt^k.
PulsarAmplitudeParams Amp
'Amplitude-parameters': h0, cosi, phi0, psi
PulsarDopplerParams Doppler
'Phase-evolution parameters': {skypos, fkdot, orbital params }
Straightforward vector type of N PulsarParams structs.
PulsarParams * data
array of pulsar-signal parameters
Contains a "spin-range", ie spins and corresponding bands at a given (SSB) reference GPS-time .
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Struct holding the results of comparing two floating-point vectors (real-valued or complex),...