Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
GeneratePPNInspiralTest.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 * \author Creighton, T. D.
22 * \file
23 * \ingroup GeneratePPNInspiral_h
24 *
25 * \brief Generates a parametrized post-Newtonian inspiral waveform.
26 *
27 * ### Usage ###
28 *
29 * \code
30 * GeneratePPNInspiralTest [-m m1 m2] [-r dist] [-i inc phii] [-f f_min f_max]
31 * [-t dt] [-w deltat] [-p order] [-d debuglevel] [-o outfile]
32 * \endcode
33 *
34 * ### Description ###
35 *
36 * This program generates the amplitude, phase, and frequency of a
37 * post-Newtonian inspiral waveform as functions of time. The following
38 * option flags are accepted:
39 * <ul>
40 * <li>[<tt>-m</tt>] Sets the binary masses to \c m1 and \c m2
41 * solar massses (default values: \f$1.4M_\odot\f$).</li>
42 * <li>[<tt>-r</tt>] Sets the binary system distance to \c dist kpc
43 * (default value: 8.5kpc).</li>
44 * <li>[<tt>-i</tt>] Sets the inclination and \e initial phase
45 * angles to \c inc and \c phii degrees (default values:
46 * 0\ degrees).</li>
47 * <li>[<tt>-f</tt>] Sets the initial and final wave frequencies to
48 * \c f_min and \c f_max Hz (default values: 40Hz and 500Hz).</li>
49 * <li>[<tt>-t</tt>] Sets the waveform sampling interval to \c dt
50 * seconds (default value: 0.01s).</li>
51 * <li>[<tt>-w</tt>] Generates actual waveforms rather than phase and
52 * amplitude functions, sampled at intervals of \c deltat seconds (no
53 * default).</li>
54 * <li>[<tt>-p</tt>] Sets the post\f${}^{n/2}\f$-Newtonian order to
55 * \f$n=\f$\c order (default value: \f$n=4\f$).</li>
56 * <li>[<tt>-d</tt>] Sets the debug level to \c debuglevel (default
57 * value:\ 0).</li>
58 * <li>[<tt>-o</tt>] Sets the output filename to \c outfile (by
59 * default no output is produced).</li>
60 * </ul>
61 *
62 * ### Algorithm ###
63 *
64 * This program simply parses the command line, sets the appropriate
65 * fields of a \c PPNParamStruc, and passes it in to
66 * <tt>LALGeneratePPNInspiral()</tt>. No maximum waveform length is
67 * specified; the function will allocate as much data as necessary.
68 *
69 * If the <tt>-w</tt> \e and <tt>-o</tt> options are given, the
70 * amplitude, phase, and frequency are generated as above, but are then
71 * resampled at intervals \c deltat to generate actual wave output.
72 *
73 */
74
75/** \name Error Codes */
76/** @{ */
77#define GENERATEPPNINSPIRALTESTC_ENORM 0 /**< Normal exit */
78#define GENERATEPPNINSPIRALTESTC_ESUB 1 /**< Subroutine failed */
79#define GENERATEPPNINSPIRALTESTC_EARG 2 /**< Error parsing arguments */
80#define GENERATEPPNINSPIRALTESTC_EVAL 3 /**< Input argument out of valid range */
81#define GENERATEPPNINSPIRALTESTC_EFILE 4 /**< Could not open file */
82#define GENERATEPPNINSPIRALTESTC_EPRINT 5 /**< Wrote past end of message string */
83/** @} */
84
85/** \cond DONT_DOXYGEN */
86#define GENERATEPPNINSPIRALTESTC_MSGENORM "Normal exit"
87#define GENERATEPPNINSPIRALTESTC_MSGESUB "Subroutine failed"
88#define GENERATEPPNINSPIRALTESTC_MSGEARG "Error parsing arguments"
89#define GENERATEPPNINSPIRALTESTC_MSGEVAL "Input argument out of valid range"
90#define GENERATEPPNINSPIRALTESTC_MSGEFILE "Could not open file"
91#define GENERATEPPNINSPIRALTESTC_MSGEPRINT "Wrote past end of message string"
92
93#include <math.h>
94#include <stdlib.h>
95#include <lal/Date.h>
96#include <lal/LALStdio.h>
97#include <lal/LALStdlib.h>
98#include <lal/LALConstants.h>
99#include <lal/AVFactories.h>
100#include <lal/SeqFactories.h>
101#include <lal/SimulateCoherentGW.h>
102#include <lal/GeneratePPNInspiral.h>
103
104/* Default parameter settings. */
105#define EPOCH (315187200000000000LL) /* about Jan. 1, 1990 */
106#define M1 (1.4)
107#define M2 (1.4)
108#define DIST (8.5)
109#define INC (0.0)
110#define PHI (0.0)
111#define FMIN (40.0)
112#define FMAX (500.0)
113#define DT (0.01)
114#define ORDER (4)
115
116/* Usage format string. */
117#define USAGE "Usage: %s [-m m1 m2] [-r dist] [-i inc phii]\n\t[-f f_min f_max] [-t dt] [-w deltat] [-p order] [-d debuglevel] [-o outfile]\n"
118
119/* Maximum output message length. */
120#define MSGLENGTH (1024)
121
122
123/* Macros for printing errors and testing subroutines. */
124#define ERROR( code, msg, statement ) \
125do \
126if ( lalDebugLevel & LALERROR ) \
127{ \
128 LALPrintError( "Error[0] %d: program %s, file %s, line %d, %s\n" \
129 " %s %s\n", (code), *argv, __FILE__, \
130 __LINE__, "$Id$", \
131 statement ? statement : "", (msg) ); \
132} \
133while (0)
134
135#define INFO( statement ) \
136do \
137if ( lalDebugLevel & LALINFO ) \
138{ \
139 LALPrintError( "Info[0]: program %s, file %s, line %d, %s\n" \
140 " %s\n", *argv, __FILE__, __LINE__, \
141 "$Id$", (statement) ); \
142} \
143while (0)
144
145#define WARNING( statement ) \
146do \
147if ( lalDebugLevel & LALWARNING ) \
148{ \
149 LALPrintError( "Warning[0]: program %s, file %s, line %d, %s\n" \
150 " %s\n", *argv, __FILE__, __LINE__, \
151 "$Id$", (statement) ); \
152} \
153while (0)
154
155#define SUB( func, statusptr ) \
156do \
157if ( (func), (statusptr)->statusCode ) \
158{ \
159 ERROR( GENERATEPPNINSPIRALTESTC_ESUB, \
160 GENERATEPPNINSPIRALTESTC_MSGESUB, \
161 "Function call \"" #func "\" failed:" ); \
162 return GENERATEPPNINSPIRALTESTC_ESUB; \
163} \
164while (0)
165
166#define CHECKVAL( val, lower, upper ) \
167do \
168if ( ( (val) < (lower) ) || ( (val) > (upper) ) ) \
169{ \
170 ERROR( GENERATEPPNINSPIRALTESTC_EVAL, \
171 GENERATEPPNINSPIRALTESTC_MSGEVAL, \
172 "Value of " #val " out of range:" ); \
173 LALPrintError( #val " = %f, range = [%f,%f]\n", (REAL8)(val), \
174 (REAL8)(lower), (REAL8)(upper) ); \
175 return GENERATEPPNINSPIRALTESTC_EVAL; \
176} \
177while (0)
178
179
180int
181main(int argc, char **argv)
182{
183 /* Command-line parsing variables. */
184 int arg; /* command-line argument counter */
185 static LALStatus stat; /* status structure */
186 CHAR *outfile = NULL; /* name of outfile */
187 REAL4 m1 = M1, m2 = M2; /* binary masses */
188 REAL4 dist = DIST; /* binary distance */
189 REAL4 inc = 0.0, phii = 0.0; /* inclination and coalescence phase */
190 REAL4 f_min = FMIN, f_max=FMAX; /* start and stop frequencies */
191 REAL8 dt = DT; /* sampling interval */
192 REAL8 deltat = 0.0; /* wave sampling interval */
193 INT4 order = ORDER; /* PN order */
194
195 /* Other variables. */
196 UINT4 i; /* index */
197 CHAR message[MSGLENGTH]; /* signal generation output message */
198 PPNParamStruc params; /* input parameters */
199 CoherentGW waveform; /* output waveform */
200 FILE *fp; /* output file pointer */
201
202
203 /*******************************************************************
204 * ARGUMENT PARSING (arg stores the current position) *
205 *******************************************************************/
206
207 arg = 1;
208 while ( arg < argc ) {
209 /* Parse mass option. */
210 if ( !strcmp( argv[arg], "-m" ) ) {
211 if ( argc > arg + 2 ) {
212 arg++;
213 m1 = atof( argv[arg++] );
214 m2 = atof( argv[arg++] );
215 }else{
217 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
218 LALPrintError( USAGE, *argv );
220 }
221 }
222 /* Parse distance option. */
223 else if ( !strcmp( argv[arg], "-r" ) ) {
224 if ( argc > arg + 1 ) {
225 arg++;
226 dist = atof( argv[arg++] );
227 }else{
229 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
230 LALPrintError( USAGE, *argv );
232 }
233 }
234 /* Parse angles option. */
235 else if ( !strcmp( argv[arg], "-i" ) ) {
236 if ( argc > arg + 2 ) {
237 arg++;
238 inc = atof( argv[arg++] )*LAL_PI/180.0;
239 phii = atof( argv[arg++] )*LAL_PI/180.0;
240 }else{
242 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
243 LALPrintError( USAGE, *argv );
245 }
246 }
247 /* Parse frequency option. */
248 else if ( !strcmp( argv[arg], "-f" ) ) {
249 if ( argc > arg + 2 ) {
250 arg++;
251 f_min = atof( argv[arg++] );
252 f_max = atof( argv[arg++] );
253 }else{
255 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
256 LALPrintError( USAGE, *argv );
258 }
259 }
260 /* Parse sampling time option. */
261 else if ( !strcmp( argv[arg], "-t" ) ) {
262 if ( argc > arg + 1 ) {
263 arg++;
264 dt = atof( argv[arg++] );
265 }else{
267 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
268 LALPrintError( USAGE, *argv );
270 }
271 }
272 /* Parse waveform sampling time option. */
273 else if ( !strcmp( argv[arg], "-w" ) ) {
274 if ( argc > arg + 1 ) {
275 arg++;
276 deltat = atof( argv[arg++] );
277 }else{
279 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
280 LALPrintError( USAGE, *argv );
282 }
283 }
284 /* Parse PN order option. */
285 else if ( !strcmp( argv[arg], "-p" ) ) {
286 if ( argc > arg + 1 ) {
287 arg++;
288 order = atoi( argv[arg++] );
289 }else{
291 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
292 LALPrintError( USAGE, *argv );
294 }
295 }
296 /* Parse output file option. */
297 else if ( !strcmp( argv[arg], "-o" ) ) {
298 if ( argc > arg + 1 ) {
299 arg++;
300 outfile = argv[arg++];
301 }else{
303 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
304 LALPrintError( USAGE, *argv );
306 }
307 }
308 /* Parse debug level option. */
309 else if ( !strcmp( argv[arg], "-d" ) ) {
310 if ( argc > arg + 1 ) {
311 arg++;
312 }else{
314 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
315 LALPrintError( USAGE, *argv );
317 }
318 }
319 /* Check for unrecognized options. */
320 else if ( argv[arg][0] == '-' ) {
322 GENERATEPPNINSPIRALTESTC_MSGEARG, 0 );
323 LALPrintError( USAGE, *argv );
325 }
326 } /* End of argument parsing loop. */
327
328 /* Make sure that values won't crash the system or anything. */
329 CHECKVAL( order, -1, 5 );
330 CHECKVAL( dt, LAL_REAL4_MIN, LAL_REAL4_MAX );
331 CHECKVAL( deltat, 0.0, LAL_REAL4_MAX );
332
333
334 /*******************************************************************
335 * INPUT SETUP *
336 *******************************************************************/
337
338 /* Fixed parameters. */
339 params.position.latitude = params.position.longitude = 0.0;
340 params.position.system = COORDINATESYSTEM_EQUATORIAL;
341 params.psi = 0.0;
342 params.lengthIn = 0;
343
344 /* Variable parameters. */
345 XLALINT8NSToGPS( &(params.epoch), EPOCH );
346 params.deltaT = dt;
347 params.mTot = m1 + m2;
348 params.eta = m1*m2/( params.mTot*params.mTot );
349 params.inc = inc;
350 params.phi = 0.0;
351 params.d = dist*LAL_PC_SI*1.0e3;
352 params.fStartIn = f_min;
353 params.fStopIn = f_max;
354
355 /* PPN parameter. */
356 params.ppn = NULL;
357 SUB( LALSCreateVector( &stat, &(params.ppn), order + 1 ), &stat );
358 params.ppn->data[0] = 1.0;
359 if ( order > 0 )
360 params.ppn->data[1] = 0.0;
361 for ( i = 2; i <= (UINT4)( order ); i++ )
362 params.ppn->data[i] = 1.0;
363
364 /* Output parameters. */
365 memset( &waveform, 0, sizeof(CoherentGW) );
366
367
368 /*******************************************************************
369 * OUTPUT GENERATION *
370 *******************************************************************/
371
372 /* Generate waveform. */
373 SUB( LALGeneratePPNInspiral( &stat, &waveform, &params ), &stat );
374
375 /* Print termination information. */
376 snprintf( message, MSGLENGTH, "%d: %s", params.termCode,
377 params.termDescription );
378 INFO( message );
379
380 /* Print coalescence phase.
381 snprintf( message, MSGLENGTH,
382 "Waveform ends %.3f cycles before coalescence",
383 -waveform.phi->data->data[waveform.phi->data->length-1]
384 / (REAL4)( LAL_TWOPI ) ); */
385 {
386 INT4 code = sprintf( message,
387 "Waveform ends %.3f cycles before coalescence",
388 -waveform.phi->data->data[waveform.phi->data->length
389 -1]
390 / (REAL4)( LAL_TWOPI ) );
391 if ( code >= MSGLENGTH || code < 0 ) {
393 GENERATEPPNINSPIRALTESTC_MSGEPRINT, 0 );
395 }
396 }
397 INFO( message );
398
399 /* Check if sampling interval was too large. */
400 if ( params.dfdt > 2.0 ) {
401 snprintf( message, MSGLENGTH,
402 "Waveform sampling interval is too large:\n"
403 "\tmaximum df*dt = %f", params.dfdt );
404 WARNING( message );
405 }
406
407 /* Renormalize phase. */
408 phii -= waveform.phi->data->data[0];
409 for ( i = 0; i < waveform.phi->data->length; i++ )
410 waveform.phi->data->data[i] += phii;
411
412 /* Write output. */
413 if ( outfile ) {
414 if ( ( fp = fopen( outfile, "w" ) ) == NULL ) {
416 GENERATEPPNINSPIRALTESTC_MSGEFILE, outfile );
418 }
419
420 /* Amplitude and phase functions: */
421 if ( deltat == 0.0 ) {
422 REAL8 t = 0.0; /* time */
423 for ( i = 0; i < waveform.a->data->length; i++, t += dt )
424 fprintf( fp, "%f %.3f %10.3e %10.3e %10.3e\n", t,
425 waveform.phi->data->data[i],
426 waveform.f->data->data[i],
427 waveform.a->data->data[2*i],
428 waveform.a->data->data[2*i+1] );
429 }
430
431 /* Waveform: */
432 else {
433 REAL8 t = 0.0;
434 REAL8 x = 0.0;
435 REAL8 dx = deltat/dt;
436 REAL8 xMax = waveform.a->data->length - 1;
437 REAL8 *phiData = waveform.phi->data->data;
438 REAL4 *fData = waveform.f->data->data;
439 REAL4 *aData = waveform.a->data->data;
440 for ( ; x < xMax; x += dx, t += deltat ) {
441 UINT4 j = floor( x );
442 REAL8 frac = x - j;
443 REAL8 p = frac*phiData[j+1] + ( 1.0 - frac )*phiData[j];
444 REAL8 f = frac*fData[j+1] + ( 1.0 - frac )*fData[j];
445 REAL8 ap = frac*aData[2*j+2] + ( 1.0 - frac )*aData[2*j];
446 REAL8 ac = frac*aData[2*j+3] + ( 1.0 - frac )*aData[2*j+1];
447
448 fprintf( fp, "%f %.3f %10.3e %10.3e %10.3e\n", t, p, f,
449 ap*cos( p ), ac*sin( p ) );
450 fflush( fp );
451 }
452 }
453
454 fclose( fp );
455 }
456
457
458 /*******************************************************************
459 * CLEANUP *
460 *******************************************************************/
461
462 SUB( LALSDestroyVector( &stat, &(params.ppn) ), &stat );
463 SUB( LALSDestroyVectorSequence( &stat, &(waveform.a->data) ),
464 &stat );
465 SUB( LALSDestroyVector( &stat, &(waveform.f->data) ), &stat );
466 SUB( LALDDestroyVector( &stat, &(waveform.phi->data) ), &stat );
467 LALFree( waveform.a );
468 LALFree( waveform.f );
469 LALFree( waveform.phi );
470
472 INFO( GENERATEPPNINSPIRALTESTC_MSGENORM );
474}
475/** \endcond */
#define SUB(func, statusptr)
#define ERROR(code, msg, statement)
#define INFO(statement)
#define WARNING(statement)
#define GENERATEPPNINSPIRALTESTC_EFILE
Could not open file.
#define GENERATEPPNINSPIRALTESTC_EARG
Error parsing arguments.
#define GENERATEPPNINSPIRALTESTC_EPRINT
Wrote past end of message string.
#define GENERATEPPNINSPIRALTESTC_ENORM
Normal exit.
#define USAGE
void LALCheckMemoryLeaks(void)
#define LALFree(p)
#define fprintf
double i
void LALSDestroyVectorSequence(LALStatus *status, REAL4VectorSequence **vectorSequence)
void LALGeneratePPNInspiral(LALStatus *stat, CoherentGW *output, PPNParamStruc *params)
Computes a parametrized post-Newtonian inspiral waveform.
#define LAL_REAL4_MAX
#define LAL_PI
#define LAL_TWOPI
#define LAL_PC_SI
#define LAL_REAL4_MIN
double REAL8
char CHAR
uint32_t UINT4
int32_t INT4
float REAL4
int LALPrintError(const char *fmt,...)
COORDINATESYSTEM_EQUATORIAL
void LALDDestroyVector(LALStatus *, REAL8Vector **)
void LALSDestroyVector(LALStatus *, REAL4Vector **)
void LALSCreateVector(LALStatus *, REAL4Vector **, UINT4)
LIGOTimeGPS * XLALINT8NSToGPS(LIGOTimeGPS *epoch, INT8 ns)
p
x
REAL4TimeVectorSeries * a
REAL8TimeSeries * phi
REAL4TimeSeries * f
This structure stores the parameters for constructing a restricted post-Newtonian waveform.
REAL4Sequence * data
REAL4VectorSequence * data
REAL4 * data
REAL8Sequence * data
REAL8 * data
LIGOTimeGPS epoch
double psi
FILE * fp
double f_min
double f_max
int main(int argc, char **argv)
Definition: version.c:32
char * outfile