Commit 5156aabe authored by Reinhard Prix's avatar Reinhard Prix

BinarySSBTimesTest: test for XLALAdd[Multi]BinaryTimes() functions

  - simple comparison in a "ScoX1-like" binary test-example
    between old LALGet[Multi]Binarytimes() functions and the new
    XLALAdd[Multi]BinaryTimes()
  - the old LAL functions have been moved into the test-file, as they
    are not used anywhere else any more.
  - refs #669, #670
Original: 7893ca0766b9db689ca64eca2f9d6968ebb34b17
parent cb722609
/*
* Copyright (C) 2013 Reinhard Prix
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
#include <config.h>
#include <math.h>
#include <sys/times.h>
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif
#ifdef HAVE_GETOPT_H
#include <getopt.h>
#endif
#include <lal/ComputeFstat.h>
#include <lal/LALInitBarycenter.h>
#include <lal/FindRoot.h>
/** \author Reinhard Prix
* \file
* \ingroup ComputeFstat_h
* \brief Tests for XLALAdd[Multi]BinaryTimes()
*
* We simply compare the results to the old+obsolete LAL functions LALGet[Multi]Binarytimes(),
* which have been moved here, and only serve for this comparison.
*
* Sky-location is picked at random each time, which allows a minimal
* Monte-Carlo validation by simply running this script repeatedly.
*
*/
// ----- local defines
#define EA_ACC 1E-9 /* the timing accuracy of LALGetBinaryTimes in seconds */
// ----- macros
#define GPS2REAL8(gps) (1.0 * (gps).gpsSeconds + 1.e-9 * (gps).gpsNanoSeconds )
// ----- global variables
extern int lalDebugLevel;
static REAL8 p,q,r; /* binary time delay coefficients (need to be global so that the LAL root finding procedure can see them) */
// initializers
LALStatus empty_LALStatus;
/* ----- internal prototypes ---------- */
static void LALGetBinarytimes (LALStatus *, SSBtimes *tBinary, const SSBtimes *tSSB, const DetectorStateSeries *DetectorStates, const BinaryOrbitParams *binaryparams, LIGOTimeGPS refTime);
static void LALGetMultiBinarytimes (LALStatus *, MultiSSBtimes **multiBinary, const MultiSSBtimes *multiSSB, const MultiDetectorStateSeries *multiDetStates, const BinaryOrbitParams *binaryparams, LIGOTimeGPS refTime);
static void EccentricAnomoly(LALStatus *status, REAL8 *tr, REAL8 lE, void *tr0);
int XLALCompareSSBtimes ( REAL8 *err_DeltaT, REAL8 *err_Tdot, const SSBtimes *t1, const SSBtimes *t2 );
int XLALCompareMultiSSBtimes ( REAL8 *err_DeltaT, REAL8 *err_Tdot, const MultiSSBtimes *m1, const MultiSSBtimes *m2 );
/* ----- function definitions ---------- */
int
main ( int argc, char *argv[] )
{
LALStatus status = empty_LALStatus;
/* read user input */
lalDebugLevel = 0;
int opt;
while ((opt = getopt( argc, argv, "v:" )) != -1) {
switch (opt) {
case 'v': /* set lalDebugLevel */
lalDebugLevel = atoi( optarg );
break;
default:
XLAL_ERROR ( XLAL_EINVAL, "Invalid commandline-option '%c'\n", opt );
break;
} // switch
}
REAL8 startTimeREAL8 = 714180733;
REAL8 duration = 180000; /* 50 hours */
REAL8 Tsft = 1800; /* assume 30min SFTs */
char earthEphem[] = DATADIR "earth00-19-DE200.dat.gz";
char sunEphem[] = DATADIR "sun00-19-DE200.dat.gz";
//REAL8 tolerance = 2e-10; /* same algorithm, should be basically identical results */
LIGOTimeGPS startTime, refTime;
XLALGPSSetREAL8 ( &startTime, startTimeREAL8 );
refTime = startTime;
// init random-generator
struct tms buf;
UINT4 seed = times(&buf);
srand ( seed );
XLALPrintInfo ("seed used = %d\n", seed );
// pick skyposition at random ----- */
SkyPosition skypos;
skypos.longitude = LAL_TWOPI * (1.0 * rand() / ( RAND_MAX + 1.0 ) ); // alpha uniform in [0, 2pi)
skypos.latitude = LAL_PI_2 - acos ( 1 - 2.0 * rand()/RAND_MAX ); // sin(delta) uniform in [-1,1]
skypos.system = COORDINATESYSTEM_EQUATORIAL;
// pick binary orbital parameters:
// somewhat inspired by Sco-X1 parameters from S2-paper (PRD76, 082001 (2007), gr-qc/0605028)
// but with a more extreme eccentricity, and random argp
REAL8 argp = LAL_TWOPI * (1.0 * rand() / ( RAND_MAX + 1.0 ) ); // uniform in [0, 2pi)
BinaryOrbitParams orbit;
XLALGPSSetREAL8 ( &orbit.tp, 731163327 ); // time of observed periapsis passage (in SSB)
orbit.argp = argp; // argument of periapsis (radians)
orbit.asini = 1.44; // projected, normalized orbital semi-major axis (s) */
orbit.ecc = 1e-2; // relatively large value, for better testing
orbit.period = 68023; // period (s) : about ~18.9h
// ----- step 0: prepare test-case input for calling the BinarySSB-functions
// setup detectors
const char *sites[3] = { "H1", "L1", "V1" };
UINT4 numDetectors = sizeof( sites ) / sizeof ( sites[0] );
MultiLALDetector *multiIFO = XLALCreateMultiLALDetector ( numDetectors );
XLAL_CHECK ( multiIFO != NULL, XLAL_EFUNC, "XLALCreateMultiLALDetector(%d) failed.\n", numDetectors );
for ( UINT4 X = 0; X < numDetectors; X ++ )
{
LALDetector *det = XLALGetSiteInfo ( sites[X] );
XLAL_CHECK ( det != NULL, XLAL_EFUNC, "XLALGetSiteInfo ('%s') failed for detector X=%d\n", sites[X], X );
multiIFO->data[X] = (*det); // struct copy
XLALFree ( det );
}
// load ephemeris
EphemerisData *edat = XLALInitBarycenter ( earthEphem, sunEphem );
XLAL_CHECK ( edat != NULL, XLAL_EFUNC, "XLALInitBarycenter('%s','%s') failed\n", earthEphem, sunEphem );
// setup multi-timeseries
MultiLIGOTimeGPSVector *multiTS;
XLAL_CHECK ( (multiTS = XLALCalloc ( 1, sizeof(*multiTS))) != NULL, XLAL_ENOMEM );
XLAL_CHECK ( (multiTS->data = XLALCalloc (numDetectors, sizeof(*multiTS->data))) != NULL, XLAL_ENOMEM );
multiTS->length = numDetectors;
for ( UINT4 X = 0; X < numDetectors; X ++ )
{
multiTS->data[X] = XLALMakeTimestamps ( startTime, duration, Tsft );
XLAL_CHECK ( multiTS->data[X] != NULL, XLAL_EFUNC, "XLALMakeTimestamps() failed.\n");
} /* for X < numIFOs */
// generate detector-states
MultiDetectorStateSeries *multiDetStates = XLALGetMultiDetectorStates ( multiTS, multiIFO, edat, 0 );
XLAL_CHECK ( multiDetStates != NULL, XLAL_EFUNC, "XLALGetMultiDetectorStates() failed.\n");
// generate isolated-NS SSB times
MultiSSBtimes *multiSSBIn = XLALGetMultiSSBtimes ( multiDetStates, skypos, refTime, SSBPREC_RELATIVISTICOPT );
XLAL_CHECK ( multiSSBIn != NULL, XLAL_EFUNC, "XLALGetMultiSSBtimes() failed.\n");
// ----- step 1: compute reference-result using old LALGetMultiBinarytimes()
MultiSSBtimes *multiBinary_ref = NULL;
LALGetMultiBinarytimes (&status, &(multiBinary_ref), multiSSBIn, multiDetStates, &orbit, refTime );
XLAL_CHECK ( status.statusCode == 0, XLAL_EFAILED, "LALGetMultiBinarytimes() failed with status = %d : '%s'\n", status.statusCode, status.statusDescription );
// ----- step 2: compute test-result using new XLALAddMultiBinaryTimes()
MultiSSBtimes *multiBinary_test = NULL;
XLAL_CHECK ( XLALAddMultiBinaryTimes ( &multiBinary_test, multiSSBIn, &orbit ) == XLAL_SUCCESS, XLAL_EFUNC );
// ----- step 3: compare results
REAL8 err_DeltaT, err_Tdot;
REAL8 tolerance = 1e-10;
int ret = XLALCompareMultiSSBtimes ( &err_DeltaT, &err_Tdot, multiBinary_ref, multiBinary_test );
XLAL_CHECK ( ret == XLAL_SUCCESS, XLAL_EFUNC, "XLALCompareMultiSSBtimes() failed.\n");
XLALPrintWarning ( "INFO: err(DeltaT) = %g, err(Tdot) = %g\n", err_DeltaT, err_Tdot );
XLAL_CHECK ( err_DeltaT < tolerance, XLAL_ETOL, "error(DeltaT) = %g exceeds tolerance of %g\n", err_DeltaT, tolerance );
XLAL_CHECK ( err_Tdot < tolerance, XLAL_ETOL, "error(Tdot) = %g exceeds tolerance of %g\n", err_Tdot, tolerance );
// ---- step 4: clean-up memory
XLALDestroyEphemerisData ( edat );
XLALDestroyMultiSSBtimes ( multiBinary_test );
XLALDestroyMultiSSBtimes ( multiBinary_ref );
XLALDestroyMultiSSBtimes ( multiSSBIn );
XLALDestroyMultiTimestamps ( multiTS );
XLALDestroyMultiDetectorStateSeries ( multiDetStates );
XLALDestroyMultiLALDetector ( multiIFO );
// check for memory-leaks
LALCheckMemoryLeaks();
return XLAL_SUCCESS;
} // main()
/**
* Compare 2 MultiSSBtimes vectors and return the maximal difference in
* DeltaT in s, and in Tdot (dimensionless)
*/
int
XLALCompareMultiSSBtimes ( REAL8 *err_DeltaT, REAL8 *err_Tdot, const MultiSSBtimes *m1, const MultiSSBtimes *m2 )
{
XLAL_CHECK ( err_DeltaT != NULL, XLAL_EINVAL );
XLAL_CHECK ( err_Tdot != NULL, XLAL_EINVAL );
XLAL_CHECK ( m1 != NULL, XLAL_EINVAL );
XLAL_CHECK ( m2 != NULL, XLAL_EINVAL );
XLAL_CHECK ( m1->length == m2->length, XLAL_EINVAL, "Number of detectors differ %d != %d\n", m1->length, m2->length );
REAL8 max_DeltaT = 0;
REAL8 max_Tdot = 0;
UINT4 numDet = m1->length;
for ( UINT4 X = 0; X < numDet; X ++ )
{
REAL8 DeltaT_X, Tdot_X;
XLAL_CHECK ( XLALCompareSSBtimes ( &DeltaT_X, &Tdot_X, m1->data[X], m2->data[X] ) == XLAL_SUCCESS, XLAL_EFUNC );
max_DeltaT = fmax ( max_DeltaT, DeltaT_X );
max_Tdot = fmax ( max_Tdot, Tdot_X );
} // for X < numDet
(*err_DeltaT) = max_DeltaT;
(*err_Tdot) = max_Tdot;
return XLAL_SUCCESS;
} // XLALCompareMultiSSBtimes()
/**
* Compare 2 SSBtimes vectors and return the maximal difference in
* DeltaT in s, and in Tdot (dimensionless)
*/
int
XLALCompareSSBtimes ( REAL8 *err_DeltaT, REAL8 *err_Tdot, const SSBtimes *t1, const SSBtimes *t2 )
{
XLAL_CHECK ( err_DeltaT != NULL, XLAL_EINVAL );
XLAL_CHECK ( err_Tdot != NULL, XLAL_EINVAL );
XLAL_CHECK ( t1 != NULL, XLAL_EINVAL );
XLAL_CHECK ( t2 != NULL, XLAL_EINVAL );
XLAL_CHECK ( (t1->DeltaT != NULL) && (t1->Tdot != NULL), XLAL_EINVAL );
XLAL_CHECK ( XLALGPSDiff ( &(t1->refTime), &(t2->refTime) ) == 0, XLAL_ETOL,
"Different reference-times: %f != %f\n", XLALGPSGetREAL8(&(t1->refTime)), XLALGPSGetREAL8(&(t2->refTime)));
UINT4 numSteps = t1->DeltaT->length;
XLAL_CHECK ( t1->Tdot->length == numSteps, XLAL_EINVAL );
XLAL_CHECK ( t2->DeltaT->length == numSteps, XLAL_EINVAL );
XLAL_CHECK ( t2->Tdot->length == numSteps, XLAL_EINVAL );
REAL8 max_DeltaT = 0;
REAL8 max_Tdot = 0;
for ( UINT4 i = 0; i < numSteps; i ++ )
{
REAL8 DeltaT_i = fabs ( t1->DeltaT->data[i] - t2->DeltaT->data[i] );
REAL8 Tdot_i = fabs ( t1->Tdot->data[i] - t2->Tdot->data[i] );
max_DeltaT = fmax ( max_DeltaT, DeltaT_i );
max_Tdot = fmax ( max_Tdot, Tdot_i );
} // for i < numSteps
(*err_DeltaT) = max_DeltaT;
(*err_Tdot) = max_Tdot;
return XLAL_SUCCESS;
} // XLALCompareSSBtimes()
// ---------- obsolete LAL functions LALGet[Multi]Binarytimes() kept here for comparison purposes
/** For a given OrbitalParams, calculate the time-differences
* \f$\Delta T_\alpha\equiv T(t_\alpha) - T_0\f$, and their
* derivatives \f$Tdot_\alpha \equiv d T / d t (t_\alpha)\f$.
*
* \note The return-vectors \a DeltaT and \a Tdot must be allocated already
* and have the same length as the input time-series \a DetStates.
*
*/
void
LALGetBinarytimes (LALStatus *status, /**< pointer to LALStatus structure */
SSBtimes *tBinary, /**< [out] DeltaT_alpha = T(t_alpha) - T_0; and Tdot(t_alpha) */
const SSBtimes *tSSB, /**< [in] DeltaT_alpha = T(t_alpha) - T_0; and Tdot(t_alpha) */
const DetectorStateSeries *DetectorStates, /**< [in] detector-states at timestamps t_i */
const BinaryOrbitParams *binaryparams, /**< [in] source binary orbit parameters */
LIGOTimeGPS refTime /**< SSB reference-time T_0 of pulsar-parameters */
)
{
UINT4 numSteps, i;
REAL8 refTimeREAL8;
REAL8 Porb; /* binary orbital period */
REAL8 asini; /* the projected orbital semimajor axis */
REAL8 e,ome ; /* the eccentricity, one minus eccentricity */
REAL8 sinw,cosw; /* the sin and cos of the argument of periapsis */
REAL8 tSSB_now; /* the SSB time at the midpoint of each SFT in REAL8 form */
REAL8 fracorb; /* the fraction of orbits completed since current SSB time */
REAL8 E; /* the eccentric anomoly */
DFindRootIn input; /* the input structure for the root finding procedure */
REAL8 acc; /* the accuracy in radians of the eccentric anomoly computation */
INITSTATUS(status);
ATTATCHSTATUSPTR (status);
ASSERT (DetectorStates, status, COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
numSteps = DetectorStates->length; /* number of timestamps */
ASSERT (tSSB, status,COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
ASSERT (tSSB->DeltaT, status,COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
ASSERT (tSSB->Tdot, status, COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
ASSERT (tBinary, status,COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
ASSERT (tBinary->DeltaT, status,COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
ASSERT (tBinary->Tdot, status, COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
ASSERT (tSSB->DeltaT->length == numSteps, status, COMPUTEFSTATC_EINPUT, COMPUTEFSTATC_MSGEINPUT);
ASSERT (tSSB->Tdot->length == numSteps, status, COMPUTEFSTATC_EINPUT, COMPUTEFSTATC_MSGEINPUT);
ASSERT (tBinary->DeltaT->length == numSteps, status, COMPUTEFSTATC_EINPUT, COMPUTEFSTATC_MSGEINPUT);
ASSERT (tBinary->Tdot->length == numSteps, status, COMPUTEFSTATC_EINPUT, COMPUTEFSTATC_MSGEINPUT);
/* convenience variables */
Porb = binaryparams->period;
e = binaryparams->ecc;
asini = binaryparams->asini;
sinw = sin(binaryparams->argp);
cosw = cos(binaryparams->argp);
ome = 1.0 - e;
refTimeREAL8 = GPS2REAL8(refTime);
/* compute p, q and r coeeficients */
p = (LAL_TWOPI/Porb)*cosw*asini*sqrt(1.0-e*e);
q = (LAL_TWOPI/Porb)*sinw*asini;
r = (LAL_TWOPI/Porb)*sinw*asini*ome;
/* Calculate the required accuracy for the root finding procedure in the main loop */
acc = LAL_TWOPI*(REAL8)EA_ACC/Porb; /* EA_ACC is defined above and represents the required timing precision in seconds (roughly) */
/* loop over the SFTs */
for (i=0; i < numSteps; i++ )
{
/* define SSB time for the current SFT midpoint */
tSSB_now = refTimeREAL8 + (tSSB->DeltaT->data[i]);
/* define fractional orbit in SSB frame since periapsis (enforce result 0->1) */
/* the result of fmod uses the dividend sign hence the second procedure */
{
REAL8 temp = fmod((tSSB_now - GPS2REAL8(binaryparams->tp)),Porb)/(REAL8)Porb;
fracorb = temp - (REAL8)floor(temp);
}
/* compute eccentric anomaly using a root finding procedure */
input.function = EccentricAnomoly; /* This is the name of the function we must solve to find E */
input.xmin = 0.0; /* We know that E will be found between 0 and 2PI */
input.xmax = LAL_TWOPI;
input.xacc = acc; /* The accuracy of the root finding procedure */
/* expand domain until a root is bracketed */
LALDBracketRoot(status->statusPtr,&input,&fracorb);
/* bisect domain to find eccentric anomoly E corresponding to the SSB time of the midpoint of this SFT */
LALDBisectionFindRoot(status->statusPtr,&E,&input,&fracorb);
/* use our value of E to compute the additional binary time delay */
tBinary->DeltaT->data[i] = tSSB->DeltaT->data[i] - ( asini*sinw*(cos(E)-e) + asini*cosw*sqrt(1.0-e*e)*sin(E) );
/* combine with Tdot (dtSSB_by_dtdet) -> dtbin_by_dtdet */
tBinary->Tdot->data[i] = tSSB->Tdot->data[i] * ( (1.0 - e*cos(E))/(1.0 + p*cos(E) - q*sin(E)) );
} /* for i < numSteps */
DETATCHSTATUSPTR (status);
RETURN(status);
} /* LALGetBinarytimes() */
/** For a given set of binary parameters we solve the following function for
* the eccentric anomoly E
*/
static void EccentricAnomoly(LALStatus *status,
REAL8 *tr,
REAL8 lE,
void *tr0
)
{
INITSTATUS(status);
ASSERT(tr0,status, 1, "Null pointer");
/* this is the function relating the observed time since periapse in the SSB to the true eccentric anomoly E */
*tr = *(REAL8 *)tr0*(-1.0) + (lE + (p*sin(lE)) + q*(cos(lE) - 1.0) + r)/(REAL8)LAL_TWOPI;
RETURN(status);
}
/** Multi-IFO version of LALGetBinarytimes().
* Get all binary-timings for all input detector-series.
*
*/
void
LALGetMultiBinarytimes (LALStatus *status, /**< pointer to LALStatus structure */
MultiSSBtimes **multiBinary, /**< [out] SSB-timings for all input detector-state series */
const MultiSSBtimes *multiSSB, /**< [in] SSB-timings for all input detector-state series */
const MultiDetectorStateSeries *multiDetStates, /**< [in] detector-states at timestamps t_i */
const BinaryOrbitParams *binaryparams, /**< [in] source binary orbit parameters */
LIGOTimeGPS refTime /**< SSB reference-time T_0 for SSB-timing */
)
{
UINT4 X, numDetectors;
MultiSSBtimes *ret = NULL;
INITSTATUS(status);
ATTATCHSTATUSPTR (status);
/* check input */
ASSERT (multiDetStates, status,COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
ASSERT (multiDetStates->length, status,COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
ASSERT (multiSSB, status,COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
ASSERT ( *multiBinary == NULL, status,COMPUTEFSTATC_ENONULL, COMPUTEFSTATC_MSGENONULL);
ASSERT (multiSSB != NULL, status,COMPUTEFSTATC_ENULL, COMPUTEFSTATC_MSGENULL);
numDetectors = multiDetStates->length;
if ( ( ret = LALCalloc( 1, sizeof( *ret ) )) == NULL ) {
ABORT (status, COMPUTEFSTATC_EMEM, COMPUTEFSTATC_MSGEMEM);
}
ret->length = numDetectors;
if ( ( ret->data = LALCalloc ( numDetectors, sizeof ( *ret->data ) )) == NULL ) {
LALFree ( ret );
ABORT (status, COMPUTEFSTATC_EMEM, COMPUTEFSTATC_MSGEMEM);
}
for ( X=0; X < numDetectors; X ++ )
{
SSBtimes *BinarytimesX = NULL;
UINT4 numStepsX = multiDetStates->data[X]->length;
ret->data[X] = LALCalloc ( 1, sizeof ( *(ret->data[X]) ) );
BinarytimesX = ret->data[X];
BinarytimesX->DeltaT = XLALCreateREAL8Vector ( numStepsX );
if ( (BinarytimesX->Tdot = XLALCreateREAL8Vector ( numStepsX )) == NULL ) {
XLALPrintError ("\nOut of memory!\n\n");
goto failed;
}
/* printf("calling LALGetBinarytimes\n"); */
LALGetBinarytimes (status->statusPtr, BinarytimesX, multiSSB->data[X], multiDetStates->data[X], binaryparams, refTime);
/* printf("finished LALGetBinarytimes\n"); */
if ( status->statusPtr->statusCode )
{
XLALPrintError ( "\nCall to LALGetBinarytimes() has failed ... \n\n");
goto failed;
}
// NOTE! It seems the old LAL-function *did not* set the reference time of the output binary-SSB vectors
// so we 'fix' this here retro-actively, in order to facilitate comparison with the new XLAL function:
BinarytimesX->refTime = multiSSB->data[X]->refTime;
} /* for X < numDet */
goto success;
failed:
/* free all memory allocated so far */
XLALDestroyMultiSSBtimes ( ret );
ABORT ( status, -1, "LALGetMultiBinarytimes failed" );
success:
(*multiBinary) = ret;
DETATCHSTATUSPTR (status);
RETURN(status);
} /* LALGetMultiBinarytimes() */
......@@ -6,6 +6,7 @@ XMLTEST = LALPulsarXMLTest
endif
check_PROGRAMS = \
BinarySSBTimesTest \
ExtrapolatePulsarSpinsTest \
FitToPulsarTest \
FlatMeshTest \
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment