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
pulsar_frequency_evolution.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2009, 2016 Matt Pitkin
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20/* Calculate the frequency evolution of a pulsar over time */
21
22#include "config.h"
23#include <stdlib.h>
24#include <stdio.h>
25#include <math.h>
26
27#include <lal/LALStdlib.h>
28#include <lal/LALAtomicDatatypes.h>
29#include <lal/LALDatatypes.h>
30#include <lal/AVFactories.h>
31#include <lal/LALgetopt.h>
32#include <lal/SkyCoordinates.h>
33#include <lal/DetectorSite.h>
34#include <lal/DetResponse.h>
35#include <lal/ReadPulsarParFile.h>
36#include <lal/BinaryPulsarTiming.h>
37#include <lal/LALString.h>
38#include <lal/SFTfileIO.h>
39#include <lal/LALBarycenter.h>
40#include <lal/LALInitBarycenter.h>
41#include <lal/LALConstants.h>
42#include <lal/XLALError.h>
43
44#include <getopt.h>
45
46#include "ppe_models.h"
47#include "ppe_utils.h"
48
49#define CODEUSAGE \
50"Usage: %s [options]\n\n"\
51" --help display this message\n"\
52" --detector IFOs to use e.g. H1\n"\
53" --par-file pulsar parameter (.par) file (full path)\n"\
54" --start start time GPS\n"\
55" --timespan time span to calculate over (seconds)\n"\
56" --deltat length of each step (seconds)\n"\
57" --output-dir output directory\n"\
58" --harmonic the frequency harmonic to output [default: 2]\n"\
59"\n"
60
61typedef struct tagInputParams {
62 CHAR *det;
64
69
72
73void get_input_args( InputParams *inputParams, INT4 argc, CHAR *argv[] );
74
77 REAL8Vector *dfsolar, REAL8Vector *dfbinary, REAL8Vector *dftotal );
78
79int main( int argc, char *argv[] )
80{
81 InputParams inputs;
82
83 FILE *fp = NULL;
84
86
87 UINT4 i = 0, npoints = 0;
88
89 EphemerisData *edat = NULL;
91 LALDetector det;
92 BarycenterInput baryinput;
93
94 CHAR outputfile[256];
95
96 get_input_args( &inputs, argc, argv );
97
98 npoints = floor( ( REAL8 )inputs.timespan / ( REAL8 )inputs.deltat );
99
101
102 /* set up ephemerises */
103 det = *XLALGetSiteInfo( inputs.det ); /* just set site as LHO */
104 baryinput.site = det;
105 baryinput.site.location[0] /= LAL_C_SI;
106 baryinput.site.location[1] /= LAL_C_SI;
107 baryinput.site.location[2] /= LAL_C_SI;
108
109 CHAR *earthfile = NULL, *sunfile = NULL, *timefile = NULL;
110 TimeCorrectionType ttype;
111
112 ttype = XLALAutoSetEphemerisFiles( &earthfile, &sunfile, &timefile, params, inputs.start, inputs.start + inputs.timespan );
113
114 edat = XLALInitBarycenter( earthfile, sunfile );
116
117 /* vector to hold frequency time series and Doppler shift time series' */
118 REAL8Vector *freqs = NULL, *dopplerss = NULL, *dopplerbs = NULL, *doppler = NULL;
119 freqs = XLALCreateREAL8Vector( npoints );
120 dopplerss = XLALCreateREAL8Vector( npoints );
121 dopplerbs = XLALCreateREAL8Vector( npoints );
122 doppler = XLALCreateREAL8Vector( npoints );
123
124 /* calculate the frequency every minute for the mean values */
125 get_freq( inputs.start, ( REAL8 )inputs.deltat, inputs.freqharm, params, baryinput, edat, tdat, ttype, freqs, dopplerss, dopplerbs, doppler );
126
127 sprintf( outputfile, "%s/frequency_evolution_%s.txt", inputs.outputdir, inputs.det );
128
129 fp = fopen( outputfile, "w" );
130 for ( i = 0; i < freqs->length; i++ ) {
131 fprintf( fp, "%lf\t%.15lf\t%.10e\t%.10e\t%.10e\n", inputs.start + ( double )( i * inputs.deltat ), freqs->data[i],
132 dopplerss->data[i], dopplerbs->data[i], doppler->data[i] );
133 }
134 fclose( fp );
135
136 XLALDestroyREAL8Vector( freqs );
137 XLALDestroyREAL8Vector( dopplerss );
138 XLALDestroyREAL8Vector( dopplerbs );
139 XLALDestroyREAL8Vector( doppler );
140
142
143 return 0;
144}
145
146void get_input_args( InputParams *inputParams, INT4 argc, CHAR *argv[] )
147{
148 struct LALoption long_options[] = {
149 { "help", no_argument, 0, 'h' },
150 { "detector", required_argument, 0, 'D' },
151 { "par-file", required_argument, 0, 'P' },
152 { "start", required_argument, 0, 's' },
153 { "timespan", required_argument, 0, 't' },
154 { "deltat", required_argument, 0, 'd' },
155 { "output-dir", required_argument, 0, 'o' },
156 { "harmonic", required_argument, 0, 'f' },
157 { 0, 0, 0, 0 }
158 };
159
160 CHAR args[] = "hD:P:s:t:d:o:f:";
161 CHAR *program = argv[0];
162
163 /* set default frequency harmonic to 2 */
164 inputParams->freqharm = 2.;
165
166 while ( 1 ) {
167 INT4 option_index = 0;
168 INT4 c;
169
170 c = LALgetopt_long( argc, argv, args, long_options, &option_index );
171 if ( c == -1 ) { /* end of options */
172 break;
173 }
174
175 switch ( c ) {
176 case 0:
177 if ( long_options[option_index].flag ) {
178 break;
179 } else {
180 XLALPrintError( "Error passing option %s with argument %s\n", long_options[option_index].name, LALoptarg );
182 }
183 case 'h': /* help message */
184 fprintf( stderr, CODEUSAGE, program );
185 exit( 0 );
186 case 'D':
187 inputParams->det = XLALStringDuplicate( LALoptarg );
188 break;
189 case 'P':
190 inputParams->parfile = XLALStringDuplicate( LALoptarg );
191 break;
192 case 's':
193 inputParams->start = atoi( LALoptarg );
194 break;
195 case 't':
196 inputParams->timespan = atoi( LALoptarg );
197 break;
198 case 'd':
199 inputParams->deltat = atoi( LALoptarg );
200 break;
201 case 'o':
202 inputParams->outputdir = XLALStringDuplicate( LALoptarg );
203 break;
204 case 'f':
205 inputParams->freqharm = atof( LALoptarg );
206 break;
207 case '?':
208 XLALPrintError( "Unknown error while parsing options." );
210 default:
211 XLALPrintError( "Unknown error while parsing options." );
213 }
214 }
215
216 if ( inputParams->start < 0. ) {
217 XLALPrintError( "Input start time must be positive" );
219 }
220
221 if ( inputParams->timespan < 0. ) {
222 XLALPrintError( "Input timespan must be positive" );
224 }
225
226 if ( inputParams->deltat < 0. ) {
227 XLALPrintError( "Input time steps must be positive" );
229 }
230
231 if ( inputParams->freqharm < 0. ) {
232 XLALPrintError( "Input frequency harmonic must be positive" );
234 }
235}
236
237/* function to return a vector of the pulsar frequency for each data point */
238void get_freq( REAL8 start, REAL8 deltaT, REAL8 freqharm,
241 REAL8Vector *freqs, REAL8Vector *dfsolar, REAL8Vector *dfbinary,
242 REAL8Vector *dftotal )
243{
244 UINT4 i = 0;
245
246 REAL8 T0 = 0., DT = 0.;
247 REAL8 time0 = 0., timep1 = 0.;
248
249 EarthState earth, earth2;
250 EmissionTime emit, emit2;
251 BinaryPulsarInput binput;
252 BinaryPulsarOutput boutput, boutput2;
253
254 /* if edat is NULL then show error */
255 if ( edat == NULL ) {
256 XLALPrintError( "%s: Ephemeris does not appear to be initialised!", __func__ );
258 }
259
260 /* get right ascension and declination */
261 REAL8 ra = 0.;
262 if ( PulsarCheckParam( params, "RA" ) ) {
263 ra = PulsarGetREAL8Param( params, "RA" );
264 } else if ( PulsarCheckParam( params, "RAJ" ) ) {
265 ra = PulsarGetREAL8Param( params, "RAJ" );
266 } else {
267 XLALPrintError( "%s: No source right ascension specified!", __func__ );
269 }
270 REAL8 dec = 0.;
271 if ( PulsarCheckParam( params, "DEC" ) ) {
272 dec = PulsarGetREAL8Param( params, "DEC" );
273 } else if ( PulsarCheckParam( params, "DECJ" ) ) {
274 dec = PulsarGetREAL8Param( params, "DECJ" );
275 } else {
276 XLALPrintError( "%s: No source declination specified!", __func__ );
278 }
279
280 REAL8 pmra = PulsarGetREAL8ParamOrZero( params, "PMRA" );
281 REAL8 pmdec = PulsarGetREAL8ParamOrZero( params, "PMDEC" );
282 REAL8 pepoch = PulsarGetREAL8ParamOrZero( params, "PEPOCH" );
283 REAL8 posepoch = PulsarGetREAL8ParamOrZero( params, "POSEPOCH" );
284 REAL8 px = PulsarGetREAL8ParamOrZero( params, "PX" ); /* parallax */
285 REAL8 dist = PulsarGetREAL8ParamOrZero( params, "DIST" ); /* distance */
286
287 /* set 1/distance if parallax or distance value is given (1/sec) */
288 if ( px != 0. ) {
289 bary.dInv = px * 1e-3 * LAL_C_SI / LAL_PC_SI;
290 } else if ( dist != 0. ) {
291 bary.dInv = LAL_C_SI / ( dist * 1e3 * LAL_PC_SI );
292 } else {
293 bary.dInv = 0.;
294 }
295
296 /* set the position and frequency epochs if not already set */
297 if ( pepoch == 0. && posepoch != 0. ) {
298 pepoch = posepoch;
299 } else if ( posepoch == 0. && pepoch != 0. ) {
300 posepoch = pepoch;
301 }
302
303 /* get frequencies */
304 REAL8 taylorcoeff = 1., tmpdt = 0.;
305 const REAL8Vector *fs = PulsarGetREAL8VectorParam( params, "F" );
306
307 for ( i = 0; i < freqs->length; i++ ) {
308 LIGOTimeGPS tgps;
309
310 T0 = pepoch;
311 time0 = start + deltaT * ( double )i;
312 timep1 = time0 + 1.;
313
314 DT = time0 - T0;
315
316 bary.tgps.gpsSeconds = ( UINT8 )floor( time0 );
317 bary.tgps.gpsNanoSeconds = ( UINT8 )floor( ( fmod( time0, 1. ) * 1e9 ) );
318
319 bary.delta = dec + DT * pmdec;
320 bary.alpha = ra + DT * pmra / cos( bary.delta );
321
322 /* call barycentring routines */
323 XLALBarycenterEarthNew( &earth, &bary.tgps, edat, tdat, ttype );
324 XLALBarycenter( &emit, &bary, &earth );
325
326 /* add 1 sec so we can get doppler shift */
327 bary.tgps.gpsSeconds = ( UINT8 )floor( timep1 );
328 bary.tgps.gpsNanoSeconds = ( UINT8 )floor( ( fmod( timep1, 1. ) * 1e9 ) );
329
330 XLALBarycenterEarthNew( &earth2, &bary.tgps, edat, tdat, ttype );
331 XLALBarycenter( &emit2, &bary, &earth2 );
332
333 /* work out frequency (assuming stationary at barycentre) */
334 taylorcoeff = 1.;
335 tmpdt = DT;
336 freqs->data[i] = fs->data[0];
337 for ( UINT4 k = 1; k < fs->length; k++ ) {
338 taylorcoeff /= ( REAL8 )k;
339 freqs->data[i] += taylorcoeff * fs->data[k] * tmpdt;
340 tmpdt *= DT;
341 }
342 freqs->data[i] *= freqharm;
343
344 /* get solar system doppler shift and add it on */
345 dfsolar->data[i] = ( emit2.deltaT - emit.deltaT ) * freqs->data[i];
346
347 /* get binary system doppler shift */
348 if ( PulsarCheckParam( params, "BINARY" ) ) {
349 binput.tb = time0 + emit.deltaT;
350 XLALGPSSetREAL8( &tgps, time0 );
351 get_earth_pos_vel( &earth, edat, &tgps );
352 binput.earth = earth;
353 XLALBinaryPulsarDeltaTNew( &boutput, &binput, params );
354
355 binput.tb = timep1 + emit2.deltaT;
356 XLALGPSSetREAL8( &tgps, timep1 );
357 get_earth_pos_vel( &earth, edat, &tgps );
358 binput.earth = earth2;
359 XLALBinaryPulsarDeltaTNew( &boutput2, &binput, params );
360
361 dfbinary->data[i] = ( boutput2.deltaT - boutput.deltaT ) * freqs->data[i];
362 } else {
363 dfbinary->data[i] = 0.;
364 }
365
366 dftotal->data[i] = dfbinary->data[i] + dfsolar->data[i];
367
368 freqs->data[i] += dftotal->data[i];
369 }
370}
void XLALBinaryPulsarDeltaTNew(BinaryPulsarOutput *output, BinaryPulsarInput *input, PulsarParameters *params)
function to calculate the binary system delay using new parameter structure
const char * program
#define __func__
log an I/O error, i.e.
int k
#define c
int LALgetopt_long(int argc, char *const *argv, const char *options, const struct LALoption *long_options, int *opt_index)
char * LALoptarg
#define no_argument
#define required_argument
const REAL8Vector * PulsarGetREAL8VectorParam(const PulsarParameters *pars, const CHAR *name)
Return a REAL8Vector parameter.
REAL8 PulsarGetREAL8ParamOrZero(const PulsarParameters *pars, const CHAR *name)
Return a REAL8 parameter if it exists, otherwise return zero.
int PulsarCheckParam(const PulsarParameters *pars, const CHAR *name)
Check for the existence of the parameter name in the PulsarParameters structure.
PulsarParameters * XLALReadTEMPOParFile(const CHAR *pulsarAndPath)
Read in the parameters from a TEMPO(2) parameter file into a PulsarParameters structure.
REAL8 PulsarGetREAL8Param(const PulsarParameters *pars, const CHAR *name)
Return a REAL8 parameter.
void PulsarFreeParams(PulsarParameters *pars)
Function to free memory from pulsar parameters.
const char * name
Definition: SearchTiming.c:93
#define fprintf
double e
TimeCorrectionData * XLALInitTimeCorrections(const CHAR *timeCorrectionFile)
An XLAL interface for reading a time correction file containing a table of values for converting betw...
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
int XLALBarycenter(EmissionTime *emit, const BarycenterInput *baryinput, const EarthState *earth)
Transforms from detector arrival time in GPS (as specified in the LIGOTimeGPS structure) to pulse em...
int XLALBarycenterEarthNew(EarthState *earth, const LIGOTimeGPS *tGPS, const EphemerisData *edat, const TimeCorrectionData *tdat, TimeCorrectionType ttype)
Computes the position and orientation of the Earth, at some arrival time, but unlike XLALBarycenterEa...
TimeCorrectionType
Enumerated type denoting the time system type to be produced in the solar system barycentring routine...
Definition: LALBarycenter.h:72
#define LAL_PC_SI
uint64_t UINT8
double REAL8
char CHAR
uint32_t UINT4
int32_t INT4
char char * XLALStringDuplicate(const char *s)
const LALDetector * XLALGetSiteInfo(const CHAR *name)
Find the site geometry-information 'LALDetector' for given a detector name (or prefix).
Definition: SFTnaming.c:218
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
#define XLAL_ERROR_VOID(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_EINVAL
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
dist
int T0
void get_earth_pos_vel(EarthState *earth, EphemerisData *edat, LIGOTimeGPS *tGPS)
Get the position and velocity of the Earth at a given time.
Definition: ppe_models.c:1215
Header file for the signal models functions used in parameter estimation code for known pulsar search...
TimeCorrectionType XLALAutoSetEphemerisFiles(CHAR **efile, CHAR **sfile, CHAR **tfile, PulsarParameters *pulsar, INT4 gpsstart, INT4 gpsend)
Automatically set the solar system ephemeris file based on environment variables and data time span.
Definition: ppe_utils.c:717
Header file for the helper functions for the parameter estimation code for known pulsar searches usin...
int main(int argc, char *argv[])
void get_freq(REAL8 start, REAL8 deltaT, REAL8 freqharm, PulsarParameters *params, BarycenterInput bary, EphemerisData *edat, TimeCorrectionData *tdat, TimeCorrectionType ttype, REAL8Vector *freqs, REAL8Vector *dfsolar, REAL8Vector *dfbinary, REAL8Vector *dftotal)
#define CODEUSAGE
void get_input_args(InputParams *inputParams, INT4 argc, CHAR *argv[])
Basic input structure to LALBarycenter.c.
REAL8 alpha
Source right ascension in ICRS J2000 coords (radians).
REAL8 dInv
1/(distance to source), in 1/sec.
LIGOTimeGPS tgps
input GPS arrival time.
REAL8 delta
Source declination in ICRS J2000 coords (radians)
LALDetector site
detector site info.
structure containing the input parameters for the binary delay function
REAL8 tb
Time of arrival (TOA) at the SSB.
EarthState earth
The current Earth state (for e.g.
structure containing the output parameters for the binary delay function
REAL8 deltaT
deltaT to add to TDB in order to account for binary
Basic output structure of LALBarycenterEarth.c.
Basic output structure produced by LALBarycenter.c.
REAL8 deltaT
(TDB) - (GPS)
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
REAL8 location[3]
int * flag
INT4 gpsNanoSeconds
The PulsarParameters structure to contain a set of pulsar parameters.
REAL8 * data
This structure will contain a vector of time corrections used during conversion from TT to TDB/TCB/Te...
double deltaT