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
SimulateTaylorCWTest.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Teviet Creighton
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/**
21 * \file
22 * \ingroup GenerateTaylorCW_h
23 * \author Creighton, T. D.
24 *
25 * \brief Generates a quasiperiodic waveform.
26 *
27 * ### Usage ###
28 *
29 * \code
30 * SimulateTaylorCWTest [-s sourcefile] [-r respfile] [-l site earthfile sunfile]
31 * [-o outfile] [-t sec nsec npt dt] [-h hsec hnsec fh]
32 * [-d debuglevel]
33 * \endcode
34 *
35 * ### Description ###
36 *
37 * This program generates a Taylor-parameterized continuous waveform as a
38 * function of time. The following option flags are accepted:
39 * <ul>
40 * <li>[<tt>-s</tt>] Reads source information from the file
41 * \c sourcefile, whose format is specified below. If absent, it
42 * injects a monochromatic wave with intrinsic frequency 100~Hz, strain
43 * amplitude 1000, and zero phase at the GPS epoch, arriving from the
44 * direction RA=00h00m, dec= \f$ 0^\circ \f$ .</li>
45 * <li>[<tt>-r</tt>] Reads a detector response function from the file
46 * \c respfile, whose format is specified below. If absent, it
47 * generates raw dimensionless strain.</li>
48 * <li>[<tt>-l</tt>] Sets the detector location and orientation.
49 * \c site must be one of the following character strings:
50 * \c LHO, \c LLO, \c VIRGO, \c GEO600, \c TAMA300,
51 * or \c CIT40. \c earthfile and \c sunfile are ephemeris
52 * files of the format expected by <tt>LALInitBarycenter()</tt>. If the
53 * <tt>-l</tt> option is not specified, a stationary (barycentric) detector
54 * aligned with the wave's + polarization is assumed.</li>
55 * <li>[<tt>-o</tt>] Writes the generated time series to the file
56 * \c outfile. If absent, the routines are exercised, but no output
57 * is written.</li>
58 * <li>[<tt>-t</tt>] Sets the timing information for the generated
59 * waveform: \c sec and \c nsec are integers specifying the start
60 * time in GPS seconds and nanoseconds, \c npt is the number of time
61 * samples generated, and \c dt is the sampling interval in seconds.
62 * If absent, <tt>-t 0 0 65536 9.765625e-4</tt> is assumed.</li>
63 * <li>[<tt>-h</tt>] Performs "ideal heterodyning" (phase subtraction)
64 * on the signal: \c hsec and \c hnsec are integers specifying an
65 * epoch of zero phase subtraction in GPS seconds and nanoseconds, and
66 * \c fh is the frequency of phase subtraction in Hz. If absent, no
67 * heterodyning is performed.</li>
68 * <li>[<tt>-d</tt>] Sets the debug level to \c debuglevel. If
69 * absent, level 0 is assumed.</li>
70 * </ul>
71 *
72 * \par Format for \c sourcefile:
73 * The source file consists
74 * of any number of lines of data, each specifying a continuous waveform.
75 * Each line consists of a GPS epoch where the frequency, phase, and
76 * Taylor coefficients are defined (\c INT8 nanoseconds), followed by
77 * 7 or more whitespace-delimited \c REAL8 numbers: the + and
78 * \f$ \times \f$ wave amplitudes (dimensionless strain) and polarization angle
79 * \f$ \psi \f$ (degrees), the right ascension and declination (degrees), the
80 * initial phase (degrees) and frequency (Hz), followed by zero or more
81 * Taylor coefficients \f$ f_k \f$ (Hz \f$ {}^k \f$ ). Note that the wave amplitudes
82 * and polarization angle are read as \c REAL8, but are later cast to
83 * \c REAL4.
84 *
85 * \par Format for \c respfile:
86 * The response function \f$ R(f) \f$
87 * gives the real and imaginary components of the transformation
88 * \e from ADC output \f$ o \f$ \e to tidal strain \f$ h \f$ via
89 * \f$ \tilde{h}(f)=R(f)\tilde{o}(f) \f$ . It is inverted internally to give
90 * the detector <em>transfer function</em> \f$ T(f)=1/R(f) \f$ . The format
91 * \c respfile is a header specifying the GPS epoch \f$ t_0 \f$ at which
92 * the response was taken (\c INT8 nanoseconds), the lowest frequency
93 * \f$ f_0 \f$ at which the response is given (\c REAL8 Hz), and the
94 * frequency sampling interval \f$ \Delta f \f$ (\c REAL8 Hz):
95 *
96 * <table><tr><td>
97 * <tt># epoch = </tt> \f$ t_0 \f$ </td></tr>
98 * <tr><td><tt># f0 = </tt> \f$ f_0 \f$ </td></tr>
99 * <tr><td><tt># deltaF = </tt> \f$ \Delta f \f$
100 * </td></tr></table>
101 *
102 * followed by two columns of \c REAL4 data giving the real
103 * and imaginary components of \f$ R(f_0+k\Delta f) \f$ .
104 *
105 * \par Format for \c outfile:
106 * The ouput files generated by
107 * this program consist of a two-line header giving the GPS epoch \f$ t_0 \f$
108 * of the first time sample (\c INT8 nanoseconds) and the sampling
109 * interval \f$ \Delta t \f$ (\c REAL8 seconds):
110 *
111 * <table><tr><td>
112 * <tt># epoch = </tt> \f$ t_0 \f$ </td></tr>
113 * <tr><td><tt># deltaT = </tt> \f$ \Delta t \f$
114 * </td></tr></table>
115 *
116 * followed by a single column of \c REAL4 data
117 * representing the undigitized output of the interferometer's
118 * gravitational-wave channel.
119 *
120 * ### Algorithm ###
121 *
122 * If a \c sourcefile is specified, this program first reads the
123 * epoch field, and then calls <tt>LALSReadVector()</tt> to read the
124 * remaining fields. (If no \c sourcefile is specified, the
125 * parameters are taken from <tt>\#define</tt>d constants.) The arguments
126 * are passed to <tt>LALGenerateTaylorCW()</tt> to generate frequency and
127 * phase timeseries. The required sampling resolution for these
128 * timeseries is estimated at \f$ 0.1/\sum_k\sqrt{|kf_0f_kT^{k-1}|} \f$ , where
129 * \f$ T \f$ is the duration of the waveform, to ensure that later
130 * interpolation of the waveforms will give phases accurate to well
131 * within a radian.
132 *
133 * The output from <tt>LALGenerateTaylorCW()</tt> is then passed to
134 * <tt>LALPulsarSimulateCoherentGW()</tt> to generate an output time series. If
135 * there are multiple lines in \c sourcefile, the procedure is
136 * repeated and the new waveforms added into the output. A warning is
137 * generated if the wave frequency exceeds the Nyquist rate for the
138 * output time series.
139 *
140 */
141
142/** \name Error Codes */
143/** @{ */
144#define SIMULATETAYLORCWTESTC_ENORM 0 /**< Normal exit */
145#define SIMULATETAYLORCWTESTC_ESUB 1 /**< Subroutine failed */
146#define SIMULATETAYLORCWTESTC_EARG 2 /**< Error parsing arguments */
147#define SIMULATETAYLORCWTESTC_EVAL 3 /**< Input argument out of valid range */
148#define SIMULATETAYLORCWTESTC_EFILE 4 /**< Could not open file */
149#define SIMULATETAYLORCWTESTC_EINPUT 5 /**< Error reading file */
150#define SIMULATETAYLORCWTESTC_EMEM 6 /**< Out of memory */
151#define SIMULATETAYLORCWTESTC_EPRINT 7 /**< Wrote past end of message string */
152/** @} */
153
154/** \cond DONT_DOXYGEN */
155#define SIMULATETAYLORCWTESTC_MSGENORM "Normal exit"
156#define SIMULATETAYLORCWTESTC_MSGESUB "Subroutine failed"
157#define SIMULATETAYLORCWTESTC_MSGEARG "Error parsing arguments"
158#define SIMULATETAYLORCWTESTC_MSGEVAL "Input argument out of valid range"
159#define SIMULATETAYLORCWTESTC_MSGEFILE "Could not open file"
160#define SIMULATETAYLORCWTESTC_MSGEINPUT "Error reading file"
161#define SIMULATETAYLORCWTESTC_MSGEMEM "Out of memory"
162#define SIMULATETAYLORCWTESTC_MSGEPRINT "Wrote past end of message string"
163
164
165#include <math.h>
166#include <stdlib.h>
167#include <lal/LALStdio.h>
168#include <lal/LALStdlib.h>
169#include <lal/LALConstants.h>
170#include <lal/StreamInput.h>
171#include <lal/AVFactories.h>
172#include <lal/SeqFactories.h>
173#include <lal/VectorOps.h>
174#include <lal/DetectorSite.h>
175#include <lal/Units.h>
176#include <lal/PulsarSimulateCoherentGW.h>
177#include <lal/GenerateTaylorCW.h>
178#include <lal/LALBarycenter.h>
179#include <lal/LALInitBarycenter.h>
180
181/* Default parameter settings. */
182#define EPOCH (0LL) /* about Jan. 1, 1990 */
183#define APLUS (1000.0)
184#define ACROSS (1000.0)
185#define RA (0.0)
186#define DEC (0.0)
187#define PSI (0.0)
188#define F0 (100.0)
189#define FH (0.0)
190#define PHI0 (0.0)
191#define SEC (0)
192#define NSEC (0)
193#define HSEC (0)
194#define HNSEC (0)
195#define DT (0.0009765625)
196#define NPT (65536)
197
198/* Usage format string. */
199#define USAGE "Usage: %s [-s sourcefile] [-o outfile] [-r respfile]\n" \
200"\t[-l site earthfile sunfile] [-o outfile] [-t sec nsec npt dt]\n" \
201"\t[-t hsec hnsec fh] [-d debuglevel]\n"
202
203/* Maximum output message length. */
204#define MSGLEN (1024)
205
206/* Upper cutoff frequency for default detector response function. */
207#define FSTOP (16384.0)
208
209/* Macros for printing errors and testing subroutines. */
210#define ERROR( code, msg, statement ) \
211do \
212if ( lalDebugLevel & LALERROR ) \
213{ \
214 XLALPrintError( "Error[0] %d: program %s, file %s, line %d, %s\n" \
215 " %s %s\n", (code), *argv, __FILE__, \
216 __LINE__, "$Id$", \
217 statement ? statement : "", (msg) ); \
218} \
219while (0)
220
221#define INFO( statement ) \
222do \
223if ( lalDebugLevel & LALINFO ) \
224{ \
225 XLALPrintError( "Info[0]: program %s, file %s, line %d, %s\n" \
226 " %s\n", *argv, __FILE__, __LINE__, \
227 "$Id$", (statement) ); \
228} \
229while (0)
230
231#define WARNING( statement ) \
232do \
233if ( lalDebugLevel & LALWARNING ) \
234{ \
235 XLALPrintError( "Warning[0]: program %s, file %s, line %d, %s\n" \
236 " %s\n", *argv, __FILE__, __LINE__, \
237 "$Id$", (statement) ); \
238} \
239while (0)
240
241#define SUB( func, statusptr ) \
242do \
243if ( (func), (statusptr)->statusCode ) \
244{ \
245 ERROR( SIMULATETAYLORCWTESTC_ESUB, SIMULATETAYLORCWTESTC_MSGESUB, \
246 "Function call \"" #func "\" failed:" ); \
247 return SIMULATETAYLORCWTESTC_ESUB; \
248} \
249while (0)
250
251#define CHECKVAL( val, lower, upper ) \
252do \
253if ( ( (val) < (lower) ) || ( (val) > (upper) ) ) \
254{ \
255 ERROR( SIMULATETAYLORCWTESTC_EVAL, SIMULATETAYLORCWTESTC_MSGEVAL, \
256 "Value of " #val " out of range:" ); \
257 if ( lalDebugLevel & LALERROR ) \
258 XLALPrintError( #val " = %f, range = [%f,%f]\n", (REAL8)(val), \
259 (REAL8)(lower), (REAL8)(upper) ); \
260 return SIMULATETAYLORCWTESTC_EVAL; \
261} \
262while (0)
263
264/* A function to convert INT8 nanoseconds to LIGOTimeGPS. */
265void
266I8ToLIGOTimeGPS( LIGOTimeGPS *output, INT8 input );
267
268/* A function to compute the combination function C(a+1,b+1). */
269UINT4
270choose( UINT4 a, UINT4 b );
271
272
273int
274main( int argc, char **argv )
275{
276 /* Command-line parsing variables. */
277 int arg; /* command-line argument counter */
278 static LALStatus stat; /* status structure */
279 CHAR *sourcefile = NULL; /* name of sourcefile */
280 CHAR *respfile = NULL; /* name of respfile */
281 CHAR *outfile = NULL; /* name of outfile */
282 CHAR *earthfile = NULL; /* name of earthfile */
283 CHAR *sunfile = NULL; /* name of sunfile */
284 CHAR *site = NULL; /* name of detector site */
285 INT4 npt = NPT; /* number of output samples */
286 INT4 sec = SEC, nsec = NSEC; /* GPS epoch of output */
287 INT4 hsec = HSEC, hnsec = HNSEC; /* GPS heterodyning epoch */
288 REAL8 dt = DT; /* output sampling interval */
289 REAL8 fh = FH; /* heterodyning frequency */
290
291 /* File reading variables. */
292 FILE *fp = NULL; /* generic file pointer */
293 BOOLEAN ok = 1; /* whether input format is correct */
294 INT8 epoch = EPOCH; /* epoch stored as an INT8 */
295
296 /* Other variables. */
297 UINT4 i, j; /* generic indecies */
298 INT8 tStart, tStop; /* start and stop times for waveform */
299 PulsarDetectorResponse detector; /* the detector in question */
300 REAL4TimeSeries output; /* detector output */
301
302
303 /*******************************************************************
304 * ARGUMENT PARSING (arg stores the current position) *
305 *******************************************************************/
306
307 arg = 1;
308 while ( arg < argc ) {
309 /* Parse source file option. */
310 if ( !strcmp( argv[arg], "-s" ) ) {
311 if ( argc > arg + 1 ) {
312 arg++;
313 sourcefile = argv[arg++];
314 } else {
316 SIMULATETAYLORCWTESTC_MSGEARG, 0 );
317 XLALPrintError( USAGE, *argv );
319 }
320 }
321 /* Parse response file option. */
322 else if ( !strcmp( argv[arg], "-r" ) ) {
323 if ( argc > arg + 1 ) {
324 arg++;
325 respfile = argv[arg++];
326 } else {
328 SIMULATETAYLORCWTESTC_MSGEARG, 0 );
329 XLALPrintError( USAGE, *argv );
331 }
332 }
333 /* Parse output file option. */
334 else if ( !strcmp( argv[arg], "-o" ) ) {
335 if ( argc > arg + 1 ) {
336 arg++;
337 outfile = argv[arg++];
338 } else {
340 SIMULATETAYLORCWTESTC_MSGEARG, 0 );
341 XLALPrintError( USAGE, *argv );
343 }
344 }
345 /* Parse detector location option. */
346 else if ( !strcmp( argv[arg], "-l" ) ) {
347 if ( argc > arg + 3 ) {
348 arg++;
349 site = argv[arg++];
350 earthfile = argv[arg++];
351 sunfile = argv[arg++];
352 } else {
354 SIMULATETAYLORCWTESTC_MSGEARG, 0 );
355 XLALPrintError( USAGE, *argv );
357 }
358 }
359 /* Parse output timing option. */
360 else if ( !strcmp( argv[arg], "-t" ) ) {
361 if ( argc > arg + 4 ) {
362 arg++;
363 sec = atoi( argv[arg++] );
364 nsec = atoi( argv[arg++] );
365 npt = atoi( argv[arg++] );
366 dt = atof( argv[arg++] );
367 } else {
369 SIMULATETAYLORCWTESTC_MSGEARG, 0 );
370 XLALPrintError( USAGE, *argv );
372 }
373 }
374 /* Parse heterodyning option. */
375 else if ( !strcmp( argv[arg], "-h" ) ) {
376 if ( argc > arg + 4 ) {
377 arg++;
378 hsec = atoi( argv[arg++] );
379 hnsec = atoi( argv[arg++] );
380 fh = atof( argv[arg++] );
381 } else {
383 SIMULATETAYLORCWTESTC_MSGEARG, 0 );
384 XLALPrintError( USAGE, *argv );
386 }
387 }
388 /* Parse debug level option. */
389 else if ( !strcmp( argv[arg], "-d" ) ) {
390 if ( argc > arg + 1 ) {
391 arg++;
392 } else {
394 SIMULATETAYLORCWTESTC_MSGEARG, 0 );
395 XLALPrintError( USAGE, *argv );
397 }
398 }
399 /* Check for unrecognized options. */
400 else if ( argv[arg][0] == '-' ) {
402 SIMULATETAYLORCWTESTC_MSGEARG, 0 );
403 XLALPrintError( USAGE, *argv );
405 }
406 } /* End of argument parsing loop. */
407
408 /* Make sure that values won't crash the system or anything. */
409 CHECKVAL( dt, LAL_REAL8_MIN, LAL_REAL8_MAX );
410 CHECKVAL( npt, 0, 2147483647 );
411
412
413 /*******************************************************************
414 * SETUP *
415 *******************************************************************/
416
417 /* Set up output structure and wave start and stop times. */
418 {
419 INT8 epochOut = nsec; /* output epoch */
420 epochOut += 1000000000LL * sec;
421 tStart = epochOut - 1000000000LL;
422 memset( &output, 0, sizeof( REAL4TimeSeries ) );
423 I8ToLIGOTimeGPS( &( output.epoch ), epochOut );
424 output.deltaT = dt;
425 snprintf( output.name, LALNameLength, "Taylor CW waveform" );
426 SUB( LALSCreateVector( &stat, &( output.data ), npt ), &stat );
427 memset( output.data->data, 0, npt * sizeof( REAL4 ) );
428 tStop = epochOut + 1000000000LL * ( INT8 )( dt * npt + 1.0 );
429 output.f0 = fh;
430 }
431 /* Adjust wave start and stop times so that they will almost
432 certainly cover the output timespan even after barycentring. */
433 if ( site ) {
434 tStart -= ( INT8 )( ( 1.1e9 ) * LAL_AU_SI / LAL_C_SI );
435 tStop += ( INT8 )( ( 1.1e9 ) * LAL_AU_SI / LAL_C_SI );
436 }
437
438
439 /* Set up detector structure. */
440 memset( &detector, 0, sizeof( PulsarDetectorResponse ) );
443 if ( !( transfer ) ) {
445 SIMULATETAYLORCWTESTC_MSGEMEM, 0 );
447 }
448 detector.heterodyneEpoch.gpsSeconds = hsec;
449 detector.heterodyneEpoch.gpsNanoSeconds = hnsec;
450 memset( transfer, 0, sizeof( COMPLEX8FrequencySeries ) );
451 if ( respfile ) {
452 REAL4VectorSequence *resp = NULL; /* response as vector sequence */
453 COMPLEX8Vector *response = NULL; /* response as complex vector */
454 COMPLEX8Vector *unity = NULL; /* vector of complex 1's */
455
456 if ( ( fp = fopen( respfile, "r" ) ) == NULL ) {
458 SIMULATETAYLORCWTESTC_MSGEFILE, respfile );
460 }
461
462 /* Read header. */
463 ok &= ( fscanf( fp, "# epoch = %" LAL_INT8_FORMAT "\n", &epoch ) == 1 );
464 I8ToLIGOTimeGPS( &( transfer->epoch ), epoch );
465 ok &= ( fscanf( fp, "# f0 = %lf\n", &( transfer->f0 ) )
466 == 1 );
467 ok &= ( fscanf( fp, "# deltaF = %lf\n",
468 &( transfer->deltaF ) ) == 1 );
469 if ( !ok ) {
471 SIMULATETAYLORCWTESTC_MSGEINPUT, respfile );
473 }
474
475 /* Read and convert body to a COMPLEX8Vector. */
476 SUB( LALSReadVectorSequence( &stat, &resp, fp ), &stat );
477 fclose( fp );
478 if ( resp->vectorLength != 2 ) {
480 SIMULATETAYLORCWTESTC_MSGEINPUT, respfile );
482 }
483 SUB( LALCCreateVector( &stat, &response, resp->length ), &stat );
484 memcpy( response->data, resp->data, 2 * resp->length * sizeof( REAL4 ) );
485 SUB( LALSDestroyVectorSequence( &stat, &resp ), &stat );
486
487 /* Convert response function to a transfer function. */
488 SUB( LALCCreateVector( &stat, &unity, response->length ), &stat );
489 for ( i = 0; i < response->length; i++ ) {
490 unity->data[i] = 1.0;
491 }
492 SUB( LALCCreateVector( &stat, &( transfer->data ),
493 response->length ), &stat );
494 SUB( LALCCVectorDivide( &stat, transfer->data, unity,
495 response ), &stat );
496 SUB( LALCDestroyVector( &stat, &response ), &stat );
497 SUB( LALCDestroyVector( &stat, &unity ), &stat );
498 }
499 /* No response file, so generate a unit response. */
500 else {
501 I8ToLIGOTimeGPS( &( transfer->epoch ), EPOCH );
502 transfer->f0 = 0.0;
503 transfer->deltaF = FSTOP;
504 SUB( LALCCreateVector( &stat, &( transfer->data ), 2 ),
505 &stat );
506 transfer->data->data[0] = 1.0;
507 transfer->data->data[1] = 1.0;
508 }
509 LALDetector *lsite = NULL;
510 EphemerisData *lephem = NULL;
511 if ( site ) {
512 /* Set detector location. */
513 if ( !strcmp( site, "LHO" ) ) {
515 } else if ( !strcmp( site, "LLO" ) ) {
517 } else if ( !strcmp( site, "VIRGO" ) ) {
519 } else if ( !strcmp( site, "GEO600" ) ) {
521 } else if ( !strcmp( site, "TAMA300" ) ) {
523 } else if ( !strcmp( site, "CIT40" ) ) {
525 } else {
527 SIMULATETAYLORCWTESTC_MSGEVAL, "Unrecognized site:" );
528 if ( lalDebugLevel & LALERROR ) {
529 XLALPrintError( "%s", site );
530 }
532 }
533 lsite = ( LALDetector * )LALMalloc( sizeof( LALDetector ) );
534 if ( !( lsite ) ) {
536 SIMULATETAYLORCWTESTC_MSGEMEM, 0 );
538 }
539 *( lsite ) = lalCachedDetectors[i];
540
541 /* Read ephemerides. */
542 XLAL_CHECK_MAIN( ( lephem = XLALInitBarycenter( earthfile, sunfile ) ) != NULL, XLAL_EFUNC );
543 }
544
545
546 /* Set up units for the above structures. */
547 output.sampleUnits = lalADCCountUnit;
548 if ( XLALUnitDivide( &( transfer->sampleUnits ),
549 &lalADCCountUnit, &lalStrainUnit ) == NULL ) {
550 return LAL_EXLAL;
551 }
552
553 detector.transfer = transfer;
554 detector.site = lsite;
555 detector.ephemerides = lephem;
556
557 /*******************************************************************
558 * OUTPUT GENERATION *
559 *******************************************************************/
560
561 /* Open sourcefile. */
562 if ( sourcefile ) {
563 if ( ( fp = fopen( sourcefile, "r" ) ) == NULL ) {
565 SIMULATETAYLORCWTESTC_MSGEFILE, sourcefile );
567 }
568 }
569
570 /* For each line in the sourcefile... */
571 while ( ok ) {
572 TaylorCWParamStruc params; /* wave generation parameters */
573 PulsarCoherentGW waveform; /* amplitude and phase structure */
574 REAL4TimeSeries signalvec; /* GW signal */
575 CHAR message[MSGLEN]; /* warning/info messages */
576
577 /* Initialize output structures. */
578 memset( &waveform, 0, sizeof( PulsarCoherentGW ) );
579 signalvec = output;
580 signalvec.data = NULL;
581
582 /* Read and convert input line. */
583 memset( &params, 0, sizeof( TaylorCWParamStruc ) );
584 I8ToLIGOTimeGPS( &( params.epoch ), tStart );
585 params.position.system = COORDINATESYSTEM_EQUATORIAL;
586 if ( sourcefile ) {
587 REAL8Vector *input = NULL; /* input parameters */
588 ok &= ( fscanf( fp, "%" LAL_INT8_FORMAT, &epoch ) == 1 );
589 if ( ok ) {
590 SUB( LALDReadVector( &stat, &input, fp, 1 ), &stat );
591 ok &= ( input->length > 6 );
592 }
593 if ( ok ) {
594 params.aPlus = input->data[0];
595 params.aCross = input->data[1];
596 params.psi = LAL_PI * input->data[2] / 180.0;
597 params.position.longitude = LAL_PI * input->data[3] / 180.0;
598 params.position.latitude = LAL_PI * input->data[4] / 180.0;
599 params.phi0 = LAL_PI * input->data[5] / 180.0;
600 params.f0 = input->data[6];
601 if ( input->length > 7 ) {
602 SUB( LALDCreateVector( &stat, &( params.f ),
603 input->length - 7 ), &stat );
604 for ( i = 0; i < input->length - 7; i++ ) {
605 params.f->data[i] = input->data[i + 7];
606 }
607 }
608 }
609 if ( input ) {
610 SUB( LALDDestroyVector( &stat, &input ), &stat );
611 }
612 } else {
613 params.aPlus = APLUS;
614 params.aCross = ACROSS;
615 params.psi = PSI;
616 params.position.longitude = RA;
617 params.position.latitude = DEC;
618 params.phi0 = PHI0;
619 params.f0 = F0;
620 }
621
622 /* Adjust frequency and spindown terms to the actual wave start
623 time. */
624 if ( params.f ) {
625 UINT4 length = params.f->length; /* number of coefficients */
626 REAL8 t = ( 1.0e-9 ) * ( REAL4 )( tStart - epoch ); /* time shift */
627 REAL8 tN = 1.0; /* t raised to various powers */
628 REAL8 fFac = 1.0; /* fractional change in frequency */
629 REAL8 tFac = 1.0; /* time integral of fFac */
630 REAL8 *fData = params.f->data; /* pointer to coeficients */
631 for ( i = 0; i < length; i++ ) {
632 REAL8 tM = 1.0; /* t raised to various powers */
633 fFac += fData[i] * ( tN *= t );
634 tFac += fData[i] * tN / ( i + 2.0 );
635 for ( j = i + 1; j < length; j++ ) {
636 fData[i] += choose( j + 1, i + 1 ) * fData[j] * ( tM *= t );
637 }
638 }
639 params.phi0 += LAL_TWOPI * params.f0 * t * tFac;
640 params.f0 *= fFac;
641 for ( i = 0; i < length; i++ ) {
642 fData[i] /= fFac;
643 }
644 }
645
646 if ( ok ) {
647 REAL4 *sigData, *outData; /* pointers to time series data */
648 REAL4 t = ( 1.0e-9 ) * ( REAL4 )( tStop - tStart ); /* duration */
649 REAL4 dtInv = 0.0; /* sampling rate for waveform */
650
651 /* Set remaining parameters, and generate waveform. */
652 if ( params.f ) {
653 REAL4 tN = 1.0; /* t raised to various powers */
654 for ( i = 0; i < params.f->length; i++ ) {
655 dtInv += sqrt( ( i + 1 ) * fabs( params.f->data[i] ) * tN );
656 tN *= t;
657 }
658 dtInv *= 10.0 * sqrt( fabs( params.f0 ) );
659 }
660 if ( dtInv < 1.0 / t ) {
661 params.deltaT = t;
662 params.length = 2;
663 } else {
664 params.deltaT = 1.0 / dtInv;
665 params.length = ( UINT4 )( t * dtInv ) + 2;
666 }
667 SUB( LALGenerateTaylorCW( &stat, &waveform, &params ), &stat );
668 if ( params.dfdt > 2.0 ) {
669 /* snprintf() can't seem to print floating-point formats.
670 snprintf( message, MSGLEN,
671 "Waveform sampling interval is too large:\n"
672 "\tmaximum df*dt = %f", params.dfdt );
673 */
674 INT4 code = sprintf( message,
675 "Waveform sampling interval is too large:\n"
676 "\tmaximum df*dt = %f", params.dfdt );
677 if ( code >= MSGLEN || code < 0 ) {
679 SIMULATETAYLORCWTESTC_MSGEPRINT, 0 );
681 }
682 WARNING( message );
683 }
684 SUB( LALSCreateVector( &stat, &( signalvec.data ), npt ), &stat );
685 SUB( LALPulsarSimulateCoherentGW( &stat, &signalvec, &waveform,
686 &detector ), &stat );
687 if ( params.f ) {
688 SUB( LALDDestroyVector( &stat, &( params.f ) ), &stat );
689 }
690
691 /* Inject waveform into output. */
692 sigData = signalvec.data->data;
693 outData = output.data->data;
694 for ( i = 0; i < ( UINT4 )( npt ); i++ ) {
695 outData[i] += sigData[i];
696 }
697
698 /* Clean up memory. */
699 SUB( LALSDestroyVectorSequence( &stat, &( waveform.a->data ) ),
700 &stat );
701 SUB( LALSDestroyVector( &stat, &( waveform.f->data ) ),
702 &stat );
703 SUB( LALDDestroyVector( &stat, &( waveform.phi->data ) ),
704 &stat );
705 LALFree( waveform.a );
706 LALFree( waveform.f );
707 LALFree( waveform.phi );
708 SUB( LALSDestroyVector( &stat, &( signalvec.data ) ), &stat );
709 }
710
711 /* Inject only one signal if there is no sourcefile. */
712 if ( !sourcefile ) {
713 ok = 0;
714 }
715 }
716
717 /* Input file is exhausted (or has a badly-formatted line ). */
718 if ( sourcefile ) {
719 fclose( fp );
720 }
721
722
723 /*******************************************************************
724 * CLEANUP *
725 *******************************************************************/
726
727 /* Print output file. */
728 if ( outfile ) {
729 if ( ( fp = fopen( outfile, "w" ) ) == NULL ) {
731 SIMULATETAYLORCWTESTC_MSGEFILE, outfile );
733 }
734 epoch = 1000000000LL * ( INT8 )( output.epoch.gpsSeconds );
735 epoch += ( INT8 )( output.epoch.gpsNanoSeconds );
736 fprintf( fp, "# epoch = %" LAL_INT8_FORMAT "\n", epoch );
737 fprintf( fp, "# deltaT = %23.16e\n", output.deltaT );
738 for ( i = 0; i < output.data->length; i++ ) {
739 fprintf( fp, "%16.9e\n", ( REAL4 )( output.data->data[i] ) );
740 }
741 fclose( fp );
742 }
743
744 /* Destroy remaining memory. */
745 SUB( LALSDestroyVector( &stat, &( output.data ) ), &stat );
746 SUB( LALCDestroyVector( &stat, &( transfer->data ) ),
747 &stat );
748 LALFree( transfer );
749 if ( site ) {
750 XLALDestroyEphemerisData( lephem );
751 LALFree( lsite );
752 }
753
754 /* Done! */
756 INFO( SIMULATETAYLORCWTESTC_MSGENORM );
758}
759
760
761/* A function to convert INT8 nanoseconds to LIGOTimeGPS. */
762void
763I8ToLIGOTimeGPS( LIGOTimeGPS *output, INT8 input )
764{
765 INT8 s = input / 1000000000LL;
766 output->gpsSeconds = ( INT4 )( s );
767 output->gpsNanoSeconds = ( INT4 )( input - 1000000000LL * s );
768 return;
769}
770
771
772/* A function to compute the combination function C(a,b). */
773UINT4
774choose( UINT4 a, UINT4 b )
775{
776 UINT4 numer = 1;
777 UINT4 denom = 1;
778 UINT4 lal_index = b + 1;
779 while ( --lal_index ) {
780 numer *= a - b + lal_index;
781 denom *= lal_index;
782 }
783 return numer / denom;
784}
785/** \endcond */
LALDetectorIndexLLODIFF
LALDetectorIndexGEO600DIFF
LALDetectorIndexTAMA300DIFF
LALDetectorIndexVIRGODIFF
LALDetectorIndexCIT40DIFF
LALDetectorIndexLHODIFF
#define SUB(func, statusptr)
#define ERROR(code, msg, statement)
#define INFO(statement)
static UINT4 choose(UINT4 a, UINT4 b)
#define PHI0
#define PSI
int j
void LALCheckMemoryLeaks(void)
#define LALMalloc(n)
#define LALFree(p)
#define LAL_EXLAL
const WeaveSearchTimingDenominator denom
Definition: SearchTiming.c:95
#define SIMULATETAYLORCWTESTC_EINPUT
Error reading file.
#define SIMULATETAYLORCWTESTC_EARG
Error parsing arguments.
#define SIMULATETAYLORCWTESTC_ENORM
Normal exit.
#define SIMULATETAYLORCWTESTC_EFILE
Could not open file.
#define SIMULATETAYLORCWTESTC_EVAL
Input argument out of valid range.
#define SIMULATETAYLORCWTESTC_EMEM
Out of memory.
#define SIMULATETAYLORCWTESTC_EPRINT
Wrote past end of message string.
#define fscanf
#define fprintf
int main(int argc, char **argv)
void LALSDestroyVectorSequence(LALStatus *status, REAL4VectorSequence **vectorSequence)
const LALDetector lalCachedDetectors[LAL_NUM_DETECTORS]
void LALGenerateTaylorCW(LALStatus *stat, PulsarCoherentGW *output, TaylorCWParamStruc *params)
Computes a Taylor-parametrized continuous waveform.
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
#define LAL_TWOPI
unsigned char BOOLEAN
double REAL8
int64_t INT8
char CHAR
uint32_t UINT4
int32_t INT4
float REAL4
LALNameLength
LALERROR
void LALPulsarSimulateCoherentGW(LALStatus *stat, REAL4TimeSeries *output, PulsarCoherentGW *CWsignal, PulsarDetectorResponse *detector)
Computes the response of a detector to a coherent gravitational wave.
static const INT4 a
COORDINATESYSTEM_EQUATORIAL
void LALDReadVector(LALStatus *status, REAL8Vector **vector, FILE *stream, BOOLEAN strict)
void LALSReadVectorSequence(LALStatus *status, REAL4VectorSequence **sequence, FILE *stream)
const LALUnit lalStrainUnit
const LALUnit lalADCCountUnit
LALUnit * XLALUnitDivide(LALUnit *output, const LALUnit *unit1, const LALUnit *unit2)
void LALCCreateVector(LALStatus *, COMPLEX8Vector **, UINT4)
void LALDCreateVector(LALStatus *, REAL8Vector **, UINT4)
void LALCDestroyVector(LALStatus *, COMPLEX8Vector **)
void LALDDestroyVector(LALStatus *, REAL8Vector **)
void LALSDestroyVector(LALStatus *, REAL4Vector **)
void LALSCreateVector(LALStatus *, REAL4Vector **, UINT4)
void LALCCVectorDivide(LALStatus *status, COMPLEX8Vector *out, const COMPLEX8Vector *in1, const COMPLEX8Vector *in2)
#define XLAL_CHECK_MAIN(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_EFUNC
message
waveform
DEC
RA
#define F0
COMPLEX8Sequence * data
COMPLEX8 * data
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
This structure stores a representation of a plane gravitational wave propagating from a particular po...
This structure contains information required to determine the response of a detector to a gravitation...
REAL4Sequence * data
REAL4 * data
REAL8 * data
This structure stores the parameters for constructing a gravitational waveform with a Taylor-polynomi...
LIGOTimeGPS epoch
double psi
enum @4 site
enum @1 epoch