32 #include <lal/LALStdlib.h>
33 #include <lal/FileIO.h>
34 #include <lal/DetectorSite.h>
35 #include <lal/LALError.h>
36 #include <lal/ConfigFile.h>
37 #include <lal/LogPrintf.h>
39 #include <lal/DopplerFullScan.h>
42 #define MIN(x,y) (x < y ? x : y)
43 #define MAX(x,y) (x > y ? x : y)
49 #define UNUSED __attribute__ ((unused))
61 typedef struct tagREAL8VectorList {
63 struct tagREAL8VectorList *
next;
64 struct tagREAL8VectorList *
prev;
105 DopplerFullScanState *
111 DopplerFullScanState *thisScan;
114 thisScan->gridType = init->
gridType;
122 switch ( thisScan->gridType ) {
129 "NOTE: gridType={metric,4,spin-square,spin-age-brk} only work for refTime (%f) == startTime (%f)!\n",
138 switch ( thisScan->gridType ) {
160 "\nGRID_SPINDOWN_{SQUARE,AGEBRK}: This option currently restricts the reference time to be the same as the start time.\n" );
164 "\nGRID_SPINDOWN_{SQUARE,AGEBRK}: gsl_vector_alloc failed\n" );
172 XLAL_CHECK_NULL( sky.numVertices == 1,
XLAL_EINVAL,
"\nGRID_SPINDOWN_{SQUARE,AGEBRK}: This option can only handle a single sky position.\n" );
178 if ( sky.vertices ) {
220 gsl_matrix_set_identity( metric );
229 gsl_matrix_free( metric );
253 DopplerFullScanState **scan,
283 if ( spinRange == NULL ) {
287 if ( scan == NULL ) {
291 ( *spinRange ) = scan->spinRange;
317 skyScanInit.gridType = init->
gridType;
322 skyScanInit.obsDuration = init->
Tspan;
324 skyScanInit.Detector = init->
Detector;
326 skyScanInit.skyGridFile = init->
gridFile;
331 scan->factoredScan = fscan;
361 nSpins[
i] = floor( scan->spinRange.fkdotBand[
i] / fscan->
skyScan.
dfkdot[
i] ) + 1.0;
367 scan->numTemplates = nTot;
368 XLALPrintInfo(
"Template grid: nSky x nFreq x nf1dot = %.0f x %.0f x %.0f = %.0f \n", nSky, nSpins[0], nSpins[1], nTot );
383 if ( ! scan->numTemplates ) {
384 switch ( scan->gridType ) {
400 LogPrintf(
LOG_NORMAL,
"template counting not implemented yet for gridType=%d!\n", scan->gridType );
406 return scan->numTemplates;
421 if (
pos == NULL || scan == NULL ) {
425 XLALPrintError(
"\nCalled XLALNextDopplerPos() on un-initialized DopplerFullScanState !\n\n" );
435 pos->refTime = scan->spinRange.refTime;
439 switch ( scan->gridType ) {
452 pos->fkdot[0] = scan->thisGridPoint->entry.data[0];
453 pos->Alpha = scan->thisGridPoint->entry.data[1];
454 pos->Delta = scan->thisGridPoint->entry.data[2];
455 pos->fkdot[1] = scan->thisGridPoint->entry.data[3];
456 pos->fkdot[2] = scan->thisGridPoint->entry.data[4];
457 pos->fkdot[3] = scan->thisGridPoint->entry.data[5];
460 if ( ( scan->thisGridPoint = scan->thisGridPoint->next ) == NULL ) {
497 XLALPrintError(
"\nGRID_SPINDOWN_{SQUARE,AGEBRK}: XLALNextLatticeTilingPoint() failed\n" );
504 pos->Alpha = gsl_vector_get( scan->spindownTilingPoint, 0 );
505 pos->Delta = gsl_vector_get( scan->spindownTilingPoint, 1 );
506 pos->fkdot[0] = gsl_vector_get( scan->spindownTilingPoint, 2 );
508 pos->fkdot[
i] = gsl_vector_get( scan->spindownTilingPoint,
i + 2 );
549 if (
pos == NULL || scan == NULL ) {
553 if ( ( fscan = scan->factoredScan ) == NULL ) {
557 range = &( scan->spinRange );
564 fkdotMax[3] =
range->fkdot[3] +
range->fkdotBand[3];
565 fkdotMax[2] =
range->fkdot[2] +
range->fkdotBand[2];
566 fkdotMax[1] =
range->fkdot[1] +
range->fkdotBand[1];
567 fkdotMax[0] =
range->fkdot[0] +
range->fkdotBand[0];
571 if ( nextPos.fkdot[0] > fkdotMax[0] ) {
572 nextPos.fkdot[0] =
range->fkdot[0];
574 if ( nextPos.fkdot[1] > fkdotMax[1] ) {
575 nextPos.fkdot[1] =
range->fkdot[1];
577 if ( nextPos.fkdot[2] > fkdotMax[2] ) {
578 nextPos.fkdot[2] =
range->fkdot[2];
580 if ( nextPos.fkdot[3] > fkdotMax[3] ) {
581 nextPos.fkdot[3] =
range->fkdot[3];
627 }
while ( (
ptr = next ) != NULL );
638 if ( scan == NULL ) {
642 if ( scan->factoredScan ) {
647 if ( scan->covering ) {
651 if ( scan->spindownTiling ) {
654 if ( scan->spindownTilingItr ) {
657 if ( scan->spindownTilingPoint ) {
658 gsl_vector_free( scan->spindownTilingPoint );
661 if ( scan->skyRegion.vertices ) {
662 XLALFree( scan->skyRegion.vertices );
679 if ( (
head == NULL ) || ( entry == NULL ) ) {
685 while (
ptr->next ) {
691 if ( ( newElement =
LALCalloc( 1,
sizeof( *newElement ) ) ) == NULL ) {
702 ptr->next = newElement;
755 XLAL_ERROR(
XLAL_EINVAL,
"\nnon-zero input spinRanges currently not supported! fkdot[%d] = %g, fkdotBand[%d] = %g\n\n",
779 while ( ! feof(
fp ) && ! ferror(
fp ) && fgets( line,
sizeof( line ),
fp ) != NULL ) {
782 if ( line[0] ==
'#' || line[0] ==
'%' ) {
787 REAL8 Freq, Alpha, Delta,
f1dot, f2dot, f3dot;
789 &Freq, &Alpha, &Delta, &
f1dot, &f2dot, &f3dot ) ) {
790 XLALPrintError(
"ERROR: Failed to parse 6 REAL8's from line %d in grid-file '%s'\n\n", numTemplates + 1, init->
gridFile );
798 alphaMin =
fmin( Alpha, alphaMin );
799 deltaMin =
fmin( Delta, deltaMin );
800 FreqMin =
fmin( Freq, FreqMin );
802 f2dotMin =
fmin( f2dot, f2dotMin );
803 f3dotMin =
fmin( f3dot, f3dotMin );
805 alphaMax =
fmax( Alpha, alphaMax );
806 deltaMax =
fmax( Delta, deltaMax );
807 FreqMax =
fmax( Freq, FreqMax );
809 f2dotMax =
fmax( f2dot, f2dotMax );
810 f3dotMax =
fmax( f3dot, f3dotMax );
813 entry->
data[0] = Freq;
814 entry->
data[1] = Alpha;
815 entry->
data[2] = Delta;
817 entry->
data[4] = f2dot;
818 entry->
data[5] = f3dot;
830 if ( ferror(
fp ) ) {
840 CHAR *skyRegionString = NULL;
849 scan->spinRange.fkdot[0] = FreqMin;
850 scan->spinRange.fkdotBand[0] = FreqMax - FreqMin;
852 scan->spinRange.fkdot[1] = f1dotMin;
853 scan->spinRange.fkdotBand[1] = f1dotMax - f1dotMin;
855 scan->spinRange.fkdot[2] = f2dotMin;
856 scan->spinRange.fkdotBand[2] = f2dotMax - f2dotMin;
858 scan->spinRange.fkdot[3] = f3dotMin;
859 scan->spinRange.fkdotBand[3] = f3dotMax - f3dotMin;
861 scan->numTemplates = numTemplates;
862 scan->covering =
head.next;
863 scan->thisGridPoint = scan->covering;
865 XLALPrintInfo(
"Template grid: nTot = %.0f\n", 1.0 * numTemplates );
866 XLALPrintInfo(
"Spanned ranges: Freq in [%g, %g], f1dot in [%g, %g], f2dot in [%g, %g], f3dot in [%g, %g]\n",
867 FreqMin, FreqMax, f1dotMin, f1dotMax, f2dotMin, f2dotMax, f3dotMin, f3dotMax );
880 const size_t dim UNUSED,
881 const gsl_matrix *
cache UNUSED,
882 const gsl_vector *point
890 const double freq = gsl_vector_get( point,
info->freq_dim );
898 LatticeTiling *tiling,
899 const size_t freq_dimension,
900 const size_t f1dot_dimension,
902 const double min_braking,
903 const double max_braking
918 info_lower.freq_dim = info_upper.freq_dim = freq_dimension;
919 info_lower.scale = -1.0 / ( ( min_braking - 1.0 ) * age );
920 info_upper.scale = -1.0 / ( ( max_braking - 1.0 ) * age );
935 const size_t dim UNUSED,
936 const gsl_matrix *
cache UNUSED,
937 const gsl_vector *point
945 const double freq = gsl_vector_get( point,
info->freq_dim );
946 const double f1dot = gsl_vector_get( point,
info->f1dot_dim );
954 LatticeTiling *tiling,
955 const size_t freq_dimension,
956 const size_t f1dot_dimension,
957 const size_t f2dot_dimension,
958 const double min_braking,
959 const double max_braking
974 info_lower.freq_dim = info_upper.freq_dim = freq_dimension;
975 info_lower.f1dot_dim = info_upper.f1dot_dim = f1dot_dimension;
976 info_lower.scale = min_braking;
977 info_upper.scale = max_braking;
static double F2DotBrakingBound(const void *data, const size_t dim UNUSED, const gsl_matrix *cache UNUSED, const gsl_vector *point)
int nextPointInFactoredGrid(PulsarDopplerParams *pos, DopplerFullScanState *scan)
return current grid-point and step forward one template in 'factored' grids (sky x f0dot x f1dot ....
int XLALInitFactoredGrid(DopplerFullScanState *scan, const DopplerFullScanInit *init)
Initialize Doppler-scanner to emulate an old-style factored template grid: 'sky x f0dot x f1dot x f2d...
void XLALDestroyDopplerFullScan(DopplerFullScanState *scan)
Destroy the a full DopplerFullScanState structure.
REAL8 XLALNumDopplerTemplates(DopplerFullScanState *scan)
Return (and compute if not done so yet) the total number of Doppler templates of the DopplerScan scan...
int XLALSetLatticeTilingF2DotBrakingBound(LatticeTiling *tiling, const size_t freq_dimension, const size_t f1dot_dimension, const size_t f2dot_dimension, const double min_braking, const double max_braking)
Set a second spindown bound derived from braking indices.
static void XLALREAL8VectorListDestroy(REAL8VectorList *head)
static REAL8VectorList * XLALREAL8VectorListAddEntry(REAL8VectorList *head, const REAL8Vector *entry)
void InitDopplerFullScan(LALStatus *status, DopplerFullScanState **scan, const DopplerFullScanInit *init)
int XLALGetDopplerSpinRange(PulsarSpinRange *spinRange, const DopplerFullScanState *scan)
Return the spin-range spanned by the given template bank stored in the opaque DopplerFullScanState.
int XLALLoadFullGridFile(DopplerFullScanState *scan, const DopplerFullScanInit *init)
load a full multi-dim template grid from the file init->gridFile, the file-format is: lines of 6 colu...
int XLALNextDopplerPos(PulsarDopplerParams *pos, DopplerFullScanState *scan)
Function to step through the full template grid point by point.
DopplerFullScanState * XLALInitDopplerFullScan(const DopplerFullScanInit *init)
Set up a full multi-dimensional grid-scan.
int XLALSetLatticeTilingF1DotAgeBrakingBound(LatticeTiling *tiling, const size_t freq_dimension, const size_t f1dot_dimension, const double age, const double min_braking, const double max_braking)
Set a first spindown bound derived from spindown age and braking indices.
static double F1DotAgeBrakingBound(const void *data, const size_t dim UNUSED, const gsl_matrix *cache UNUSED, const gsl_vector *point)
CHAR * XLALSkySquare2String(REAL8 Alpha, REAL8 Delta, REAL8 AlphaBand, REAL8 DeltaBand)
parse a 'classical' sky-square (Alpha, Delta, AlphaBand, DeltaBand) into a "SkyRegion"-string of the ...
int XLALParseSkyRegionString(SkyRegion *region, const CHAR *input)
parse a skyRegion-string into a SkyRegion structure: the expected string-format is " (longitude1,...
int XLALInitDopplerSkyScan(DopplerSkyScanState *skyScan, const DopplerSkyScanInit *init)
Initialize the Doppler sky-scanner.
void XLALDestroyDopplerSkyScan(DopplerSkyScanState *skyScan)
Destroy the DopplerSkyScanState structure.
#define DOPPLERSCANH_MSGENONULL
#define DOPPLERSCANH_ENONULL
#define DOPPLERSCANH_MSGEXLAL
#define DOPPLERSCANH_EINPUT
#define DOPPLERSCANH_MSGEINPUT
#define DOPPLERSCANH_EXLAL
@ STATE_READY
initialized and ready
@ STATE_FINISHED
all templates have been read
@ STATE_IDLE
not initialized yet
DopplerGridType
different types of grids:
@ GRID_METRIC
generate grid using a 2D sky-metric
@ GRID_SPINDOWN_SQUARE
spindown tiling for a single sky position and square parameter space
@ GRID_FLAT
"flat" sky-grid: fixed step-size (dAlpha,dDelta)
@ GRID_FILE_SKYGRID
read skygrid from a file
@ GRID_FILE_FULLGRID
load the full D-dim grid from a file
@ GRID_METRIC_SKYFILE
'hybrid' grid-construction: use skygrid from file, metric for others
@ GRID_ISOTROPIC
approximately isotropic sky-grid
@ GRID_SPINDOWN_AGEBRK
spindown tiling for a single sky position and non-square parameter space defined by spindown age and ...
#define ABORT(statusptr, code, mesg)
#define ASSERT(assertion, statusptr, code, mesg)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define XLAL_INIT_DECL(var,...)
int XLALSetLatticeTilingBound(LatticeTiling *tiling, const size_t dim, const LatticeTilingBound func, const size_t data_len, const void *data_lower, const void *data_upper)
Set a parameter-space bound on a dimension of the lattice tiling.
UINT8 XLALTotalLatticeTilingPoints(const LatticeTilingIterator *itr)
Return the total number of points covered by the lattice tiling iterator.
int XLALSetLatticeTilingConstantBound(LatticeTiling *tiling, const size_t dim, const double bound1, const double bound2)
Set a constant lattice tiling parameter-space bound, given by the minimum and maximum of the two supp...
int XLALSetTilingLatticeAndMetric(LatticeTiling *tiling, const TilingLattice lattice, const gsl_matrix *metric, const double max_mismatch)
Set the tiling lattice, parameter-space metric, and maximum prescribed mismatch.
LatticeTilingIterator * XLALCreateLatticeTilingIterator(const LatticeTiling *tiling, const size_t itr_ndim)
Create a new lattice tiling iterator.
void XLALDestroyLatticeTilingIterator(LatticeTilingIterator *itr)
Destroy a lattice tiling iterator.
void XLALDestroyLatticeTiling(LatticeTiling *tiling)
Destroy a lattice tiling.
int XLALNextLatticeTilingPoint(LatticeTilingIterator *itr, gsl_vector *point)
Advance lattice tiling iterator, and optionally return the next point in point.
LatticeTiling * XLALCreateLatticeTiling(const size_t ndim)
Create a new lattice tiling.
int XLALSetLatticeTilingPadding(LatticeTiling *tiling, const size_t dim, const double lower_bbox_pad, const double upper_bbox_pad, const int lower_intp_pad, const int upper_intp_pad, const int find_bound_extrema)
Control the padding of lattice tiling parameter-space bounds in the given dimension.
@ TILING_LATTICE_ANSTAR
An-star ( ) lattice.
void void LogPrintfVerbatim(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
void LogPrintf(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
int XLALSpindownMetric(gsl_matrix *metric, double Tspan)
Frequency and frequency derivative components of the metric, suitable for a directed search with only...
REAL8 PulsarSpins[PULSAR_MAX_SPINS]
Typedef for fixed-size array holding GW frequency and derivatives fk = d^k Freq/dt^k|(tau_ref)
#define PULSAR_MAX_SPINS
maximal number of spin-parameters (Freq + spindowns) we can handle
void XLALNormalizeSkyPosition(double *restrict longitude, double *restrict latitude)
COORDINATESYSTEM_EQUATORIAL
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
#define XLAL_ERROR_NULL(...)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK_NULL(assertion,...)
#define XLAL_PRINT_DEPRECATION_WARNING(replacement)
int XLALGPSCmp(const LIGOTimeGPS *t0, const LIGOTimeGPS *t1)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
REAL8 XLALGPSDiff(const LIGOTimeGPS *t1, const LIGOTimeGPS *t0)
Structure describing a region in paramter-space (a,d,f,f1dot,..).
LALPulsarMetricType metricType
which metric to use if GRID_METRIC
const CHAR * gridFile
filename for sky-grid or full-grid if GRID_FILE_SKYGRID or GRID_FILE_FULLGRID
REAL8 metricMismatch
for GRID_METRIC and GRID_ISOTROPIC
REAL8 extraArgs[3]
extra grid-specific setup arguments
const LALDetector * Detector
Current detector.
LIGOTimeGPS startTime
start-time of the observation
REAL8 Tspan
total time spanned by observation
DopplerRegion searchRegion
Doppler-space region to be covered + scanned.
BOOLEAN projectMetric
project metric on f=const subspace
PulsarDopplerParams stepSizes
user-settings for stepsizes if GRID_FLAT
const EphemerisData * ephemeris
ephemeris-data for numerical metrics
DopplerGridType gridType
which type of grid to generate
PulsarSpins fkdotBand
spin-intervals
CHAR * skyRegionString
sky-region string '(a1,d1), (a2,d2), ..'
PulsarSpins fkdot
reference time of definition of Doppler parameters
struct tagDopplerSkyGrid * next
initialization-structure passed to InitDopplerSkyScan()
this structure reflects the current state of a DopplerSkyScan
UINT4 numSkyGridPoints
how many skygrid-points
PulsarSpins dfkdot
fixed-size steps in spins
DopplerSkyGrid * skyNode
pointer to current grid-node in skygrid
scan_state_t state
idle, ready or finished
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
PulsarSpins fkdot
Intrinsic spins: [Freq, f1dot, f2dot, ...
REAL8 Delta
Sky position: DEC (latitude) in equatorial coords and radians.
LIGOTimeGPS refTime
Reference time of pulsar parameters (in SSB!)
REAL8 Alpha
Sky position: RA (longitude) in equatorial coords and radians.
Contains a "spin-range", ie spins and corresponding bands at a given (SSB) reference GPS-time .
doubly linked list of REAL8-vectors (physical vectors)
struct tagREAL8VectorList * next
struct tagREAL8VectorList * prev
structure describing a polygon-region in the sky
PulsarDopplerParams thisPoint
current doppler-position of the scan
DopplerSkyScanState skyScan
keep track of sky-grid stepping
--— internal [opaque] type to store the state of a FULL multidimensional grid-scan --—
REAL8VectorList * covering
multi-dimensional covering
SkyRegion skyRegion
sky-range covered by template bank
factoredGridScan_t * factoredScan
only used to emulate FACTORED grids sky x Freq x f1dot
REAL8VectorList * thisGridPoint
pointer to current grid-point
INT2 state
idle, ready or finished
LatticeTiling * spindownTiling
spindown lattice tiling
DopplerGridType gridType
what type of grid are we dealing with
PulsarSpinRange spinRange
spin-range covered by template bank
REAL8 numTemplates
total number of templates in the grid
LatticeTilingIterator * spindownTilingItr
iterator over spindown lattice tiling
gsl_vector * spindownTilingPoint
current point in spindown lattice tiling