BinarySSBTimesTest.c 20.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
/*
 * 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>

24
#include <lal/LALPulsarVCSInfo.h>
25
#include <lal/ComputeFstat.h>
26
#include <lal/LALgetopt.h>
27 28
#include <lal/LALInitBarycenter.h>
#include <lal/FindRoot.h>
29
#include <lal/UserInput.h>
30

31 32
/**
 * \author Reinhard Prix
33
 * \file
34
 * \ingroup lalpulsar_coh
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
 * \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
52
static REAL8 A,B;          /* binary time delay coefficients (need to be global so that the LAL root finding procedure can see them) */
53

54 55 56 57 58
// local types
typedef struct
{
  INT4 randSeed;	/**< allow user to specify random-number seed for reproducible noise-realizations */
} UserInput_t;
59

60 61 62 63 64 65 66 67 68
/* Type defining the orbital parameters of a binary pulsar */
typedef struct tagBinaryOrbitParams {
  LIGOTimeGPS tp;         /* time of observed periapsis passage (in SSB) */
  REAL8 argp;             /* argument of periapsis (radians) */
  REAL8 asini;            /* projected, normalized orbital semi-major axis (s) */
  REAL8 ecc;              /* orbital eccentricity */
  REAL8 period;           /* orbital period (sec) */
} BinaryOrbitParams;

69 70 71
/* ----- 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);
72
static void EccentricAnomoly(LALStatus *status, REAL8 *tr, REAL8 lE, void *x0);
73 74 75 76 77 78 79 80 81

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[] )
{
82 83
  LALStatus XLAL_INIT_DECL(status);
  UserInput_t XLAL_INIT_DECL(uvar_s);
84 85
  UserInput_t *uvar = &uvar_s;

86 87 88
  struct tms buf;
  uvar->randSeed = times(&buf);

89
  // ---------- register all our user-variable ----------
90
  XLALRegisterUvarMember(   randSeed,             INT4, 's', OPTIONAL, "Specify random-number seed for reproducible noise.");
91 92

  /* read cmdline & cfgfile  */
93
  BOOLEAN should_exit = 0;
94
  XLAL_CHECK( XLALUserVarReadAllInput( &should_exit, argc, argv, lalPulsarVCSInfoList ) == XLAL_SUCCESS, XLAL_EFUNC );
95 96
  if ( should_exit ) {
    exit (1);
97 98
  }

99 100
  srand ( uvar->randSeed );

101 102 103
  REAL8 startTimeREAL8 	= 714180733;
  REAL8 duration 	= 180000;	/* 50 hours */
  REAL8 Tsft 		= 1800;		/* assume 30min SFTs */
104 105
  char earthEphem[] 	= TEST_PKG_DATA_DIR "earth00-19-DE200.dat.gz";
  char sunEphem[]   	= TEST_PKG_DATA_DIR "sun00-19-DE200.dat.gz";
106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134

  //REAL8 tolerance = 2e-10;	/* same algorithm, should be basically identical results */

  LIGOTimeGPS startTime, refTime;
  XLALGPSSetREAL8 ( &startTime, startTimeREAL8 );
  refTime = startTime;

  // 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] );

135 136
  MultiLALDetector multiIFO;
  multiIFO.length = numDetectors;
137 138 139 140
  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 );
141
      multiIFO.sites[X] = (*det);	 // struct copy
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
      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 ++ )
    {
158
      multiTS->data[X] = XLALMakeTimestamps ( startTime, duration, Tsft, 0 );
159 160 161 162
      XLAL_CHECK ( multiTS->data[X] != NULL, XLAL_EFUNC, "XLALMakeTimestamps() failed.\n");
    } /* for X < numIFOs */

  // generate detector-states
163
  MultiDetectorStateSeries *multiDetStates = XLALGetMultiDetectorStates ( multiTS, &multiIFO, edat, 0 );
164 165 166 167 168 169 170 171 172 173 174 175 176
  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;
177 178
  PulsarDopplerParams doppler;
  memset(&doppler, 0, sizeof(doppler));
179
  doppler.fkdot[0] = 100;
180 181 182 183 184 185
  doppler.tp = orbit.tp;
  doppler.argp = orbit.argp;
  doppler.asini = orbit.asini;
  doppler.ecc = orbit.ecc;
  doppler.period = orbit.period;
  XLAL_CHECK ( XLALAddMultiBinaryTimes ( &multiBinary_test, multiSSBIn, &doppler ) == XLAL_SUCCESS, XLAL_EFUNC );
186 187 188

  // ----- step 3: compare results
  REAL8 err_DeltaT, err_Tdot;
189
  REAL8 tolerance = 1e-9;
190 191 192 193 194 195 196 197 198
  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
199
  XLALDestroyUserVars();
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290
  XLALDestroyEphemerisData ( edat );
  XLALDestroyMultiSSBtimes ( multiBinary_test );
  XLALDestroyMultiSSBtimes ( multiBinary_ref );
  XLALDestroyMultiSSBtimes ( multiSSBIn );
  XLALDestroyMultiTimestamps ( multiTS );
  XLALDestroyMultiDetectorStateSeries ( multiDetStates );

  // 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

291 292 293 294 295 296 297 298 299 300 301 302 303 304
#define COMPUTEFSTATC_ENULL 		1
#define COMPUTEFSTATC_ENONULL 		2
#define COMPUTEFSTATC_EINPUT   		3
#define COMPUTEFSTATC_EMEM   		4
#define COMPUTEFSTATC_EXLAL		5
#define COMPUTEFSTATC_EIEEE		6

#define COMPUTEFSTATC_MSGENULL 		"Arguments contained an unexpected null pointer"
#define COMPUTEFSTATC_MSGENONULL 	"Output pointer is non-NULL"
#define COMPUTEFSTATC_MSGEINPUT   	"Invalid input"
#define COMPUTEFSTATC_MSGEMEM   	"Out of memory. Bad."
#define COMPUTEFSTATC_MSGEXLAL		"XLAL function call failed"
#define COMPUTEFSTATC_MSGEIEEE		"Floating point failure"

305 306 307 308
/**
 * 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$.
309
 *
310 311
 * \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.
312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362
 *
 */
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 */
363 364 365 366
  A  = (LAL_TWOPI/Porb)*cosw*asini*sqrt(1.0-e*e) - e;
  B  = (LAL_TWOPI/Porb)*sinw*asini;
  REAL8 tp = GPS2REAL8(binaryparams->tp);
  REAL8 Tp = tp  + asini*sinw*ome;
367 368 369 370 371 372 373 374 375 376 377 378 379

  /* 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 */
380 381 382
      REAL8 temp = fmod((tSSB_now - Tp),Porb)/(REAL8)Porb;
      fracorb = temp - (REAL8)floor(temp);
      REAL8 x0 = fracorb * LAL_TWOPI;
383 384 385 386 387 388 389 390

      /* 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 */
391
      LALDBracketRoot(status->statusPtr,&input,&x0);
392 393

      /* bisect domain to find eccentric anomoly E corresponding to the SSB time of the midpoint of this SFT */
394
      LALDBisectionFindRoot(status->statusPtr,&E,&input,&x0);
395 396 397 398 399

      /* 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 */
400
      tBinary->Tdot->data[i] = tSSB->Tdot->data[i] * ( (1.0 - e*cos(E))/(1.0 + A*cos(E) - B*sin(E)) );
401 402 403 404 405 406 407 408 409

    } /* for i < numSteps */


  DETATCHSTATUSPTR (status);
  RETURN(status);

} /* LALGetBinarytimes() */

410 411 412
/**
 * For a given set of binary parameters we solve the following function for
 * the eccentric anomoly E
413 414 415 416
 */
static void EccentricAnomoly(LALStatus *status,
			     REAL8 *tr,
			     REAL8 lE,
417
			     void *x0
418 419 420
			     )
{
  INITSTATUS(status);
421
  ASSERT(x0,status, 1, "Null pointer");
422 423

  /* this is the function relating the observed time since periapse in the SSB to the true eccentric anomoly E */
424
  (*tr) = - (*(REAL8 *)x0) + (lE + A*sin(lE) + B*(cos(lE) - 1.0));
425 426 427 428

  RETURN(status);
}

429 430
/**
 * Multi-IFO version of LALGetBinarytimes().
431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507
 * 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() */