LAL: StackMetricTest.c Source File

00001 /*
00002 *  Copyright (C) 2007 Jolien Creighton, Teviet Creighton
00003 *
00004 *  This program is free software; you can redistribute it and/or modify
00005 *  it under the terms of the GNU General Public License as published by
00006 *  the Free Software Foundation; either version 2 of the License, or
00007 *  (at your option) any later version.
00008 *
00009 *  This program is distributed in the hope that it will be useful,
00010 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 *  GNU General Public License for more details.
00013 *
00014 *  You should have received a copy of the GNU General Public License
00015 *  along with with program; see the file COPYING. If not, write to the
00016 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
00017 *  MA  02111-1307  USA
00018 */
00019 
00020 /****************************** <lalVerbatim file="StackMetricTestCV">
00021 Author: Creighton, T. D.
00022 $Id: StackMetricTest.c,v 1.4 2007/06/08 14:41:52 bema Exp $
00023 **************************************************** </lalVerbatim> */
00024 
00025 /********************************************************** <lalLaTeX>
00026 
00027 \subsection{Program \texttt{StackMetricTest.c}}
00028 \label{ss:StackMetricTest.c}
00029 
00030 Computes the parameter space metric for a coherent pulsar search.
00031 
00032 \subsubsection*{Usage}
00033 \begin{verbatim}
00034 StackMetricTest [-p n dt t0] [-l lat lon] [-d debuglevel]
00035                 [ra dec f0 [f1 [...]]]
00036 \end{verbatim}
00037 
00038 \subsubsection*{Description}
00039 
00040 This test program computes the stack search metric for a particular
00041 location in parameter space.  The following option flags are accepted:
00042 \begin{itemize}
00043 \item[\texttt{-p}] Sets the search parameters: the number of stacks
00044 \verb@n@, the length of each stack \verb@dt@ (in seconds), and the
00045 start time of the first stack \verb@t0@ (in seconds of GPS time).
00046 \item[\texttt{-l}] Sets the detector latitude to \verb@lat@ (in
00047 radians north from the equator) and longitude to \verb@lon@ (in
00048 radians east of the prime meridian).
00049 \item[\texttt{-d}] Sets the debug level to \verb@debuglevel@.
00050 \end{itemize}
00051 Once any of the above options are processed, any remaining
00052 command-line options are taken to be the parameter-space location of
00053 the source: its right ascension \verb@ra@ (in radians), its
00054 declination \verb@dec@ (in radians), its frequency \verb@f0@ (in Hz),
00055 and zero or more spindown parameters \verb@f@$k$ (in $\mathrm{Hz}^k$),
00056 all evaluated at the start time of the search \verb@t0@.  If any (or
00057 all) of the command-line arguments are missing, they will be set from
00058 \verb@#define@d defaults.
00059 
00060 \subsubsection*{Exit codes}
00061 ****************************************** </lalLaTeX><lalErrTable> */
00062 #define STACKMETRICTESTC_ENORM 0
00063 #define STACKMETRICTESTC_ESUB  1
00064 #define STACKMETRICTESTC_EARG  2
00065 #define STACKMETRICTESTC_EVAL  3
00066 
00067 #define STACKMETRICTESTC_MSGENORM "Normal exit"
00068 #define STACKMETRICTESTC_MSGESUB  "Subroutine failed"
00069 #define STACKMETRICTESTC_MSGEARG  "Error parsing arguments"
00070 #define STACKMETRICTESTC_MSGEVAL  "Input argument out of valid range"
00071 /******************************************** </lalErrTable><lalLaTeX>
00072 
00073 \subsubsection*{Algorithm}
00074 
00075 \subsubsection*{Uses}
00076 \begin{verbatim}
00077 lalDebugLevel
00078 LALPrintError()                 LALCheckMemoryLeaks()
00079 LALMalloc()                     LALFree()
00080 LALDCreateVector()              LALDDestroyVector()
00081 LALStackMetric()                LALProjectMetric()
00082 LALTBaryPtolemaic()             LALDTBaryPtolemaic()
00083 LALTSpin()                      LALDTSpin()
00084 LALDTComp()
00085 LALGetEarthTimes()
00086 \end{verbatim}
00087 
00088 \subsubsection*{Notes}
00089 
00090 \vfill{\footnotesize\input{StackMetricTestCV}}
00091 
00092 ******************************************************* </lalLaTeX> */
00093 
00094 #include <math.h>
00095 #include <stdlib.h>
00096 #include <lal/LALStdlib.h>
00097 #include <lal/LALConstants.h>
00098 #include <lal/AVFactories.h>
00099 #include <lal/StackMetric.h>
00100 #include <lal/PulsarTimes.h>
00101 
00102 NRCSID(STACKMETRICTESTC,"$Id: StackMetricTest.c,v 1.4 2007/06/08 14:41:52 bema Exp $");
00103 
00104 /* Default parameter settings. */
00105 int lalDebugLevel=0;
00106 #define NSTACKS 1
00107 #define STACKLENGTH 100000.0  /* arbitrary */
00108 #define STARTTIME 0.0         /* arbitrary */
00109 #define LATITUDE 0.91188      /* GEO600 location */
00110 #define LONGITUDE 0.17142     /* GEO600 location */
00111 #define RIGHTASCENSION 5.0309 /* a known pulsar */
00112 #define DECLINATION 0.27925   /* a known pulsar */
00113 #define FREQUENCY 1000.0      /* arbitrary */
00114 
00115 /* Usage format string. */
00116 #define USAGE "Usage: %s [-p n dt t0] [-l lat lon] [-d debuglevel]\n\
00117                        [ra dec f0 [f1 [...]]]\n"
00118 
00119 /* Input error checking: Some accepted parameter ranges. */
00120 #define NMAX  10000 /* 1 <= number of stacks <= NMAX */
00121 #define DTMAX  3e10 /* 1/f_0 < stack length <= DTMAX */
00122 #define F0MAX  1e4  /* 0 < f_0 <= FOMAX */
00123 #define TAUMIN 1e4  /* |f_k| <= TAUMIN^(-k) */
00124 /* |latitude| and |declination| are <= LAL_PI,
00125    |longitude| and |right ascension| are <= LAL_TWOPI */
00126 
00127 /* Macros for printing errors and testing subroutines. */
00128 #define ERROR( code, msg, statement )                                \
00129 do                                                                   \
00130 if ( lalDebugLevel & LALERROR )                                      \
00131 {                                                                    \
00132   LALPrintError( "Error[0] %d: program %s, file %s, line %d, %s\n"   \
00133                  "        %s %s\n", (code), *argv, __FILE__,         \
00134                  __LINE__, STACKMETRICTESTC, statement ? statement : \
00135                  "", (msg) );                                        \
00136 }                                                                    \
00137 while (0)
00138 
00139 #define INFO( statement )                                            \
00140 do                                                                   \
00141 if ( lalDebugLevel & LALINFO )                                       \
00142 {                                                                    \
00143   LALPrintError( "Info[0]: program %s, file %s, line %d, %s\n"       \
00144                  "        %s\n", *argv, __FILE__, __LINE__,          \
00145                  STACKMETRICTESTC, (statement) );                    \
00146 }                                                                    \
00147 while (0)
00148 
00149 #define SUB( func, statusptr )                                       \
00150 do                                                                   \
00151 if ( (func), (statusptr)->statusCode )                               \
00152 {                                                                    \
00153   ERROR( STACKMETRICTESTC_ESUB, STACKMETRICTESTC_MSGESUB,            \
00154          "Function call \"" #func "\" failed:" );                    \
00155   return STACKMETRICTESTC_ESUB;                                      \
00156 }                                                                    \
00157 while (0)
00158 
00159 #define CHECKVAL( val, lower, upper )                                \
00160 do                                                                   \
00161 if ( ( (val) <= (lower) ) || ( (val) > (upper) ) )                   \
00162 {                                                                    \
00163   ERROR( STACKMETRICTESTC_EVAL, STACKMETRICTESTC_MSGESUB,            \
00164          "Value of " #val " out of range:" );                        \
00165   LALPrintError( #val " = %f, range = (%f,%f]\n", (REAL8)(val),      \
00166                  (REAL8)(lower), (REAL8)(upper) );                   \
00167   return STACKMETRICTESTC_EVAL;                                      \
00168 }                                                                    \
00169 while (0)
00170 
00171 /* A global pointer for debugging. */
00172 #ifndef NDEBUG
00173 char *lalWatch;
00174 #endif
00175 
00176 #include <lal/LALStdio.h>
00177 #define MAXLEN 1024
00178 int
00179 fprintderr( FILE *fp, REAL8 x, REAL8 dx );
00180 
00181 int
00182 main(int argc, char **argv)
00183 {
00184   static LALStatus stat;
00185   INT4  arg;
00186   UINT4 i, j, k;
00187   UINT4 nSpin = 0;
00188   UINT4 nSky = 2;
00189   UINT4 n = NSTACKS;
00190   REAL8 dt = STACKLENGTH;
00191   REAL8 t0 = STARTTIME;
00192   REAL8 lat = LATITUDE;
00193   REAL8 lon = LONGITUDE;
00194   REAL8 ra = RIGHTASCENSION;
00195   REAL8 dec = DECLINATION;
00196   REAL8 f0 = FREQUENCY;
00197   REAL8 *spin = NULL;
00198   REAL8Vector *lambda = NULL;
00199   REAL8Vector *metric = NULL;
00200   static LIGOTimeGPS start;
00201   static PulsarTimesParamStruc spinParams;
00202   static PulsarTimesParamStruc baryParams;
00203   static PulsarTimesParamStruc compParams;
00204   static MetricParamStruc params;
00205 
00206   /* Parse argument list.  arg stores the current position. */
00207   arg = 1;
00208   while ( arg < argc ) {
00209     /* Parse search parameter option. */
00210     if ( !strcmp( argv[arg], "-p" ) ) {
00211       if ( argc > arg + 3 ) {
00212         arg++;
00213         n = atoi( argv[arg++] );
00214         dt = atof( argv[arg++] );
00215         t0 = atof( argv[arg++] );
00216       }else{
00217         ERROR( STACKMETRICTESTC_EARG, STACKMETRICTESTC_MSGEARG, 0 );
00218         LALPrintError( USAGE, *argv );
00219         return STACKMETRICTESTC_EARG;
00220       }
00221     }
00222     /* Parse detector position option. */
00223     else if ( !strcmp( argv[arg], "-l" ) ) {
00224       if ( argc > arg + 2 ) {
00225         arg++;
00226         lat = atof( argv[arg++] );
00227         lon = atof( argv[arg++] );
00228       }else{
00229         ERROR( STACKMETRICTESTC_EARG, STACKMETRICTESTC_MSGEARG, 0 );
00230         LALPrintError( USAGE, *argv );
00231         return STACKMETRICTESTC_EARG;
00232       }
00233     }
00234     /* Parse debug level option. */
00235     else if ( !strcmp( argv[arg], "-d" ) ) {
00236       if ( argc > arg + 1 ) {
00237         arg++;
00238         lalDebugLevel = atoi( argv[arg++] );
00239       }else{
00240         ERROR( STACKMETRICTESTC_EARG, STACKMETRICTESTC_MSGEARG, 0 );
00241         LALPrintError( USAGE, *argv );
00242         return STACKMETRICTESTC_EARG;
00243       }
00244     }
00245     /* Check for unrecognized options. */
00246     else if ( argv[arg][0] == '-' ) {
00247       ERROR( STACKMETRICTESTC_EARG, STACKMETRICTESTC_MSGEARG, 0 );
00248       LALPrintError( USAGE, *argv );
00249       return STACKMETRICTESTC_EARG;
00250     }
00251     /* Parse remaining options. */
00252     else {
00253       if ( argc > arg + 2 ) {
00254         ra = atof( argv[arg++] );
00255         dec = atof( argv[arg++] );
00256         f0 = atof( argv[arg++] );
00257         nSpin = argc - arg;
00258         if ( nSpin )
00259           spin = (REAL8 *)LALMalloc( nSpin*sizeof(REAL8) );
00260         j = 0;
00261         while ( arg < argc )
00262           spin[j++] = atof( argv[arg++] );
00263       }else{
00264         ERROR( STACKMETRICTESTC_EARG, STACKMETRICTESTC_MSGEARG, 0 );
00265         LALPrintError( USAGE, *argv );
00266         return STACKMETRICTESTC_EARG;
00267       }
00268     }
00269   } /* End of argument parsing loop. */
00270 
00271   /* Do error trapping on input parameters. */
00272   if ( lalDebugLevel & LALERROR ) {
00273     CHECKVAL( n, 0, NMAX );
00274     CHECKVAL( dt, 1.0/f0, DTMAX );
00275     CHECKVAL( lat, -LAL_PI, LAL_PI );
00276     CHECKVAL( lon, -LAL_TWOPI, LAL_TWOPI );
00277     CHECKVAL( ra, -LAL_TWOPI, LAL_TWOPI );
00278     CHECKVAL( dec, -LAL_PI, LAL_PI );
00279     CHECKVAL( f0, 0.0, F0MAX );
00280     for ( j = 0; j < nSpin; j++ ) {
00281       REAL4 inverseAge = pow( spin[j], 1.0/( j + 1 ) );
00282       CHECKVAL( inverseAge, -1.0/TAUMIN, 1.0/TAUMIN );
00283     }
00284   }
00285 
00286   /* Set up start time. */
00287   start.gpsSeconds = (INT4)t0;
00288   start.gpsNanoSeconds = (INT4)( 1.0e9*( t0 - start.gpsSeconds ) );
00289   t0 = 0.0;
00290 
00291   /* Set up constant parameters for barycentre transformation. */
00292   baryParams.epoch = start;
00293   baryParams.latitude = lat;
00294   baryParams.longitude = lon;
00295   SUB( LALGetEarthTimes( &stat, &baryParams ), &stat );
00296 
00297   /* Set up constant parameters for spindown transformation. */
00298   spinParams.epoch = start;
00299   spinParams.t0 = t0;
00300 
00301   /* Set up constant parameters for composed transformation. */
00302   compParams.epoch = start;
00303   compParams.t1 = LALTBaryPtolemaic;
00304   compParams.t2 = LALTSpin;
00305   compParams.dt1 = LALDTBaryPtolemaic;
00306   compParams.dt2 = LALDTSpin;
00307   compParams.constants1 = &baryParams;
00308   compParams.constants2 = &spinParams;
00309   compParams.nArgs = 2;
00310 
00311   /* Set up constant parameters for metric calculation. */
00312   params.dtCanon = LALDTComp;
00313   params.constants = &compParams;
00314   /* To ignore spindown, uncomment the following:
00315   params.dtCanon = LALDTBaryPtolemaic;
00316   params.constants = &baryParams;
00317   nSpin = 0; */
00318   /* To ignore detector motion, uncomment the following:
00319   params.dtCanon = LALDTSpin;
00320   params.constants = &spinParams;
00321   nSky = 0; */
00322 
00323   params.start = t0;
00324   params.deltaT = dt;
00325   params.n = n;
00326   params.errors = 1;
00327 
00328   /* Set up the parameter list. */
00329   SUB( LALDCreateVector( &stat, &lambda, nSpin + nSky + 1 ), &stat );
00330   lambda->data[0] = f0;
00331   if ( nSky ) {
00332     lambda->data[1] = ra;
00333     lambda->data[2] = dec;
00334   }
00335   if ( nSpin ) {
00336     memcpy( lambda->data + nSky + 1, spin, nSpin*sizeof(REAL8) );
00337     LALFree( spin );
00338   }
00339 
00340   /* Compute the metric, and free lambda. */
00341   n = nSpin + nSky + 1;
00342   if ( params.errors )
00343     SUB( LALDCreateVector( &stat, &metric, n*(n+1) ), &stat );
00344   else
00345     SUB( LALDCreateVector( &stat, &metric, (n*(n+1)) >> 1 ), &stat );
00346   SUB( LALStackMetric( &stat, metric, lambda, &params ), &stat );
00347   SUB( LALDDestroyVector( &stat, &lambda ), &stat );
00348 
00349   /* Print the metric, project it, and print it again. */
00350   for ( i=0, k=0; i<n; i++ ){
00351     if ( params.errors )
00352       for ( j=0; j<=i; j++, k+=2 )
00353         fprintf( stdout, "%8.1e(%7.1e) ", metric->data[k],
00354                  metric->data[k+1] );
00355     else
00356       for ( j=0; j<=i; j++, k++ )
00357         fprintf( stdout, "%10.3e ", metric->data[k] );
00358     fprintf( stdout, "\n" );
00359   }
00360   SUB( LALProjectMetric( &stat, metric, params.errors ), &stat );
00361   fprintf( stdout, "\n" );
00362   for ( i=0, k=0; i<n; i++ ){
00363     if ( params.errors )
00364       for ( j=0; j<=i; j++, k+=2 )
00365         fprintf( stdout, "%8.1e(%7.1e) ", metric->data[k],
00366                  metric->data[k+1] );
00367     else
00368       for ( j=0; j<=i; j++, k++ )
00369         fprintf( stdout, "%10.3e ", metric->data[k] );
00370     fprintf( stdout, "\n" );
00371   }
00372 
00373   if ( nSpin == 0 ) {
00374     if ( params.errors ) {
00375       REAL8 det = metric->data[4]*metric->data[10]
00376         - metric->data[8]*metric->data[8];
00377       REAL8 ddet = fabs( metric->data[4]*metric->data[11] )
00378         + fabs( metric->data[5]*metric->data[10] )
00379         + fabs( metric->data[8]*metric->data[9] )*2.0;
00380       fprintderr( stdout, det, ddet );
00381       fprintf( stdout, "\n" );
00382     } else {
00383       REAL8 det = metric->data[2]*metric->data[5]
00384         - metric->data[4]*metric->data[4];
00385       fprintf( stdout, "%10.3e\n", det );
00386     }
00387   }
00388 
00389   /* Free the metric, and exit. */
00390   SUB( LALDDestroyVector( &stat, &metric ), &stat );
00391   LALCheckMemoryLeaks();
00392   INFO( STACKMETRICTESTC_MSGENORM );
00393   return STACKMETRICTESTC_ENORM;
00394 }
00395 
00396 
00397 int
00398 fprintderr( FILE *fp, REAL8 x, REAL8 dx ) {
00399   CHAR format[MAXLEN]; /* format string for fprintf() */
00400   INT4 gsd = 0;        /* place index of greatest significant digit */
00401   INT4 lsd = 0;        /* place index of least significant digit */
00402   REAL8 norm;          /* normalization factor */
00403 
00404   /* Compute gsd, lsd, y, and dy. */
00405   if ( dx < LAL_REAL8_EPS*fabs( x ) )
00406     dx = 0.0;
00407   if ( dx > 0.0 ) {
00408     REAL8 lsdd = log( 0.5*dx )/log( 10.0 );
00409     if ( lsdd >= 0.0 )
00410       lsd = (INT4)( lsdd );
00411     else
00412       lsd = (INT4)( lsdd ) - 1;
00413   }
00414   if ( x != 0.0 ) {
00415     REAL8 gsdd = log( fabs( x ) )/log( 10.0 );
00416     if ( gsdd >= 0.0 )
00417       gsd = (INT4)( gsdd );
00418     else
00419       gsd = (INT4)( gsdd ) - 1;
00420   }
00421 
00422   /* If x is zero, format is determined entirely by dx. */
00423   if ( x == 0.0 ) {
00424     if ( dx <= 0.0 )
00425       return fprintf( fp, "0" );
00426     if ( abs( lsd ) > 3 ) {
00427       norm = pow( 10.0, -lsd );
00428       return fprintf( fp, "( 0 +/- %.0f )e%+i", dx*norm, lsd );
00429     }
00430     if ( lsd <= 0 ) {
00431       LALSnprintf( format, MAXLEN, "%%.%if +/- %%.%if", -lsd, -lsd );
00432       return fprintf( fp, format, 0.0, dx );
00433     }
00434     norm = pow( 10.0, -lsd );
00435     LALSnprintf( format, MAXLEN, "0 +/- %%.0f%%0%ii", lsd );
00436     return fprintf( fp, format, dx*norm, 0 );
00437   }
00438 
00439   /* If number is exact to 8-byte precision, print it as such. */
00440   if ( dx <= 0.0 ) {
00441     if ( abs( gsd ) > 3 )
00442       return fprintf( fp, "%.16e", x );
00443     LALSnprintf( format, MAXLEN, "%%.%if", 16 - gsd );
00444     return fprintf( fp, format, x );
00445   }
00446 
00447   /* Otherwise, format depends on x and dx. */
00448   if ( gsd < lsd )
00449     gsd = lsd;
00450   if ( lsd > 3 || gsd < -3 ) {
00451     norm = pow( 10.0, -gsd );
00452     LALSnprintf( format, MAXLEN, "( %%.%if +/- %%.%if )e%+i",
00453                  gsd - lsd, gsd - lsd, gsd );
00454     return fprintf( fp, format, x*norm, dx*norm );
00455   }
00456   if ( lsd <= 0 ) {
00457     LALSnprintf( format, MAXLEN, "%%.%if +/- %%.%if", -lsd, -lsd );
00458     return fprintf( fp, format, x, dx );
00459   }
00460   norm = pow( 10.0, -lsd );
00461   LALSnprintf( format, MAXLEN, "%%.0f%%0%ii +/- %%.0f%%0%ii", lsd,
00462                lsd );
00463   return fprintf( fp, format, x*norm, 0, dx*norm, 0 );
00464 }