23#include <lal/XLALError.h>
24#include <lal/LALBarycenter.h>
25#include <lal/LALInitBarycenter.h>
26#include <lal/LogPrintf.h>
27#include <lal/CWMakeFakeData.h>
28#include <lal/LALConstants.h>
29#include <lal/ExtrapolatePulsarSpins.h>
30#include <lal/ComputeFstat.h>
31#include <lal/DetectorStates.h>
32#include <lal/LFTandTSutils.h>
33#include <lal/LALString.h>
34#include <lal/UserInput.h>
35#include <lal/LALPulsarVCSInfo.h>
77main(
int argc,
char *argv[] )
90 uvar->Delta[0] = -
LAL_PI / 2;
91 uvar->Delta[1] =
LAL_PI / 2;
94 uvar->f1dot[0] = -3
e-9;
97 uvar->f2dot[1] = 1
e-16;
98 uvar->FreqResolution[0] = 0.1;
99 uvar->FreqResolution[1] = 1;
100 uvar->numFreqBins[0] = 1000;
101 uvar->numFreqBins[1] = 100000;
102 uvar->Tseg[0] = 10 * 3600;
103 uvar->Tseg[1] = 250 * 3600;
104 uvar->orbitasini[0] = 0;
105 uvar->orbitasini[1] = 0;
106 uvar->orbitPeriod[0] = 0;
107 uvar->orbitPeriod[1] = 0;
108 uvar->orbitEcc[0] = 0;
109 uvar->orbitEcc[1] = 0;
110 uvar->orbitArgp[0] = 0;
111 uvar->orbitArgp[1] = 0;
112 uvar->orbitTp[0].gpsSeconds = 0;
113 uvar->orbitTp[1].gpsSeconds = 0;
115 uvar->numSegments = 90;
117 uvar->startTime.gpsSeconds = 711595934;
119 uvar->sharedWorkspace = 1;
121 uvar->perSegmentSFTs = 1;
130 uvar->outputInfo = NULL;
178 if ( have_orbitasini || have_orbitEcc || have_orbitPeriod || have_orbitArgp || have_orbitTp ) {
180 XLAL_CHECK( ( uvar->orbitasini[1] == 0 ) || ( have_orbitPeriod && ( uvar->orbitPeriod[0] > 0 ) && have_orbitTp ),
XLAL_EINVAL,
"If orbitasini>0 then we also need 'orbitPeriod>0' and 'orbitTp'\n" );
185 CHAR *logstring = NULL;
217 injectSqrtSX.sqrtSn[X] = 1;
223 optionalArgs.
Dterms = uvar->Dterms;
225 FILE *timingLogFILE = NULL;
226 FILE *timingParFILE = NULL;
227 if ( uvar->outputInfo != NULL ) {
228 char *parFname = NULL;
232 XLAL_CHECK_MAIN( ( timingLogFILE = fopen( uvar->outputInfo,
"ab" ) ) != NULL,
XLAL_ESYS,
"Failed to open '%s' for appending\n", uvar->outputInfo );
233 XLAL_CHECK_MAIN( ( timingParFILE = fopen( parFname,
"ab" ) ) != NULL,
XLAL_ESYS,
"Failed to open '%s' for appending\n", parFname );
235 fprintf( timingLogFILE,
"%s\n", logstring );
236 fprintf( timingParFILE,
"%s\n", logstring );
237 fprintf( timingParFILE,
"%%%%%8s %20s %20s %20s %20s %20s %20s %20s %20s %12s %20s %20s %20s %20s %20s\n",
238 "Nseg",
"Tseg",
"Freq",
"FreqBand",
"dFreq",
"f1dot",
"f2dot",
"Alpha",
"Delta",
"memUsageMB",
"asini",
"period",
"ecc",
"argp",
"tp" );
244#define drawFromREAL8Range(range) (range[0] + (range[1] - range[0]) * rand() / RAND_MAX )
245#define drawFromINT4Range(range) (range[0] + (INT4)round(1.0*(range[1] - range[0]) * rand() / RAND_MAX) )
246#define drawFromEPOCHRange(epoch, range) XLALGPSSetREAL8 ( epoch, (XLALGPSGetREAL8(&range[0]) + round(1.0*(XLALGPSGetREAL8(&range[1]) - XLALGPSGetREAL8(&range[0])) * rand() / RAND_MAX) ))
248 for (
INT4 i = 0;
i < uvar->numTrials;
i ++ ) {
250 Tseg_i = (
UINT4 ) uvar->Tsft * ceil( Tseg_i / uvar->Tsft );
253 for (
INT4 l = 0;
l < uvar->numSegments;
l ++ ) {
254 startTime_l->
data[
l] = (
l == 0 ) ? uvar->startTime : endTime_l->
data[
l - 1];
266 LIGOTimeGPS refTime = { uvar->startTime.gpsSeconds + 0.5 * uvar->numSegments * Tseg_i, 0 };
267 spinRange_i.refTime = refTime;
273 Doppler_i.refTime = refTime;
276 sDeltaRange[0] = sin( uvar->Delta[0] );
277 sDeltaRange[1] = sin( uvar->Delta[1] );
279 memcpy( &Doppler_i.fkdot, &spinRange_i.fkdot,
sizeof( Doppler_i.fkdot ) );;
289 REAL8 dFreq_i = FreqResolution_i / Tseg_i;
290 REAL8 FreqBand_i = numFreqBins_i * dFreq_i;
294 fprintf( stderr,
"trial %d/%d: Tseg = %.1f d, numSegments = %d, Alpha = %.2f rad, Delta = %.2f rad, Freq = %.6f Hz, f1dot = %.1e Hz/s, f2dot = %.1e Hz/s^2, R = %.2f, numFreqBins = %d, asini = %.2f, period = %.2f, ecc = %.2f, argp = %.2f, tp=%"LAL_GPS_FORMAT" [dFreq = %.2e Hz, FreqBand = %.2e Hz]\n",
295 i + 1, uvar->numTrials, Tseg_i / 86400.0, uvar->numSegments, Doppler_i.Alpha, Doppler_i.Delta, Doppler_i.fkdot[0], Doppler_i.fkdot[1], Doppler_i.fkdot[2], FreqResolution_i, numFreqBins_i, Doppler_i.asini, Doppler_i.period, Doppler_i.ecc, Doppler_i.argp,
LAL_GPS_PRINT( Doppler_i.tp ), dFreq_i, FreqBand_i );
297 spinRange_i.fkdotBand[0] = FreqBand_i;
298 REAL8 minCoverFreq_il, maxCoverFreq_il;
300 if ( ! uvar->perSegmentSFTs ) {
304 for (
INT4 l = 0;
l < uvar->numSegments;
l ++ ) {
305 if ( uvar->sharedWorkspace &&
l > 0 ) {
311 if ( uvar->perSegmentSFTs ) {
319 for (
INT4 l = 0;
l < uvar->numSegments;
l ++ ) {
325 for (
INT4 l = 0;
l < uvar->numSegments;
l ++ ) {
329 if ( timingLogFILE != NULL ) {
335 REAL8 memUsage = memEnd - memBase;
337 fprintf( stderr,
"%-15s: memoryUsage = %6.1f MB\n", FmethodName, memUsage );
339 if ( timingParFILE != NULL ) {
340 fprintf( timingParFILE,
"%10d %20d %20.16g %20.16g %20.16g %20.16g %20.16g %20.16g %20.16g %12g %20.16g %20.16g %20.16g %20.16g %"LAL_GPS_FORMAT"\n",
341 uvar->numSegments, Tseg_i, Doppler_i.fkdot[0], FreqBand_i, dFreq_i, Doppler_i.fkdot[1], Doppler_i.fkdot[2], Doppler_i.Alpha, Doppler_i.Delta, memUsage, Doppler_i.asini, Doppler_i.period, Doppler_i.ecc, Doppler_i.argp,
LAL_GPS_PRINT( Doppler_i.tp )
349 if ( timingLogFILE != NULL ) {
350 fclose( timingLogFILE );
352 if ( timingParFILE != NULL ) {
353 fclose( timingParFILE );
379 FILE *
file = fopen(
"/proc/self/status",
"r" );
382 if (
file == NULL ) {
386 if ( strncmp(
line,
"VmRSS:", 6 ) == 0 ) {
392 return ( result / 1024.0 );
398 int i = strlen(
line );
399 const char *
p =
line;
400 while ( *p < '0' || *p >
'9' ) {
int main(int argc, char *argv[])
REAL8 XLALGetCurrentHeapUsageMB(void)
static int parseLine(char *line)
#define drawFromEPOCHRange(epoch, range)
#define drawFromINT4Range(range)
#define drawFromREAL8Range(range)
void LALCheckMemoryLeaks(void)
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
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.
const UserChoices * XLALFstatMethodChoices(void)
Return pointer to a static array of all (available) FstatMethodType choices.
void XLALDestroyFstatInputVector(FstatInputVector *inputs)
Free all memory associated with a FstatInputVector structure.
int XLALAppendFstatTiming2File(const FstatInput *input, FILE *fp, BOOLEAN printHeader)
FstatInputVector * XLALCreateFstatInputVector(const UINT4 length)
Create a FstatInputVector of the given length, for example for setting up F-stat searches over severa...
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.
void XLALDestroyFstatResults(FstatResults *Fstats)
Free all memory associated with a FstatResults structure.
@ FMETHOD_RESAMP_BEST
Resamp: best guess of the fastest available implementation
@ FSTATQ_2F
Compute multi-detector .
@ FSTATQ_2F_PER_DET
Compute for each detector.
#define LAL_GPS_PRINT(gps)
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)
char char * XLALStringDuplicate(const char *s)
int char * XLALStringAppend(char *s, const char *append)
char * XLALVCSInfoString(const LALVCSInfoList vcs_list, const int verbose, const char *prefix)
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
MultiLIGOTimeGPSVector * XLALMakeMultiTimestamps(LIGOTimeGPS tStart, REAL8 Tspan, REAL8 Tsft, REAL8 Toverlap, UINT4 numDet)
Same as XLALMakeTimestamps() just for several detectors, additionally specify the number of detectors...
SFTCatalog * XLALMultiAddToFakeSFTCatalog(SFTCatalog *catalog, const LALStringVector *detectors, const MultiLIGOTimeGPSVector *timestamps)
Multi-detector and multi-timestamp wrapper of XLALAddToFakeSFTCatalog().
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
void XLALDestroyMultiTimestamps(MultiLIGOTimeGPSVector *multiTS)
Destroy a MultiLIGOTimeGPSVector timestamps vector.
void XLALDestroyTimestampVector(LIGOTimeGPSVector *vect)
De-allocate a LIGOTimeGPSVector.
LALStringVector * XLALCreateStringVector(const CHAR *str1,...)
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_MAIN(assertion,...)
p
RUN ANALYSIS SCRIPT ###.
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.
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 ...
BOOLEAN collectTiming
a flag to turn on/off the collection of F-stat-method-specific timing-data
UINT4 Dterms
Number of Dirichlet kernel terms, used by some Demod methods; see FstatMethodType.
FstatMethodType FstatMethod
Method to use for computing the -statistic.
XLALComputeFstat() computed results structure.
A vector of 'timestamps' of type LIGOTimeGPS.
LIGOTimeGPS * data
array of timestamps
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...
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
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()