Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
CWSignalBandTest.c
Go to the documentation of this file.
1/*
2 *
3 * Copyright (C) 2018 Reinhard Prix
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with with program; see the file COPYING. If not, write to the
17 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18 * MA 02110-1301 USA
19 */
20
21/*********************************************************************************/
22/**
23 * \author Reinhard Prix
24 * \file
25 * \brief Test for XLALCWSignal[Covering]Band().
26 *
27 */
28#include <math.h>
29
30#include <lal/Date.h>
31#include <lal/AVFactories.h>
32#include <lal/PulsarDataTypes.h>
33#include <lal/ExtrapolatePulsarSpins.h>
34#include <lal/LALInitBarycenter.h>
35#include <lal/LogPrintf.h>
36
38
39int main( void )
40{
41
42 UINT4 seed = 1;
43 srand( seed );
44 UINT4 Ntrials = 100;
45
46 // test functions predicting CW signal bands
47 // 1) XLALCWSignalCoveringBand(): total covering band for all signals over the sky and for given spin ranges
48 // 2a) XLALSingleCWSignalCoveringBand(): actual covering band for a specific signal
49 // 2b) XLALPrepareSingleCWSignalCoveringBand(): helper function compute detStates + 'worst' skyposition for widest Doppler band
50 LIGOTimeGPS t0 = { 1206527521, 0};
51 REAL8 Tspan = 1e6;
53 XLALGPSAdd( &t1, Tspan );
54 LIGOTimeGPS refTime = t0;
55 XLALGPSAdd( &t1, -Tspan ); // reftime outside of search-band for generality
56
57 REAL8 asini = 0;
58 REAL8 period = 0;
59 REAL8 ecc = 0;
60 LIGOTimeGPS tp = {0, 0};
61 REAL8 argp = 0;
62
63 PulsarSpinRange spinRange;
65 fkdot[0] = 1000;
66 fkdot[1] = -3e-9;
67 fkdot[2] = 0;
68 XLAL_CHECK_MAIN( XLALInitPulsarSpinRangeFromSpins( &spinRange, &refTime, fkdot, fkdot ) == XLAL_SUCCESS, XLAL_EFUNC );
69
70 FILE *fp = stderr;
71 if ( lalDebugLevel & LALINFOBIT ) {
72 XLAL_CHECK_MAIN( ( fp = fopen( "test_SignalCoveringBand.dat", "wb" ) ) != NULL, XLAL_EFUNC );
73 }
74
75 // ----- 1) get covering band for all signals over the sky for given spin-range
76 REAL8 minCoverFreq, maxCoverFreq;
77 REAL8 tic = XLALGetCPUTime();
78 XLAL_CHECK_MAIN( XLALCWSignalCoveringBand( &minCoverFreq, &maxCoverFreq, &t0, &t1, &spinRange, asini, period, ecc ) == XLAL_SUCCESS, XLAL_EFUNC );
79 REAL8 time_XLALCWSignalCoveringBand = XLALGetCPUTime() - tic;
80 REAL8 coveringFreqBand = maxCoverFreq - minCoverFreq;
81 fprintf( fp, "coveringFreqBand = %g;\n", coveringFreqBand );
82
83 // ----- 2a+b) get covering band for specific signals drawn randomly over the sky, check they satisfy 'width<=worst' and lie within total coverage band
84 DetectorStateSeries *detStates;
85 SkyPosition skypos_maxdoppler;
86 REAL8 dT = 1800;
87 const LALDetector *detector;
88 XLAL_CHECK_MAIN( ( detector = XLALGetSiteInfo( "H1" ) ) != NULL, XLAL_EFUNC );
90 XLAL_CHECK_MAIN( ( edat = XLALInitBarycenter( TEST_PKG_DATA_DIR "earth00-40-DE405.dat.gz", TEST_PKG_DATA_DIR "sun00-40-DE405.dat.gz" ) ) != NULL, XLAL_EFUNC );
91 tic = XLALGetCPUTime();
92 XLAL_CHECK_MAIN( ( detStates = XLALPrepareCWSignalBand( &skypos_maxdoppler, t0, Tspan, dT, detector, edat ) ) != NULL, XLAL_EFUNC );
93 REAL8 time_XLALPrepareCWSignalBand = XLALGetCPUTime() - tic;
94 REAL8 minFreq, maxFreq;
96 doppler.refTime = refTime;
97 doppler.Alpha = skypos_maxdoppler.longitude;
98 doppler.Delta = skypos_maxdoppler.latitude;
99 memcpy( doppler.fkdot, fkdot, sizeof( fkdot ) );
100 doppler.asini = asini;
101 doppler.period = period;
102 doppler.ecc = ecc;
103 doppler.tp = tp;
104 doppler.argp = argp;
105
106 XLAL_CHECK_MAIN( XLALCWSignalBand( &minFreq, &maxFreq, detStates, &doppler ) == XLAL_SUCCESS, XLAL_EFUNC );
107 REAL8 maxFreqBandWidth = maxFreq - minFreq;
108
109 fprintf( fp, "Alpha_max = %g; Delta_max = %g; minFreq = %g; maxFreq = %g;\nmaxFreqBandWidth = %g;\n",
110 doppler.Alpha, doppler.Delta, minFreq, maxFreq, maxFreqBandWidth );
111 XLAL_CHECK_MAIN( maxFreqBandWidth <= coveringFreqBand, XLAL_EFAILED,
112 "(maxFreqBandWidth = %g) > (coveringFreqBand = %g)\n",
113 maxFreqBandWidth, coveringFreqBand );
114
115 fprintf( fp, "%%%% Alpha Delta minFreq maxFreq FreqBandWidth\n" );
116 fprintf( fp, "results = [\n" );
117 REAL8 time_XLALCWSignalBand = 0;
118 for ( UINT4 i = 0; i < Ntrials; i ++ ) {
119 // draw random skyposition isotropically over the sky
120 doppler.Alpha = LAL_TWOPI * ( 1.0 * rand() / ( RAND_MAX + 1.0 ) ); // alpha uniform in [0, 2pi)
121 doppler.Delta = LAL_PI_2 - acos( 1 - 2.0 * rand() / RAND_MAX ); // sin(delta) uniform in [-1,1]
122 tic = XLALGetCPUTime();
123 XLAL_CHECK_MAIN( XLALCWSignalBand( &minFreq, &maxFreq, detStates, &doppler ) == XLAL_SUCCESS, XLAL_EFUNC );
124 time_XLALCWSignalBand += ( XLALGetCPUTime() - tic );
125 REAL8 FreqBandWidth = maxFreq - minFreq;
126 XLAL_CHECK_MAIN( FreqBandWidth <= maxFreqBandWidth, XLAL_EFAILED,
127 "Alpha = %g, Delta = %g: (FreqBandWidth = %g) > (maxFreqBandWidth = %g)\n",
128 doppler.Alpha, doppler.Delta, FreqBandWidth, maxFreqBandWidth );
129 fprintf( fp, "%10g, %10g %10.8g %10.8g %10.2e;\n", doppler.Alpha, doppler.Delta, minFreq, maxFreq, FreqBandWidth );
130 }
131 fprintf( fp, "];\n" );
132 time_XLALCWSignalBand /= Ntrials;
133
134 fprintf( fp, "time_XLALCWSignalCoveringBand = %.2e;\n", time_XLALCWSignalCoveringBand );
135 fprintf( fp, "time_XLALPrepareCWSignalBand = %.2e;\n", time_XLALPrepareCWSignalBand );
136 fprintf( fp, "time_XLALCWSignalBand = %.2e;\n", time_XLALCWSignalBand );
137
140
141 if ( fp ) {
142 fclose( fp );
143 }
144
146
147 return XLAL_SUCCESS;
148
149} /* main() */
int main(void)
int test_SignalCoveringBand(UINT4 Ntrials)
void LALCheckMemoryLeaks(void)
#define fprintf
double e
void XLALDestroyDetectorStateSeries(DetectorStateSeries *detStates)
Get rid of a DetectorStateSeries.
DetectorStateSeries * XLALPrepareCWSignalBand(SkyPosition *skypos_maxdoppler, const LIGOTimeGPS tStart, const REAL8 Tspan, const REAL8 dT, const LALDetector *detector, const EphemerisData *edat)
(Optional) Helper function for using XLALCWSignalBand(): compute DetectorStateSeries for given time-s...
int XLALInitPulsarSpinRangeFromSpins(PulsarSpinRange *range, const LIGOTimeGPS *refTime, const PulsarSpins fkdot1, const PulsarSpins fkdot2)
Initialise a PulsarSpinRange struct from two PulsarSpins structs.
int XLALCWSignalBand(REAL8 *minFreq, REAL8 *maxFreq, const DetectorStateSeries *detStates, const PulsarDopplerParams *doppler)
Determines the frequency band occupied by the frequency evolution of a given CW signal between two GP...
int XLALCWSignalCoveringBand(REAL8 *minCoverFreq, REAL8 *maxCoverFreq, const LIGOTimeGPS *time1, const LIGOTimeGPS *time2, const PulsarSpinRange *spinRange, const REAL8 binaryMaxAsini, const REAL8 binaryMinPeriod, const REAL8 binaryMaxEcc)
Determines a frequency band which covers the frequency evolution of a band of CW signals between two ...
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 LAL_PI_2
#define LAL_TWOPI
double REAL8
#define XLAL_INIT_DECL(var,...)
uint32_t UINT4
#define lalDebugLevel
LALINFOBIT
REAL8 XLALGetCPUTime(void)
REAL8 PulsarSpins[PULSAR_MAX_SPINS]
Typedef for fixed-size array holding GW frequency and derivatives fk = d^k Freq/dt^k|(tau_ref)
const LALDetector * XLALGetSiteInfo(const CHAR *name)
Find the site geometry-information 'LALDetector' for given a detector name (or prefix).
Definition: SFTnaming.c:218
#define XLAL_CHECK_MAIN(assertion,...)
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFAILED
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
double t0
Timeseries of DetectorState's, representing the detector-info at different timestamps.
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
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 .
REAL8 longitude
REAL8 latitude