Commit 6f7b430f authored by Reinhard Prix's avatar Reinhard Prix

extended XLALMakeTimestamps() to overlapping SFTs, added XLALMakeMultiTimestamps()

   - allow XLALMakeTimestamps() to explicitly generate overlapping SFT
     timestamps, and make sure end-time is 'covered' only once
     (replacing some ugly code from mfdv4 that previously dealt with that)
   - added trivial multi-IFO generalization XLALMakeMultiTimestamps()
   - refs #763
Original: 1b0133e8beee228b3217b6eafd4c6546dd2e842c
parent 8fc417cf
......@@ -972,9 +972,6 @@ XLALInitMakefakedata ( ConfigVars_t *cfg, UserVariables_t *uvar )
/* have we got our timestamps yet?: If not, we must get them from (start, duration) user-input */
if ( ! cfg->timestamps )
{
REAL8 tStep = uvar->Tsft - uvar->SFToverlap;
LIGOTimeGPS tStart;
REAL8 t0, tLast;
if ( !haveStart || !haveDuration ) {
XLAL_ERROR ( XLAL_EINVAL, "Need to have either --timestampsFile OR (--startTime,--duration) OR --noiseSFTs\n\n");
}
......@@ -982,22 +979,10 @@ XLALInitMakefakedata ( ConfigVars_t *cfg, UserVariables_t *uvar )
XLAL_ERROR ( XLAL_EINVAL, "--SFToverlap cannot be larger than --Tsft!\n\n");
}
LIGOTimeGPS tStart;
/* internally always use timestamps, so we generate them */
XLALGPSSetREAL8 ( &tStart, uvar->startTime );
cfg->timestamps = XLALMakeTimestamps ( tStart, uvar->duration, tStep );
XLAL_CHECK ( cfg->timestamps != NULL, XLAL_EFUNC );
/* "prune" last timestamp(s) if the one before-last also covers the end
* (this happens for overlapping SFTs with LALMakeTimestamps() as used above
*/
t0 = XLALGPSGetREAL8( &(cfg->timestamps->data[0]) );
tLast = XLALGPSGetREAL8 ( &(cfg->timestamps->data[cfg->timestamps->length - 1 ]) );
while ( tLast - t0 + uvar->Tsft > uvar->duration + 1e-6)
{
cfg->timestamps->length --;
tLast = XLALGPSGetREAL8 ( &(cfg->timestamps->data[cfg->timestamps->length - 1 ]) );
}
XLAL_CHECK ( ( cfg->timestamps = XLALMakeTimestamps ( tStart, uvar->duration, uvar->Tsft, uvar->SFToverlap )) != NULL, XLAL_EFUNC );
} /* if !cfg->timestamps */
......
......@@ -636,33 +636,17 @@ XLALInitMakefakedata ( ConfigVars_t *cfg, UserVariables_t *uvar )
/* have we got our timestamps yet?: If not, we must get them from (start, duration) user-input */
if ( ! cfg->timestamps )
{
REAL8 tStep = uvar->Tsft - uvar->SFToverlap;
LIGOTimeGPS tStart;
REAL8 t0, tLast;
if ( !haveStart || !haveDuration ) {
XLAL_ERROR ( XLAL_EINVAL, "Need to have either --timestampsFile OR (--startTime,--duration) OR --noiseSFTs\n\n");
}
if ( uvar->SFToverlap > uvar->Tsft ) {
XLAL_ERROR ( XLAL_EINVAL, "--SFToverlap cannot be larger than --Tsft!\n\n");
}
/* internally always use timestamps, so we generate them */
XLALGPSSetREAL8 ( &tStart, uvar->startTime );
cfg->timestamps = XLALMakeTimestamps ( tStart, uvar->duration, tStep );
cfg->timestamps = XLALMakeTimestamps ( tStart, uvar->duration, uvar->Tsft, uvar->SFToverlap );
XLAL_CHECK ( cfg->timestamps != NULL, XLAL_EFUNC );
/* "prune" last timestamp(s) if the one before-last also covers the end
* (this happens for overlapping SFTs with LALMakeTimestamps() as used above
*/
t0 = XLALGPSGetREAL8( &(cfg->timestamps->data[0]) );
tLast = XLALGPSGetREAL8 ( &(cfg->timestamps->data[cfg->timestamps->length - 1 ]) );
while ( tLast - t0 + uvar->Tsft > uvar->duration + 1e-6)
{
cfg->timestamps->length --;
tLast = XLALGPSGetREAL8 ( &(cfg->timestamps->data[cfg->timestamps->length - 1 ]) );
}
} /* if !cfg->timestamps */
} /* END: setup signal start + duration */
......@@ -759,7 +743,7 @@ XLALInitMakefakedata ( ConfigVars_t *cfg, UserVariables_t *uvar )
{
/* convert MJD peripase to GPS using Matt Pitkins code found at lal/packages/pulsar/src/BinaryPulsarTimeing.c */
REAL8 GPSfloat;
GPSfloat = LALTTMJDtoGPS(uvar->orbitTpSSBMJD);
GPSfloat = XLALTTMJDtoGPS(uvar->orbitTpSSBMJD);
XLALGPSSetREAL8(&(orbit->tp),GPSfloat);
}
else if ((set5 && set6) && !set7)
......@@ -802,7 +786,7 @@ XLALInitMakefakedata ( ConfigVars_t *cfg, UserVariables_t *uvar )
/* convert MJD to GPS using Matt Pitkins code found at lal/packages/pulsar/src/BinaryPulsarTimeing.c */
REAL8 GPSfloat;
GPSfloat = LALTTMJDtoGPS(uvar->refTimeMJD);
GPSfloat = XLALTTMJDtoGPS(uvar->refTimeMJD);
XLALGPSSetREAL8(&(cfg->pulsar.Doppler.refTime),GPSfloat);
}
else
......
......@@ -299,8 +299,8 @@ XLALSignalToSFTs ( const REAL4TimeSeries *signalvec, /**< input time-series */
LIGOTimeGPSVector *timestamps;
if ( params->timestamps == NULL )
{
timestamps = XLALMakeTimestamps ( tStart, duration, params->Tsft );
XLAL_CHECK_NULL ( timestamps != NULL, XLAL_EFUNC );
REAL8 Toverlap = 0;
XLAL_CHECK_NULL ( ( timestamps = XLALMakeTimestamps ( tStart, duration, params->Tsft, Toverlap )) != NULL, XLAL_EFUNC );
/* see if the last timestamp is valid (can fit a full SFT in?), if not, drop it */
LIGOTimeGPS lastTs = timestamps->data[timestamps->length-1];
if ( XLALGPSDiff ( &lastTs, &tLast ) > 0 ) // if lastTs > tLast
......
......@@ -535,7 +535,7 @@ LALMakeTimestamps ( LALStatus *status, /**< pointer to LALStatus structure */
ASSERT (*timestamps == NULL,status, SFTUTILS_ENONULL, SFTUTILS_MSGENONULL);
LIGOTimeGPSVector *ts;
ts = XLALMakeTimestamps ( tStart, duration, tStep );
ts = XLALMakeTimestamps ( tStart, duration, tStep, 0 );
if ( ts == NULL )
{
XLALPrintError ("XLALMakeTimestamps() failed with xlalErrno = %d\n", xlalErrno );
......
......@@ -43,6 +43,9 @@
/*---------- internal types ----------*/
/*---------- Global variables ----------*/
static REAL8 fudge_up = 1 + 10 * LAL_REAL8_EPS; // about ~1 + 2e-15
static REAL8 fudge_down = 1 - 10 * LAL_REAL8_EPS; // about ~1 - 2e-15
/* empty struct initializers */
const PSDVector empty_PSDVector;
const MultiPSDVector empty_MultiPSDVector;
......@@ -277,47 +280,86 @@ XLALDestroyTimestampVector ( LIGOTimeGPSVector *vect)
} /* XLALDestroyTimestampVector() */
/** Given a start-time, duration and 'stepsize' tStep, returns a list of timestamps
* covering this time-stretch.
/** Given a start-time, Tspan, Tsft and Toverlap, returns a list of timestamps
* covering this time-stretch (allowing for overlapping SFT timestamps).
*
* NOTE: boundary-handling: the returned list of timestamps are guaranteed to *cover* the
* interval [tStart, tStart+duration], assuming a each timestamp covers a length of 'tStep'
* This implies that the actual timestamps-coverage can extend up to 'tStep' beyond 'tStart+duration'.
* interval [tStart, tStart+duration], assuming a each timestamp covers a length of 'Tsft'
* This implies that the actual timestamps-coverage can extend up to 'Tsft' beyond 'tStart+duration'.
*/
LIGOTimeGPSVector *
XLALMakeTimestamps ( LIGOTimeGPS tStart, /**< GPS start-time */
REAL8 duration, /**< duration in seconds */
REAL8 tStep /**< length of one (SFT) timestretch in seconds */
XLALMakeTimestamps ( LIGOTimeGPS tStart, /**< GPS start-time */
REAL8 Tspan, /**< total duration to cover, in seconds */
REAL8 Tsft, /**< Tsft: SFT length of each timestamp, in seconds */
REAL8 Toverlap /**< time to overlap successive SFTs by, in seconds */
)
{
XLAL_CHECK_NULL ( tStep > 0, XLAL_EDOM, "Invalid non-positive input 'tStart = %g'\n", tStep );
XLAL_CHECK_NULL ( duration > 0, XLAL_EDOM, "Invalid non-positive input 'duration = %g'\n", duration );
UINT4 numSFTs = ceil( duration / tStep ); /* >= 1 !*/
XLAL_CHECK_NULL ( Tspan > 0, XLAL_EDOM );
XLAL_CHECK_NULL ( Tsft > 0, XLAL_EDOM );
XLAL_CHECK_NULL ( Toverlap >= 0, XLAL_EDOM );
XLAL_CHECK_NULL ( Toverlap < Tsft, XLAL_EDOM ); // we must actually advance
REAL8 Tstep = Tsft - Toverlap; // guaranteed > 0
UINT4 numSFTsMax = ceil ( Tspan * fudge_down / Tstep ); /* >= 1 !*/
// now we might be covering the end-time several times, if using overlapping SFTs, so
// let's trim this back down so that end-time is covered exactly once
UINT4 numSFTs = numSFTsMax;
while ( (numSFTs >= 2) && ( (numSFTs - 2) * Tstep + Tsft > Tspan) ) {
numSFTs --;
}
LIGOTimeGPSVector *ts;
XLAL_CHECK_NULL ( (ts = XLALCreateTimestampVector ( numSFTs )) != NULL, XLAL_EFUNC );
LIGOTimeGPSVector *ret;
XLAL_CHECK_NULL ( (ret = XLALCreateTimestampVector ( numSFTs )) != NULL, XLAL_EFUNC );
ts->deltaT = tStep;
ret->deltaT = Tsft;
LIGOTimeGPS tt = tStart; /* initialize to start-time */
for (UINT4 i = 0; i < numSFTs; i++)
for ( UINT4 i = 0; i < numSFTs; i++ )
{
ts->data[i] = tt;
ret->data[i] = tt;
/* get next time-stamp */
/* NOTE: we add the interval tStep successively (rounded correctly to ns each time!)
* instead of using iSFT*Tsft, in order to avoid possible ns-rounding problems
* with REAL8 intervals, which becomes critial from about 100days on...
*/
XLAL_CHECK_NULL ( XLALGPSAdd ( &tt, tStep ) != NULL, XLAL_EFUNC );
XLAL_CHECK_NULL ( XLALGPSAdd ( &tt, Tstep ) != NULL, XLAL_EFUNC );
} /* for i < numSFTs */
return ts;
return ret;
} /* XLALMakeTimestamps() */
/**
* Same as XLALMakeTimestamps() just for several detectors,
* additionally specify the number of detectors.
*/
MultiLIGOTimeGPSVector *
XLALMakeMultiTimestamps ( LIGOTimeGPS tStart, /**< GPS start-time */
REAL8 Tspan, /**< total duration to cover, in seconds */
REAL8 Tsft, /**< Tsft: SFT length of each timestamp, in seconds */
REAL8 Toverlap, /**< time to overlap successive SFTs by, in seconds */
UINT4 numDet /**< number of timestamps-vectors to generate */
)
{
XLAL_CHECK_NULL ( numDet >= 1, XLAL_EINVAL );
MultiLIGOTimeGPSVector *ret;
XLAL_CHECK_NULL ( ( ret = XLALCalloc ( 1, sizeof(*ret))) != NULL, XLAL_ENOMEM );
XLAL_CHECK_NULL ( ( ret->data = XLALCalloc ( numDet, sizeof(ret->data[0]) )) != NULL, XLAL_ENOMEM );
ret->length = numDet;
for ( UINT4 X=0; X < numDet; X ++ )
{
XLAL_CHECK_NULL ( (ret->data[X] = XLALMakeTimestamps ( tStart, Tspan, Tsft, Toverlap ) ) != NULL, XLAL_EFUNC );
} // for X < numDet
return ret;
} /* XLALMakeMultiTimestamps() */
/** Extract timstamps-vector from the given SFTVector
*/
LIGOTimeGPSVector *
......@@ -868,9 +910,6 @@ XLALExtractBandFromSFTVector ( const SFTVector *inSFTs, REAL8 fmin, REAL8 Band )
XLAL_CHECK_NULL ( (fmin >= SFTf0) && ( fmin + Band <= SFTf0 + SFTBand ), XLAL_EINVAL,
"Requested frequency-band [%f,%f] Hz not contained SFTs [%f, %f] Hz.\n", fmin, fmin + Band, SFTf0, SFTf0 + SFTBand );
REAL8 eps = 10 * LAL_REAL8_EPS; // about ~2e-15
REAL8 fudge_up = 1 + eps;
REAL8 fudge_down = 1 - eps;
UINT4 firstBin = floor ( fmin*fudge_up / df ); // round *down*, allowing for eps 'fudge'
REAL8 fmax = fmin + Band;
UINT4 lastBin = ceil ( ( fmax*fudge_down) / df ); // round *up*, allowing for eps fudge
......
......@@ -113,7 +113,8 @@ COMPLEX8Vector *XLALrefineCOMPLEX8Vector (const COMPLEX8Vector *in, UINT4 refine
SFTVector *XLALExtractBandFromSFTVector ( const SFTVector *inSFTs, REAL8 fmin, REAL8 Band );
LIGOTimeGPSVector *XLALCreateTimestampVector (UINT4 len);
LIGOTimeGPSVector *XLALMakeTimestamps ( LIGOTimeGPS tStart, REAL8 duration, REAL8 tStep );
LIGOTimeGPSVector *XLALMakeTimestamps ( LIGOTimeGPS tStart, REAL8 Tspan, REAL8 Tsft, REAL8 Toverlap );
MultiLIGOTimeGPSVector *XLALMakeMultiTimestamps ( LIGOTimeGPS tStart, REAL8 Tspan, REAL8 Tsft, REAL8 Toverlap, UINT4 numDet );
LIGOTimeGPSVector *XLALExtractTimestampsFromSFTs ( const SFTVector *sfts );
MultiLIGOTimeGPSVector *XLALExtractMultiTimestampsFromSFTs ( const MultiSFTVector *multiSFTs );
......
......@@ -140,8 +140,8 @@ XLALSimulateExactPulsarSignal ( const PulsarSignalParams *params )
/* get timestamps of timeseries plus detector-states */
REAL8 dt = 1.0 / params->samplingRate;
LIGOTimeGPSVector *timestamps = XLALMakeTimestamps ( params->startTimeGPS, params->duration, dt );
XLAL_CHECK_NULL ( timestamps != NULL, XLAL_EFUNC, "XLALMakeTimestamps() failed.\n");
LIGOTimeGPSVector *timestamps;
XLAL_CHECK_NULL ( (timestamps = XLALMakeTimestamps ( params->startTimeGPS, params->duration, dt, 0 )) != NULL, XLAL_EFUNC );
UINT4 numSteps = timestamps->length;
......
......@@ -149,7 +149,7 @@ main ( int argc, char *argv[] )
for ( UINT4 X = 0; X < numDetectors; X ++ )
{
multiTS->data[X] = XLALMakeTimestamps ( startTime, duration, Tsft );
multiTS->data[X] = XLALMakeTimestamps ( startTime, duration, Tsft, 0 );
XLAL_CHECK ( multiTS->data[X] != NULL, XLAL_EFUNC, "XLALMakeTimestamps() failed.\n");
} /* for X < numIFOs */
......
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