Loading [MathJax]/jax/output/HTML-CSS/config.js
LALPulsar 7.1.1.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
pulsar_crosscorr.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2007 Badri Krishnan
3 * Copyright (C) 2008 Christine Chung, Badri Krishnan and John Whelan
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with with program; see the file COPYING. If not, write to the
17 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18 * MA 02110-1301 USA
19 */
20
21/**
22 * \author Christine Chung, Badri Krishnan, John Whelan
23 * \date 2008
24 * \file
25 * \ingroup lalpulsar_bin_CrossCorr
26 * \brief Perform CW cross-correlation search
27 *
28 * Id: pulsar_crosscorr.c,v 1.23 2009/03/13 00:43:04 cchung Exp
29 *
30 */
31
32
33#include "config.h"
34
35#include <gsl/gsl_permutation.h>
36
37#include <lal/PulsarCrossCorr.h>
38#include <lal/DopplerScan.h>
39#include <lal/ExtrapolatePulsarSpins.h>
40#include <lal/LALString.h>
41#include <lal/LALPulsarVCSInfo.h>
42
43#include "pulsar_crosscorr.h"
44
45/* globals, constants and defaults */
46
47
48/* user input variables */
54
61REAL8 uvar_dAlpha, uvar_dDelta; /* resolution for isotropic sky-grid */
76
77
78CHAR *uvar_ephemEarth; /**< Earth ephemeris file to use */
79CHAR *uvar_ephemSun; /**< Sun ephemeris file to use */
80
81
88
89#define F0 100
90#define FBAND 1
91
92#define BLOCKSRNGMED 51
93#define MAXFILENAMELENGTH 512 /* maximum # of characters of a filename */
94
95#define DIROUT "./output/" /* output directory */
96#define FILEOUT "CrossCorr_out.dat"
97#define DEBUGOUT "estimator.dat"
98#define BASENAMEOUT "radio" /* prefix file output */
99
100#define SKYFILE "./skypatchfile"
101#define SKYREGION "allsky"
102
103#define TRUE (1==1)
104#define FALSE (1==0)
105
106#define SQUARE(x) ((x)*(x))
107#define CUBE(x) ((x)*(x)*(x))
108
109#define N_SPINDOWN_DERIVS 6
110
112
113
114int main( int argc, char *argv[] )
115{
116 /* LALStatus pointer */
117 static LALStatus status;
118
119 /* sft related variables */
120
121 SFTVector *inputSFTs = NULL;
122 SFTtype *tmpSFT = NULL;
123 REAL8FrequencySeries *tmpPSD = NULL;
124
125 SFTListElement *sftList, *sftHead = NULL, *sftTail = NULL;
126 PSDListElement *psdList, *psdHead = NULL, *psdTail = NULL;
127 REAL8ListElement *freqList, *freqHead = NULL, *freqTail = NULL;
128 REAL8ListElement *phaseList, *phaseHead = NULL, *phaseTail = NULL;
129 REAL8 deltaF_SFT, timeBase;
130 REAL8 *psi = NULL;
131 INT8 counter = 0;
132 COMPLEX8FrequencySeries *sft1 = NULL, *sft2 = NULL;
133 INT4 sameDet;
134 REAL8FrequencySeries *psd1, *psd2;
135 COMPLEX16Vector *yalpha = NULL, *ualpha = NULL;
136
137 /*estimator arrays */
138 COMPLEX16Vector *gplus = NULL, *gcross = NULL;
139 static REAL8Vector *aplussq1, *aplussq2, *acrossq1, *acrossq2;
140 REAL8Vector *galphasq, *galphare, *galphaim;
141 INT4 i;
142
143 CrossCorrAmps amplitudes;
144 CrossCorrBeamFn *beamfns1, *beamfns2;
145 CrossCorrBeamFnListElement *beamList, *beamHead = NULL, *beamTail = NULL;
146
147
148 /* information about all the ifos */
149 DetChoice detChoice;
150 LIGOTimeGPS firstTimeStamp, lastTimeStamp;
151 REAL8 tObs;
152
153 /* ephemeris */
154 EphemerisData *edat = NULL;
155
156 /* skypatch info */
157 REAL8 *skyAlpha = NULL, *skyDelta = NULL,
158 *skySizeAlpha = NULL, *skySizeDelta = NULL;
159 INT8 nSkyPatches, skyCounter = 0;
160 SkyPatchesInfo skyInfo;
161
162 /* frequency loop info */
163 INT8 nfreqLoops = 1, freqCounter = 0;
164 INT8 nParams = 0;
165 REAL8 f_current = 0.0;
166 INT8 ualphacounter = 0.0;
167
168 /* frequency derivative loop info. we can go up to f_doubledot */
169 INT8 nfdotLoops = 1, fdotCounter = 0;
170 INT8 nfddotLoops = 1, fddotCounter = 0;
171 REAL8 fdot_current = 0.0, delta_fdot = 0.0;
172 REAL8 fddot_current = 0.0, delta_fddot = 0.0;
173
174 /* frequency derivative array, if we search over q1, q2, n */
175 INT8 nq1Loops = 1, nq2Loops = 1, nnLoops = 1;
176 INT8 q1Counter = 0, q2Counter = 0, nCounter = 0;
177 REAL8 q1_current = 0.0, q2_current = 0.0, n_current = 0.0;
178 REAL8 delta_q1 = 0.0, delta_q2 = 0.0, delta_n = 0.0;
179 REAL8Vector *fdots = NULL;
180
181 INT8 paramCounter = 0;
182
183 REAL8Vector *sigmasq;
184 PulsarDopplerParams thisPoint;
185 static REAL8Vector *rho, *variance;
186 REAL8 tmpstat, freq1, fbin1, phase1, freq2, fbin2, phase2;
187 UINT4 bin1, bin2;
188 REAL8 tmpstat2, tmpstat3, tmpstat4;
189
190 REAL8 doppWings, fMin, fMax;
191
192
193 /* output files */
194 FILE *fpdebug = NULL;
195 FILE *fpCrossCorr = NULL;
196 CHAR *fnameDebug = NULL;
197 CHAR *fnameCrossCorrCand = NULL;
198
199 /* sft constraint variables */
200 LIGOTimeGPS startTimeGPS, endTimeGPS, refTime;
201
202 FILE *skytest = NULL;
203
204 SkyPosition skypos;
205
206 /* new SFT I/O data types */
207 SFTCatalog *catalog = NULL;
208 static SFTConstraints constraints;
209 INT4 sftcounter = 0, slidingcounter = 1, listLength = 0;
210
211 /* to time the code*/
212 time_t t1, t2;
213
214 /* LAL error-handler */
216
217
219
220 /* read all command line variables */
221 BOOLEAN should_exit = 0;
223 if ( should_exit ) {
224 exit( 1 );
225 }
226
227 /* very basic consistency checks on user input */
228 if ( uvar_f0 < 0 ) {
229 fprintf( stderr, "start frequency must be positive\n" );
230 exit( 1 );
231 }
232
233 if ( uvar_fBand < 0 ) {
234 fprintf( stderr, "search frequency band must be positive\n" );
235 exit( 1 );
236 }
237
239 fprintf( stderr, "search frequency resolution must be > 0\n" );
240 exit( 1 );
241 }
242
244 fprintf( stderr, "fdot and fddot are not used if useQCoeffs is set\n" );
245 exit( 1 );
246 }
247
249 fprintf( stderr, "useQCoeffs must be set in order to search over q1, q2 or braking index\n" );
250 exit( 1 );
251 }
252
253 if ( uvar_fdotBand < 0 ) {
254 fprintf( stderr, "frequency derivative band must be positive\n" );
255 exit( 1 );
256 }
257
259 fprintf( stderr, "frequency derivative resolution must be > 0\n" );
260 exit( 1 );
261 }
262
263 if ( uvar_fddotBand < 0 ) {
264 fprintf( stderr, "frequency double derivative band must be positive\n" );
265 exit( 1 );
266 }
267
268 if ( uvar_fRef < 0 ) {
269 fprintf( stderr, "reference frequency must be positive\n" );
270 exit( 1 );
271 }
272
274 fprintf( stderr, "frequency double derivative resolution must be > 0\n" );
275 exit( 1 );
276 }
277
279 fprintf( stderr, "if uvar_averagePsi = TRUE, psi will not be used\n" );
280 exit( 1 );
281 }
282
284 fprintf( stderr, "if uvar_averageIota = TRUE, cosi will not be used\n" );
285 exit( 1 );
286 }
287
288 if ( uvar_autoCorrelate ) {
289 fprintf( stderr, "autoCorrelate option is not currently usable\n" );
290 exit( 1 );
291 }
292
293
294
295 /* if debugging, and averaging over both psi & iota, OR using exact psi & iota
296 * values, print out estimator info */
299
300 fnameDebug = LALCalloc( strlen( uvar_debugOut ) + strlen( uvar_dirnameOut ) + 1, sizeof( CHAR ) );
301 if ( fnameDebug == NULL ) {
302 fprintf( stderr, "error allocating memory [pulsar_crosscorr.c %d]\n", __LINE__ );
303 return ( PULSAR_CROSSCORR_EMEM );
304 }
305
306 strcpy( fnameDebug, uvar_dirnameOut );
307 strcat( fnameDebug, uvar_debugOut );
308 if ( !( fpdebug = fopen( fnameDebug, "wb" ) ) ) {
309 fprintf( stderr, "Unable to open debugging output file '%s' for writing.\n", fnameDebug );
311 }
312
313 }
314
315 /* initialise output file name */
316
317 fnameCrossCorrCand = LALCalloc( strlen( uvar_filenameOut ) + strlen( uvar_dirnameOut ) + 1, sizeof( CHAR ) );
318 if ( fnameCrossCorrCand == NULL ) {
319 fprintf( stderr, "error allocating memory [pulsar_crosscorr.c %d]\n", __LINE__ );
320 return ( PULSAR_CROSSCORR_EMEM );
321 }
322
323 strcpy( fnameCrossCorrCand, uvar_dirnameOut );
324 strcat( fnameCrossCorrCand, uvar_filenameOut );
325
326 if ( !( fpCrossCorr = fopen( fnameCrossCorrCand, "wb" ) ) ) {
327 fprintf( stderr, "Unable to open output-file '%s' for writing.\n", fnameCrossCorrCand );
329 }
330
331
332 if ( uvar_QCoeffs ) {
333 fprintf( fpCrossCorr, "##Alpha\tDelta\tFrequency\tQ1\tQ2\tBraking Index\tNormalised Power\n" );
334 } else {
335 fprintf( fpCrossCorr, "##Alpha\tDelta\tFrequency\t Fdot \t Fddot \t Normalised Power\n" );
336 }
337
338
339
340 /* set sft catalog constraints */
341 constraints.detector = NULL;
342 constraints.timestamps = NULL;
343 constraints.minStartTime = NULL;
344 constraints.maxStartTime = NULL;
345
346 XLALGPSSet( &startTimeGPS, 0, 0 );
347 XLALGPSSet( &endTimeGPS, LAL_INT4_MAX, 0 );
348
350 XLALGPSSetREAL8( &startTimeGPS, uvar_startTime );
351
352 constraints.minStartTime = &startTimeGPS;
353 }
354
356 XLALGPSSetREAL8( &endTimeGPS, uvar_endTime );
357 constraints.maxStartTime = &endTimeGPS;
358 }
359
360 /* get sft catalog */
361 /* note that this code depends very heavily on the fact that the catalog
362 returned by XLALSFTdataFind is time sorted */
363
364 time( &t1 );
365
366 XLAL_CHECK_MAIN( ( catalog = XLALSFTdataFind( uvar_sftDir, &constraints ) ) != NULL, XLAL_EFUNC );
367 if ( ( catalog == NULL ) || ( catalog->length == 0 ) ) {
368 fprintf( stderr, "Unable to match any SFTs with pattern '%s'\n",
369 uvar_sftDir );
370 exit( 1 );
371 }
372
373 time( &t2 );
374
375 if ( uvar_timingOn ) {
376 fprintf( stderr, "Time taken to load sft catalog: %f s\n", difftime( t2, t1 ) );
377 }
378
379 /* get SFT parameters so that we can initialise search frequency resolutions */
380 /* calculate deltaF_SFT */
381 deltaF_SFT = catalog->data[0].header.deltaF; /* frequency resolution */
382 timeBase = 1.0 / deltaF_SFT; /* sft baseline */
383
384 /* catalog is ordered in time so we can get start, end time and tObs */
385 firstTimeStamp = catalog->data[0].header.epoch;
386 lastTimeStamp = catalog->data[catalog->length - 1].header.epoch;
387 tObs = XLALGPSDiff( &lastTimeStamp, &firstTimeStamp ) + timeBase;
388
389 /*set pulsar reference time */
391 XLALGPSSetREAL8( &refTime, uvar_refTime );
392 } else { /*if refTime is not set, set it to midpoint of sfts*/
393 XLALGPSSetREAL8( &refTime, ( 0.5 * tObs ) + XLALGPSGetREAL8( &firstTimeStamp ) );
394 }
395
396 /* set frequency resolution defaults if not set by user */
397 if ( !( XLALUserVarWasSet( &uvar_fResolution ) ) ) {
398 uvar_fResolution = 1 / tObs;
399 }
400
401 /*get number of frequency loops*/
402 nfreqLoops = rint( uvar_fBand / uvar_fResolution );
403 /* if we are using spindown parameters, initialise the fdots array */
404 if ( uvar_QCoeffs ) {
405
406 fdots = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
407 fdots->length = N_SPINDOWN_DERIVS;
408 fdots->data = ( REAL8 * )LALCalloc( fdots->length, sizeof( REAL8 ) );
409
410 if ( uvar_q1Band > 0 ) {
411 nq1Loops = rint( uvar_q1Band / uvar_q1Resolution );
412 }
413
414 if ( uvar_q2Band > 0 ) {
415 nq2Loops = rint( uvar_q2Band / uvar_q2Resolution );
416 }
417
418 if ( uvar_brakingindexBand > 0 ) {
420 }
421
422 delta_q1 = uvar_q1Resolution;
423
424 delta_q2 = uvar_q2Resolution;
425
427 }
428 /*otherwise just search over f0, fdot, fddot */
429 else {
430
431 if ( !( XLALUserVarWasSet( &uvar_fdotResolution ) ) ) {
432 uvar_fdotResolution = SQUARE( 1 / tObs );
433 }
434
436 uvar_fddotResolution = CUBE( 1 / tObs );
437 }
438
439 if ( uvar_fdotBand > 0 ) {
440 nfdotLoops = rint( uvar_fdotBand / uvar_fdotResolution );
441 }
442
443 if ( uvar_fddotBand > 0 ) {
444 nfddotLoops = rint( uvar_fddotBand / uvar_fddotResolution );
445 }
446
447 delta_fdot = uvar_fdotResolution;
448
449 delta_fddot = uvar_fddotResolution;
450 }
451
452 /* set up ephemeris */
454
455 /* set up skypatches */
456 if ( ( skytest = fopen( uvar_skyfile, "r" ) ) == NULL ) {
457 fprintf( stderr, "skyfile doesn't exist\n" );
458 }
459
463 &status );
464 nSkyPatches = skyInfo.numSkyPatches;
465 skyAlpha = skyInfo.alpha;
466 skyDelta = skyInfo.delta;
467 skySizeAlpha = skyInfo.alphaSize;
468 skySizeDelta = skyInfo.deltaSize;
469
470
471 /* curly As */
472 /* because we have the option of either averaging over i or not, we
473 need to calculate A_{+,x}^2 and A_xA_+ rather than the individual
474 values because <A_x> = 0 */
475 if ( uvar_averageIota ) {
476 amplitudes.Aplussq = 7.0 / 15.0;
477 amplitudes.Acrosssq = 1.0 / 3.0;
478 amplitudes.AplusAcross = 0;
479 } else {
480 amplitudes.Aplussq = SQUARE( ( 1.0 + uvar_cosi * uvar_cosi ) / 2.0 );
481 amplitudes.Acrosssq = SQUARE( uvar_cosi );
482 amplitudes.AplusAcross = ( uvar_cosi / 2.0 ) + ( CUBE( uvar_cosi ) / 2.0 );
483 }
484
485 /* polarisation angle */
486 if ( XLALUserVarWasSet( &uvar_psi ) ) {
487 psi = ( REAL8 * ) LALCalloc( 1, sizeof( REAL8 ) );
488 *psi = uvar_psi;
489 } else {
490 psi = NULL;
491 }
492
493 /* initialise output arrays */
494 if ( uvar_QCoeffs ) {
495 nParams = nSkyPatches * nfreqLoops * nq1Loops * nq2Loops * nnLoops;
496 } else {
497 nParams = nSkyPatches * nfreqLoops * nfdotLoops * nfddotLoops;
498 }
499 rho = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
500 rho->length = nParams;
501 rho->data = ( REAL8 * )LALCalloc( nParams, sizeof( REAL8 ) );
502
503 variance = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
504 variance->length = nParams;
505 variance->data = ( REAL8 * )LALCalloc( nParams, sizeof( REAL8 ) );
506
507 /*initialise debugging vectors */
508 aplussq1 = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
509 aplussq1->length = nParams;
510 aplussq1->data = ( REAL8 * )LALCalloc( nParams, sizeof( REAL8 ) );
511
512 aplussq2 = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
513 aplussq2->length = nParams;
514 aplussq2->data = ( REAL8 * )LALCalloc( nParams, sizeof( REAL8 ) );
515
516 acrossq1 = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
517 acrossq1->length = nParams;
518 acrossq1->data = ( REAL8 * )LALCalloc( nParams, sizeof( REAL8 ) );
519
520 acrossq2 = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
521 acrossq2->length = nParams;
522 acrossq2->data = ( REAL8 * )LALCalloc( nParams, sizeof( REAL8 ) );
523
524 galphasq = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
525 galphasq->length = nParams;
526 galphasq->data = ( REAL8 * )LALCalloc( nParams, sizeof( REAL8 ) );
527
528 galphare = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
529 galphare->length = nParams;
530 galphare->data = ( REAL8 * )LALCalloc( nParams, sizeof( REAL8 ) );
531
532 galphaim = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
533 galphaim->length = nParams;
534 galphaim->data = ( REAL8 * )LALCalloc( nParams, sizeof( REAL8 ) );
535
536 /*initialise detector choice*/
537 detChoice = uvar_detChoice;
538
539
540 {
541 /* block for calculating frequency range to read from SFTs */
542 /* user specifies freq and fdot range at reftime
543 we translate this range of fdots to start and endtime and find
544 the largest frequency band required to cover the
545 frequency evolution */
546 PulsarSpinRange spinRange_startTime; /**< freq and fdot range at start-time of observation */
547 PulsarSpinRange spinRange_endTime; /**< freq and fdot range at end-time of observation */
548 PulsarSpinRange spinRange_refTime; /**< freq and fdot range at the reference time */
549
550 REAL8 startTime_freqLo, startTime_freqHi, endTime_freqLo, endTime_freqHi, freqLo, freqHi;
551
552 REAL8Vector *fdotsMin = NULL;
553 REAL8Vector *fdotsMax = NULL;
554
555 UINT4 k;
556
557 fdotsMin = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
558 fdotsMin->length = N_SPINDOWN_DERIVS;
559 fdotsMin->data = ( REAL8 * )LALCalloc( fdotsMin->length, sizeof( REAL8 ) );
560
561 fdotsMax = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
562 fdotsMax->length = N_SPINDOWN_DERIVS;
563 fdotsMax->data = ( REAL8 * )LALCalloc( fdotsMax->length, sizeof( REAL8 ) );
564
565 XLAL_INIT_MEM( spinRange_startTime );
566 XLAL_INIT_MEM( spinRange_endTime );
567 XLAL_INIT_MEM( spinRange_refTime );
568
569 spinRange_refTime.refTime = refTime;
570 spinRange_refTime.fkdot[0] = uvar_f0;
571 spinRange_refTime.fkdotBand[0] = uvar_fBand;
572
573 /* this assumes that user input parameter ranges such as uvar_fBand are positive */
574 if ( uvar_QCoeffs ) {
575
578
579 /* need to add uvar_xxResolution to account for the fact that we do an extra loop in q1, q2 and n
580 * when uvar_xxBand is greater than 0 */
584
585 for ( k = 1; k < fdotsMin->length; k++ ) {
586 spinRange_refTime.fkdot[k] = fdotsMin->data[k - 1];
587 spinRange_refTime.fkdotBand[k] = fdotsMax->data[k - 1] - fdotsMin->data[k - 1];
588 }
589
590 } else {
591 spinRange_refTime.fkdot[1] = uvar_fdot;
592 spinRange_refTime.fkdotBand[1] = uvar_fdotBand;
593
594 spinRange_refTime.fkdot[2] = uvar_fddot;
595 spinRange_refTime.fkdotBand[2] = uvar_fddotBand;
596 }
597
598 XLAL_CHECK_MAIN( XLALExtrapolatePulsarSpinRange( &spinRange_startTime, &spinRange_refTime, XLALGPSDiff( &firstTimeStamp, &spinRange_refTime.refTime ) ) == XLAL_SUCCESS, XLAL_EFUNC );
599 XLAL_CHECK_MAIN( XLALExtrapolatePulsarSpinRange( &spinRange_endTime, &spinRange_refTime, XLALGPSDiff( &lastTimeStamp, &spinRange_refTime.refTime ) ) == XLAL_SUCCESS, XLAL_EFUNC );
600
601 startTime_freqLo = spinRange_startTime.fkdot[0]; /* lowest search freq at start time */
602 startTime_freqHi = startTime_freqLo + spinRange_startTime.fkdotBand[0]; /* highest search freq. at start time*/
603 endTime_freqLo = spinRange_endTime.fkdot[0];
604 endTime_freqHi = endTime_freqLo + spinRange_endTime.fkdotBand[0];
605
606 /* freqLo = min of low frequencies and freqHi = max of high frequencies */
607 freqLo = startTime_freqLo < endTime_freqLo ? startTime_freqLo : endTime_freqLo;
608 freqHi = startTime_freqHi > endTime_freqHi ? startTime_freqHi : endTime_freqHi;
609
610 /* add wings for Doppler modulation and running median block size */
611 /* remove fBand from doppWings because we are going bin-by-bin (?) */
612
613 doppWings = freqHi * VTOT ;
614 fMin = freqLo - doppWings - uvar_blocksRngMed * deltaF_SFT;
615 fMax = freqHi + doppWings + uvar_blocksRngMed * deltaF_SFT;
616
617 XLALDestroyREAL8Vector( fdotsMin );
618 XLALDestroyREAL8Vector( fdotsMax );
619 }
620
621
622 slidingcounter = 0;
623
624 time( &t1 );
625
626 fprintf( stderr, "beginning main calculations over %" LAL_INT8_FORMAT " loops:\n%" LAL_INT8_FORMAT " freq, %" LAL_INT8_FORMAT " Q1, %" LAL_INT8_FORMAT "n\n", nParams, nfreqLoops, nq1Loops, nnLoops );
627 /***********start main calculations**************/
628 /*outer loop over all sfts in catalog, so that we load only the relevant sfts each time*/
629 for ( sftcounter = 0; sftcounter < ( INT4 )catalog->length - 1; sftcounter++ ) {
630 tmpSFT = NULL;
631 tmpPSD = NULL;
632 yalpha = NULL;
633 ualpha = NULL;
634 sigmasq = NULL;
635 gplus = NULL;
636 gcross = NULL;
637
638 /* reset all search parameters */
639 paramCounter = 0;
640 counter = 0;
641 nCounter = 0;
642 q1Counter = 0;
643 q2Counter = 0;
644 freqCounter = 0;
645 fdotCounter = 0;
646 fddotCounter = 0;
647 skyCounter = 0;
648
649 /* throw away first sft from inputSFTs, and first psdvector, frequency, phase vectors, beamfns */
650 if ( sftcounter > 0 ) {
651
652 LAL_CALL( DeleteSFTHead( &status, &sftHead ), &status );
653
654 LAL_CALL( DeletePSDHead( &status, &psdHead ), &status );
655
656 LAL_CALL( DeleteREAL8Head( &status, &freqHead ), &status );
657
658 LAL_CALL( DeleteREAL8Head( &status, &phaseHead ), &status );
659
660 LAL_CALL( DeleteBeamFnHead( &status, &beamHead ), &status );
661
662
663 }
664 /* make a second sft catalog with only the ones within maxlag of the current sft*/
665 /* do all sfts with same time together */
666 while ( ( slidingcounter < ( INT4 )catalog->length ) &&
667 ( XLALGPSDiff( &catalog->data[slidingcounter].header.epoch, &catalog->data[sftcounter].header.epoch ) <= uvar_maxlag ) ) {
668
669 inputSFTs = NULL;
670
671 LAL_CALL( CopySFTFromCatalog( &status, catalog, &inputSFTs, fMin, fMax, slidingcounter ), &status );
672
673 LAL_CALL( AddSFTtoList( &status, &sftHead, &sftTail, &( inputSFTs->data[0] ) ), &status );
674
675 LAL_CALL( AddPSDtoList( &status, &psdHead, &psdTail, sftTail->sft.data->length ), &status );
676
677 tmpSFT = &( sftTail->sft );
678 tmpPSD = &( psdTail->psd );
679
681 tmpSFT, uvar_blocksRngMed, 0.0 ) == XLAL_SUCCESS, XLAL_EFUNC );
682
683 LAL_CALL( AddREAL8toList( &status, &freqHead, &freqTail ), &status );
684
685 LAL_CALL( AddREAL8toList( &status, &phaseHead, &phaseTail ), &status );
686
687 LAL_CALL( AddBeamFntoList( &status, &beamHead, &beamTail ), &status );
688
689 slidingcounter++;
690 XLALDestroySFTVector( inputSFTs );
691
692 }
693
694 listLength = slidingcounter - sftcounter;
695 if ( listLength > 1 ) {
696
697 while ( paramCounter < nParams ) {
698 f_current = uvar_f0 + ( uvar_fResolution * freqCounter );
699
700 if ( uvar_QCoeffs ) {
701
702 q1_current = uvar_q1 + ( delta_q1 * q1Counter );
703 q2_current = uvar_q2 + ( delta_q2 * q2Counter );
704 n_current = uvar_brakingindex + ( delta_n * nCounter );
705
706 skyCounter++;
707 if ( skyCounter == nSkyPatches ) {
708 skyCounter = 0;
709 nCounter++;
710 }
711 if ( nCounter == nnLoops ) {
712 nCounter = 0;
713 q2Counter++;
714 }
715 if ( q2Counter == nq2Loops ) {
716 q2Counter = 0;
717 q1Counter++;
718 }
719 if ( q1Counter == nq1Loops ) {
720 q1Counter = 0;
721 freqCounter++;
722 }
723
724
725 } else {
726
727 fdot_current = uvar_fdot + ( delta_fdot * fdotCounter );
728 fddot_current = uvar_fddot + ( delta_fddot * fddotCounter );
729
730 skyCounter++;
731 if ( skyCounter == nSkyPatches ) {
732 skyCounter = 0;
733 fddotCounter++;
734 }
735 if ( fddotCounter == nfddotLoops ) {
736 fddotCounter = 0;
737 fdotCounter++;
738 }
739 if ( fdotCounter == nfdotLoops ) {
740 fdotCounter = 0;
741 freqCounter++;
742 }
743
744 }
745
746 LAL_CALL( InitDoppParams( &status, fdots, &thisPoint, refTime, f_current, q1_current, q2_current, n_current,
747 fdot_current, fddot_current ), &status );
748
749 /* set sky positions and skypatch sizes */
750 thisPoint.Alpha = skyAlpha[skyCounter];
751 thisPoint.Delta = skyDelta[skyCounter];
752
753 /* get the amplitude modulation coefficients */
754 skypos.longitude = thisPoint.Alpha;
755 skypos.latitude = thisPoint.Delta;
757
758
759 LAL_CALL( GetBeamInfo( &status, beamHead, sftHead, freqHead, phaseHead, skypos,
760 edat, &thisPoint ), &status );
761
762 /* loop over SFT mini-list to get pairs */
763 ualphacounter = 0;
764
765 /* correlate sft pairs */
766 sftList = sftHead;
767 psdList = psdHead;
768 freqList = freqHead;
769 phaseList = phaseHead;
770 beamList = beamHead;
771
772 sft1 = &( sftList->sft );
773 psd1 = &( psdList->psd );
774 freq1 = freqList->val;
775 bin1 = ( UINT4 )ceil( ( ( freq1 - psd1->f0 ) / ( deltaF_SFT ) ) - 0.5 );
776 fbin1 = psd1->f0 + ( REAL8 )bin1 * deltaF_SFT;
777 phase1 = phaseList->val;
778 beamfns1 = &( beamList->beamfn );
779
780 /*while there are elements in the sft minilist, keep
781 going and check whether it should be paired with SFT1.
782 there is no need to check the lag as the sfts must satisfy this condition
783 already to be in the mini-list*/
784 while ( sftList->nextSFT ) {
785
786 /*if we are autocorrelating, we want the head to be paired with itself first*/
787 if ( ( sftList == sftHead ) && uvar_autoCorrelate ) {
788 sft2 = &( sftList->sft );
789 psd2 = &( psdList->psd );
790 freq2 = freqList->val;
791 bin2 = ( UINT4 )ceil( ( ( freq2 - psd2->f0 ) / ( deltaF_SFT ) ) - 0.5 );
792 fbin2 = psd2->f0 + ( REAL8 )bin2 * deltaF_SFT;
793 phase2 = phaseList->val;
794 beamfns2 = &( beamList->beamfn );
795
796 sftList = ( SFTListElement * )sftList->nextSFT;
797 psdList = ( PSDListElement * )psdList->nextPSD;
798 freqList = ( REAL8ListElement * )freqList->nextVal;
799 phaseList = ( REAL8ListElement * )phaseList->nextVal;
800 beamList = ( CrossCorrBeamFnListElement * )beamList->nextBeamfn;
801
802 } else { /*otherwise just step to the next sft*/
803
804 sftList = ( SFTListElement * )sftList->nextSFT;
805 psdList = ( PSDListElement * )psdList->nextPSD;
806 freqList = ( REAL8ListElement * )freqList->nextVal;
807 phaseList = ( REAL8ListElement * )phaseList->nextVal;
808 beamList = ( CrossCorrBeamFnListElement * )beamList->nextBeamfn;
809
810 sft2 = &( sftList->sft );
811 psd2 = &( psdList->psd );
812 freq2 = freqList->val;
813 bin2 = ( UINT4 )ceil( ( ( freq2 - psd2->f0 ) / ( deltaF_SFT ) ) - 0.5 );
814 fbin2 = psd2->f0 + ( REAL8 )bin2 * deltaF_SFT;
815 phase2 = phaseList->val;
816 beamfns2 = &( beamList->beamfn );
817 }
818 /*strcmp returns 0 if strings are equal, >0 if strings are different*/
819 sameDet = strcmp( sft1->name, sft2->name );
820 /* if they are different, set sameDet to 1 so that it will match if
821 detChoice == DIFFERENT */
822 if ( sameDet != 0 ) {
823 sameDet = 1;
824 }
825
826 /* however, if detChoice = ALL, then we want sameDet to match it */
827 if ( detChoice == ALL ) {
828 sameDet = detChoice;
829 }
830
831 /* decide whether to add this pair or not */
832 if ( sameDet == ( INT4 )detChoice ) {
833
834 /* increment the size of Y, u, sigmasq vectors by 1 */
835 yalpha = XLALResizeCOMPLEX16Vector( yalpha, 1 + ualphacounter );
836 ualpha = XLALResizeCOMPLEX16Vector( ualpha, 1 + ualphacounter );
837 sigmasq = XLALResizeREAL8Vector( sigmasq, 1 + ualphacounter );
838 gplus = XLALResizeCOMPLEX16Vector( gplus, 1 + ualphacounter );
839 gcross = XLALResizeCOMPLEX16Vector( gcross, 1 + ualphacounter );
840
841 LAL_CALL( LALCorrelateSingleSFTPair( &status, &( yalpha->data[ualphacounter] ),
842 sft1, sft2, psd1, psd2, bin1, bin2 ),
843 &status );
844
845 LAL_CALL( LALCalculateSigmaAlphaSq( &status, &sigmasq->data[ualphacounter],
846 bin1, bin2, psd1, psd2 ),
847 &status );
848 /*if we are averaging over psi and cos(iota), call the simplified
849 Ualpha function*/
851 LAL_CALL( LALCalculateAveUalpha( &status, &ualpha->data[ualphacounter],
852 phase1, phase2, fbin1, fbin2, deltaF_SFT, *beamfns1, *beamfns2,
853 sigmasq->data[ualphacounter] ),
854 &status );
855
856 } else {
857 LAL_CALL( LALCalculateUalpha( &status, &ualpha->data[ualphacounter], amplitudes,
858 phase1, phase2, fbin1, fbin2, deltaF_SFT, *beamfns1, *beamfns2,
859 sigmasq->data[ualphacounter], psi, &gplus->data[ualphacounter], &gcross->data[ualphacounter] ),
860 &status );
861 }
862 ualphacounter++;
863
864 }
865 } /*finish loop over sft pairs*/
866
867 /* calculate rho from Yalpha and Ualpha, if there were pairs */
868 if ( ualphacounter > 0 ) {
869 tmpstat = 0;
870 tmpstat2 = 0;
871 tmpstat3 = 0;
872 tmpstat4 = 0;
873 LAL_CALL( LALCalculateCrossCorrPower( &status, &tmpstat, yalpha, ualpha ),
874 &status );
875
876 rho->data[counter] += tmpstat;
877
878 /* calculate standard deviation of rho (Eq 4.6) */
879 LAL_CALL( LALNormaliseCrossCorrPower( &status, &tmpstat, ualpha, sigmasq ),
880 &status );
881
882 variance->data[counter] += tmpstat;
883
884 /*
885 for (i=0; i < (INT4)ualpha->length; i++) {
886 printf("%g %g\n", sigmasq->data[i] * ualpha->data[i].re, sigmasq->data[i] * ualpha->data[i].im);
887 }
888 */
889
890
892 LAL_CALL( LALCalculateEstimators( &status, &tmpstat, &tmpstat2, &tmpstat3, &tmpstat4, yalpha, gplus, gcross, sigmasq ), &status );
893
894 aplussq1->data[counter] += tmpstat;
895 aplussq2->data[counter] += tmpstat2;
896 acrossq1->data[counter] += tmpstat3;
897 acrossq2->data[counter] += tmpstat4;
898
899 for ( i = 0; i < ( INT4 )ualpha->length; i++ ) {
900
901 galphasq->data[counter] += SQUARE( sigmasq->data[i] * creal( ualpha->data[i] ) ) + SQUARE( sigmasq->data[i] * cimag( ualpha->data[i] ) );
902
903 galphare->data[counter] += ( sigmasq->data[i] * creal( ualpha->data[i] ) );
904 galphaim->data[counter] += -( sigmasq->data[i] * cimag( ualpha->data[i] ) );
905
906 }
907
908 }
909 } /* endif (ualphacounter > 0) */
910
912
913 for ( i = 0; i < ( INT4 )ualpha->length; i++ ) {
914
915 galphasq->data[counter] += ( SQUARE( sigmasq->data[i] * creal( ualpha->data[i] ) ) + SQUARE( sigmasq->data[i] * cimag( ualpha->data[i] ) ) );
916 galphare->data[counter] += ( sigmasq->data[i] * creal( ualpha->data[i] ) );
917 galphaim->data[counter] += -( sigmasq->data[i] * cimag( ualpha->data[i] ) );
918
919 }
920 }
921 counter++;
922 paramCounter++;
923
924 } /*endwhile*/
925
930
931 XLALDestroyREAL8Vector( sigmasq );
932
933
934
935
936 } /*end if listLength > 1 */
937 } /* finish loop over all sfts */
938
939 fprintf( stderr, "finish loop over all sfts\n" );
940
941 time( &t2 );
942
943 if ( uvar_timingOn ) {
944 fprintf( stderr, "Time taken for main loop: %f\n", difftime( t2, t1 ) );
945 }
946
947 counter = 0;
948
949 time( &t1 );
950 /* print output - all variables to file */
951
952 for ( freqCounter = 0; freqCounter < nfreqLoops; freqCounter++ ) {
953
954 f_current = uvar_f0 + ( uvar_fResolution * freqCounter );
955
956 if ( uvar_QCoeffs ) { /*if searching over q1, q2, n*/
957
958 /* Q1 loop */
959 for ( q1Counter = 0; q1Counter < nq1Loops; q1Counter++ ) {
960
961 q1_current = uvar_q1 + ( delta_q1 * q1Counter );
962
963 /* Q2 loop */
964 for ( q2Counter = 0; q2Counter < nq2Loops; q2Counter++ ) {
965
966 q2_current = uvar_q2 + ( delta_q2 * q2Counter );
967
968 /* n loop */
969 for ( nCounter = 0; nCounter < nnLoops; nCounter++ ) {
970
971 n_current = uvar_brakingindex + ( delta_n * nCounter );
972
973 for ( skyCounter = 0; skyCounter < nSkyPatches; skyCounter++ ) {
974 /* initialize Doppler parameters of the potential source */
975 thisPoint.Alpha = skyAlpha[skyCounter];
976 thisPoint.Delta = skyDelta[skyCounter];
977
978 /*normalise rho by stddev */
979//printf("raw rho %1.15g\n", rho->data[counter]);
980 rho->data[counter] = rho->data[counter] / sqrt( variance->data[counter] );
981 fprintf( fpCrossCorr, "%1.5f\t %1.5f\t %1.5f\t %e\t %e\t %e\t %1.10g\n", thisPoint.Alpha,
982 thisPoint.Delta, f_current,
983 q1_current, q2_current, n_current, rho->data[counter] );
984
986 fprintf( fpdebug, "%1.5f %e %e %g %g %g\n", f_current, sqrt( fabs( aplussq2->data[counter] / aplussq1->data[counter] ) ), sqrt( fabs( acrossq2->data[counter] / acrossq1->data[counter] ) ), galphasq->data[counter], galphare->data[counter], galphaim->data[counter] );
987 }
988
990 fprintf( fpdebug, "%1.5f %g %g %g\n", f_current, galphasq->data[counter], galphare->data[counter], galphaim->data[counter] );
991 }
992
993 counter++;
994 }
995 } /*end n loop*/
996 } /*end q2loop*/
997 } /*end q1 loop*/
998 } /*endif uvar_Qcoeffs */
999
1000 else {
1001
1002 /* frequency derivative loop */
1003 for ( fdotCounter = 0; fdotCounter < nfdotLoops; fdotCounter++ ) {
1004
1005 fdot_current = uvar_fdot + ( delta_fdot * fdotCounter );
1006
1007 /* frequency double derivative loop */
1008 for ( fddotCounter = 0; fddotCounter < nfddotLoops; fddotCounter++ ) {
1009 for ( skyCounter = 0; skyCounter < nSkyPatches; skyCounter++ ) {
1010 /* initialize Doppler parameters of the potential source */
1011 thisPoint.Alpha = skyAlpha[skyCounter];
1012 thisPoint.Delta = skyDelta[skyCounter];
1013 /*normalise rho*/
1014 rho->data[counter] = rho->data[counter] / sqrt( variance->data[counter] );
1015 fprintf( fpCrossCorr, "%1.5f\t %1.5f\t %1.5f\t %e\t %e\t %1.10f\n", thisPoint.Alpha,
1016 thisPoint.Delta, f_current,
1017 fdot_current, fddot_current, rho->data[counter] );
1018 counter++;
1019 }
1020 }
1021 }
1022 } /*endelse uvar_Qcoeffs*/
1023 }
1024
1025 time( &t2 );
1026 if ( uvar_timingOn ) {
1027 fprintf( stderr, "Time taken to write to output file: %f\n", difftime( t2, t1 ) );
1028 }
1029 /* select candidates */
1030
1031
1032
1033 fclose( fpCrossCorr );
1034
1036
1037 fclose( fpdebug );
1038
1039 }
1040 XLALDestroySFTCatalog( catalog );
1041
1043
1044 LALFree( skyAlpha );
1045 LALFree( skyDelta );
1046 LALFree( skySizeAlpha );
1047 LALFree( skySizeDelta );
1048 LALFree( fnameCrossCorrCand );
1049 LALFree( fnameDebug );
1050 XLALDestroyREAL8Vector( variance );
1052 XLALDestroyREAL8Vector( aplussq1 );
1053 XLALDestroyREAL8Vector( aplussq2 );
1054 XLALDestroyREAL8Vector( acrossq1 );
1055 XLALDestroyREAL8Vector( acrossq2 );
1056 XLALDestroyREAL8Vector( galphasq );
1057 XLALDestroyREAL8Vector( galphare );
1058 XLALDestroyREAL8Vector( galphaim );
1059
1060 if ( !uvar_averagePsi ) {
1061 LALFree( psi );
1062 }
1063
1064 if ( uvar_QCoeffs ) {
1065 XLALDestroyREAL8Vector( fdots );
1066 }
1067
1068 /*free the last few elements (if they're not already free). */
1069 if ( beamHead ) {
1070 if ( beamHead != beamTail ) {
1071 XLALFree( beamTail );
1072 }
1073 XLALFree( beamHead );
1074 }
1075
1076 if ( sftHead ) {
1077 if ( sftHead != sftTail ) {
1078 tmpSFT = &( sftTail->sft );
1079 XLALDestroySFT( tmpSFT );
1080 }
1081 tmpSFT = &( sftHead->sft );
1082 XLALDestroySFT( tmpSFT );
1083 }
1084 if ( psdHead ) {
1085 if ( psdHead != psdTail ) {
1086 XLALDestroyREAL8FrequencySeries( &( psdTail->psd ) );
1087 }
1088
1089 XLALDestroyREAL8FrequencySeries( &( psdHead->psd ) );
1090 }
1091
1092 if ( phaseHead ) {
1093 if ( phaseHead != phaseTail ) {
1094 XLALFree( phaseTail );
1095 }
1096 XLALFree( phaseHead );
1097 }
1098
1099 if ( freqHead ) {
1100 if ( freqHead != freqTail ) {
1101 XLALFree( freqTail );
1102 }
1103 XLALFree( freqHead );
1104 }
1105
1107
1109
1110 if ( lalDebugLevel ) {
1111 REPORTSTATUS( &status );
1112 }
1113
1114 return status.statusCode;
1115
1116} /* main */
1117
1118
1119/**
1120 * Set up location of skypatch centers and sizes
1121 * If user specified skyRegion then use DopplerScan function
1122 * to construct an isotropic grid. Otherwise use skypatch file.
1123 */
1124void SetUpRadiometerSkyPatches( LALStatus *status, /**< pointer to LALStatus structure */
1125 SkyPatchesInfo *out, /**< output skypatches info */
1126 CHAR *skyFileName, /**< name of skypatch file */
1127 CHAR *skyRegion, /**< skyregion (if isotropic grid is to be constructed) */
1128 REAL8 dAlpha, /**< alpha resolution (if isotropic grid is to be constructed) */
1129 REAL8 dDelta ) /**< delta resolution (if isotropic grid is to be constructed) */
1130{
1131
1132 DopplerSkyScanInit XLAL_INIT_DECL( scanInit ); /* init-structure for DopperScanner */
1133 DopplerSkyScanState XLAL_INIT_DECL( thisScan ); /* current state of the Doppler-scan */
1134 UINT4 nSkyPatches, skyCounter;
1135 PulsarDopplerParams dopplerpos;
1136
1137 INITSTATUS( status );
1139
1143
1144 if ( skyRegion ) {
1145
1146 scanInit.dAlpha = dAlpha;
1147 scanInit.dDelta = dDelta;
1148 scanInit.gridType = GRID_ISOTROPIC;
1149 scanInit.metricType = LAL_PMETRIC_NONE;
1150 scanInit.skyRegionString = ( CHAR * )LALCalloc( 1, strlen( skyRegion ) + 1 );
1151 strcpy( scanInit.skyRegionString, skyRegion );
1152 /* scanInit.Freq = usefulParams.spinRange_midTime.fkdot[0] + usefulParams.spinRange_midTime.fkdotBand[0]; */
1153
1154 /* set up the grid */
1155 TRY( InitDopplerSkyScan( status->statusPtr, &thisScan, &scanInit ), status );
1156
1157 nSkyPatches = out->numSkyPatches = thisScan.numSkyGridPoints;
1158 out->alpha = ( REAL8 * )LALCalloc( 1, nSkyPatches * sizeof( REAL8 ) );
1159 out->delta = ( REAL8 * )LALCalloc( 1, nSkyPatches * sizeof( REAL8 ) );
1160 out->alphaSize = ( REAL8 * )LALCalloc( 1, nSkyPatches * sizeof( REAL8 ) );
1161 out->deltaSize = ( REAL8 * )LALCalloc( 1, nSkyPatches * sizeof( REAL8 ) );
1162
1163 /* loop over skygrid points */
1164 XLALNextDopplerSkyPos( &dopplerpos, &thisScan );
1165
1166 skyCounter = 0;
1167 while ( thisScan.state != STATE_FINISHED ) {
1168
1169 out->alpha[skyCounter] = dopplerpos.Alpha;
1170 out->delta[skyCounter] = dopplerpos.Delta;
1171 out->alphaSize[skyCounter] = dAlpha;
1172 out->deltaSize[skyCounter] = dDelta;
1173
1174 if ( ( dopplerpos.Delta > 0 ) && ( dopplerpos.Delta < atan( 4 * LAL_PI / dAlpha / dDelta ) ) ) {
1175 out->alphaSize[skyCounter] = dAlpha * cos( dopplerpos.Delta - 0.5 * dDelta ) / cos( dopplerpos.Delta );
1176 }
1177
1178 if ( ( dopplerpos.Delta < 0 ) && ( dopplerpos.Delta > -atan( 4 * LAL_PI / dAlpha / dDelta ) ) ) {
1179 out->alphaSize[skyCounter] = dAlpha * cos( dopplerpos.Delta + 0.5 * dDelta ) / cos( dopplerpos.Delta );
1180 }
1181
1182 XLALNextDopplerSkyPos( &dopplerpos, &thisScan );
1183 skyCounter++;
1184
1185 } /* end while loop over skygrid */
1186
1187 /* free dopplerscan stuff */
1188 TRY( FreeDopplerSkyScan( status->statusPtr, &thisScan ), status );
1189 if ( scanInit.skyRegionString ) {
1190 LALFree( scanInit.skyRegionString );
1191 }
1192
1193 } else {
1194
1195 /* read skypatch info */
1196 {
1197 FILE *fpsky = NULL;
1198 INT4 r;
1199 REAL8 temp1, temp2, temp3, temp4;
1200
1202
1203 if ( ( fpsky = fopen( skyFileName, "r" ) ) == NULL ) {
1205 }
1206
1207 nSkyPatches = 0;
1208 do {
1209 r = fscanf( fpsky, "%lf%lf%lf%lf\n", &temp1, &temp2, &temp3, &temp4 );
1210 /* make sure the line has the right number of entries or is EOF */
1211 if ( r == 4 ) {
1212 nSkyPatches++;
1213 }
1214 } while ( r != EOF );
1215 rewind( fpsky );
1216
1217 out->numSkyPatches = nSkyPatches;
1218 out->alpha = ( REAL8 * )LALCalloc( nSkyPatches, sizeof( REAL8 ) );
1219 out->delta = ( REAL8 * )LALCalloc( nSkyPatches, sizeof( REAL8 ) );
1220 out->alphaSize = ( REAL8 * )LALCalloc( nSkyPatches, sizeof( REAL8 ) );
1221 out->deltaSize = ( REAL8 * )LALCalloc( nSkyPatches, sizeof( REAL8 ) );
1222
1223 for ( skyCounter = 0; skyCounter < nSkyPatches; skyCounter++ ) {
1224 r = fscanf( fpsky, "%lf%lf%lf%lf\n", out->alpha + skyCounter, out->delta + skyCounter,
1225 out->alphaSize + skyCounter, out->deltaSize + skyCounter );
1226 }
1227
1228 fclose( fpsky );
1229 } /* end skypatchfile reading block */
1230 } /* end setting up of skypatches */
1231
1233
1234 /* normal exit */
1235 RETURN( status );
1236}
1237
1239 REAL8Vector *fdots,
1240 PulsarDopplerParams *thisPoint,
1241 LIGOTimeGPS refTime,
1242 REAL8 f_current,
1243 REAL8 q1_current,
1244 REAL8 q2_current,
1245 REAL8 n_current,
1246 REAL8 fdot_current,
1247 REAL8 fddot_current )
1248{
1249
1250 INT4 i;
1251
1252 INITSTATUS( status );
1254
1255
1256 /**************** Option 1: Searching over spindown parameters ******************/
1257
1258 if ( uvar_QCoeffs ) { /*if searching over q1, q2, n*/
1259
1260
1261 CalculateFdots( status->statusPtr, fdots, f_current, q1_current, q2_current, n_current );
1262 /* initialize Doppler parameters of the potential source */
1263 thisPoint->fkdot[0] = f_current;
1264 for ( i = 1; i < PULSAR_MAX_SPINS; i++ ) {
1265 thisPoint->fkdot[i] = fdots->data[i - 1];
1266 }
1267 thisPoint->refTime = refTime;
1268 } /*endif*/
1269
1270 else { /* if searching through f, fdots instead */
1271
1272
1273 /* initialize Doppler parameters of the potential source */
1274
1275 XLAL_INIT_MEM( thisPoint->fkdot );
1276 thisPoint->fkdot[0] = f_current;
1277 thisPoint->fkdot[1] = fdot_current;
1278 thisPoint->fkdot[2] = fddot_current;
1279 thisPoint->refTime = refTime;
1280 } /*endelse*/
1281
1283
1284 /* normal exit */
1285 RETURN( status );
1286
1287}
1288
1289
1292 SFTListElement *sftHead,
1293 REAL8ListElement *freqHead,
1294 REAL8ListElement *phaseHead,
1295 SkyPosition skypos,
1297 PulsarDopplerParams *thisPoint )
1298{
1299
1300 REAL8 freq1;
1301 REAL8 phase1;
1302 REAL8Vector thisVel, thisPos;
1303 LIGOTimeGPSVector *ts = NULL;
1304 const LALDetector *det;
1305 REAL8 tOffs;
1306 AMCoeffs *AMcoef = NULL;
1307 SFTListElement *sft = NULL;
1308 REAL8ListElement *freqtmp, *phasetmp;
1310 LIGOTimeGPS tgps;
1311
1312 INITSTATUS( status );
1314
1315 freq1 = 0;
1316 phase1 = 0;
1317
1318
1319 sft = sftHead;
1320 freqtmp = freqHead;
1321 phasetmp = phaseHead;
1322 beamtmp = beamHead;
1323
1325
1326
1327 thisVel.length = 3;
1328 thisPos.length = 3;
1329 tOffs = 0.5 / sft->sft.deltaF;
1330
1331 /* get information about all detectors including velocity and
1332 timestamps */
1333 /*only have 1 element in detectorStateSeries and AMCoeffs because
1334 the detector has to be updated for every SFT */
1335 while ( sft ) {
1336
1337 DetectorStateSeries *detState = NULL;
1338
1339 /* get midpoint of sft */
1340 tgps = sft->sft.epoch;
1341 XLALGPSAdd( &tgps, tOffs );
1342
1343 det = XLALGetSiteInfo( sft->sft.name );
1344
1345 ts->data[0] = sft->sft.epoch;
1346 /* note that this function returns the velocity at the
1347 mid-time of the SFTs -- should not make any
1348 difference */
1349
1350 detState = XLALGetDetectorStates( ts, det, edat, tOffs );
1351
1352 detState->detector = *det;
1353
1354 AMcoef = XLALComputeAMCoeffs( detState, skypos );
1355 thisVel.data = detState->data[0].vDetector;
1356 thisPos.data = detState->data[0].rDetector;
1357
1358 LALGetSignalFrequencyInSFT( status->statusPtr, &freq1, &tgps, thisPoint,
1359 &thisVel );
1360
1361 LALGetSignalPhaseInSFT( status->statusPtr, &phase1, &tgps, thisPoint,
1362 &thisPos );
1363
1364 freqtmp->val = freq1;
1365 phasetmp->val = phase1;
1366
1367 /* store a and b in the CrossCorrBeamFn */
1368
1369 beamtmp->beamfn.a = ( AMcoef->a->data[0] );
1370 beamtmp->beamfn.b = ( AMcoef->b->data[0] );
1371
1372 /*
1373 printf("beam A %1.15g\n", beamtmp->beamfn.a);
1374 printf("beam B %1.15g\n", beamtmp->beamfn.b);
1375 printf("vel %1.15g %1.15g %1.15g\n", thisVel.data[0], thisVel.data[1], thisVel.data[2]);
1376 printf("pos %1.15g %1.15g %1.15g\n\n", thisPos.data[0], thisPos.data[1], thisPos.data[2]);
1377 */
1378
1379 /* clean up AMcoefs */
1380 XLALDestroyAMCoeffs( AMcoef );
1382 sft = ( SFTListElement * )sft->nextSFT;
1383 freqtmp = ( REAL8ListElement * )freqtmp->nextVal;
1384 phasetmp = ( REAL8ListElement * )phasetmp->nextVal;
1385 beamtmp = ( CrossCorrBeamFnListElement * )beamtmp->nextBeamfn;
1386 }
1387
1389 sft = NULL;
1390 freqtmp = NULL;
1391 phasetmp = NULL;
1392 beamtmp = NULL;
1393
1395
1396 RETURN( status );
1397
1398}
1399
1400
1402 SFTCatalog *catalog,
1403 SFTVector **sft,
1404 REAL8 fMin,
1405 REAL8 fMax,
1406 INT4 sftindex )
1407{
1408 SFTCatalog *slidingcat;
1410
1411 INITSTATUS( status );
1413
1414
1415 slidingcat = NULL;
1416 *sft = NULL;
1417
1419 /*check that we are loading an sensible frequency range*/
1420 if ( fMin < catalog->data[sftindex].header.f0 || fMax > ( catalog->data[sftindex].header.f0 + catalog->data[sftindex].numBins * catalog->data[sftindex].header.deltaF ) ) {
1422 }
1423
1424 if ( ( slidingcat = LALCalloc( 1, sizeof( SFTCatalog ) ) ) == NULL ) {
1426 }
1427
1428 slidingcat->length = 1;
1429 slidingcat->data = LALRealloc( slidingcat->data, 1 * sizeof( *( slidingcat->data ) ) );
1430 /* memset(&(slidingcat->data[0]), 0, sizeof(ret->data[0]));*/
1431 desc = &( slidingcat->data[0] );
1432 /* desc->locator = LALCalloc(1, sizeof(*(desc->locator)));*/
1433
1434 desc->locator = catalog->data[sftindex].locator;
1435 desc->header = catalog->data[sftindex].header;
1436 desc->comment = catalog->data[sftindex].comment;
1437 desc->numBins = catalog->data[sftindex].numBins;
1438 desc->version = catalog->data[sftindex].version;
1439 desc->crc64 = catalog->data[sftindex].crc64;
1440
1441 *sft = XLALLoadSFTs( slidingcat, fMin, fMax );
1442 /* XLALDestroySFTCatalog (slidingcat );*/
1443 XLALFree( slidingcat->data );
1444 LALFree( slidingcat );
1445
1447
1448 RETURN( status );
1449
1450}
1451
1453 SFTListElement **sftHead,
1454 SFTListElement **sftTail,
1455 SFTtype *sft )
1456{
1457
1458 SFTListElement *sftList;
1459
1460 INITSTATUS( status );
1462
1463 sftList = ( SFTListElement * )LALCalloc( 1, sizeof( SFTListElement ) );
1464
1465 XLAL_CHECK_LAL( status, XLALCopySFT( &( sftList->sft ), sft ) == XLAL_SUCCESS, XLAL_EFUNC );
1466
1467 sftList->nextSFT = NULL;
1468 if ( !( *sftHead ) ) {
1469 *sftHead = sftList;
1470 } else {
1471 ( *sftTail )->nextSFT = ( SFTListElement * )sftList;
1472 }
1473
1474 *sftTail = sftList;
1475
1477
1478 RETURN( status );
1479
1480
1481}
1482
1484 PSDListElement **psdHead,
1485 PSDListElement **psdTail,
1486 INT4 length )
1487{
1488 PSDListElement *psdList;
1489
1490 INITSTATUS( status );
1492
1493 psdList = ( PSDListElement * )LALCalloc( 1, sizeof( PSDListElement ) );
1494 psdList->psd.data = XLALCreateREAL8Sequence( length );
1495 psdList->nextPSD = NULL;
1496 if ( !( *psdHead ) ) {
1497 *psdHead = psdList;
1498 } else {
1499 ( *psdTail )->nextPSD = ( PSDListElement * )psdList;
1500 }
1501 *psdTail = psdList;
1502
1504
1505 RETURN( status );
1506
1507
1508}
1509
1511 REAL8ListElement **head,
1512 REAL8ListElement **tail )
1513{
1514 REAL8ListElement *List;
1515
1516 INITSTATUS( status );
1518
1519
1520 List = ( REAL8ListElement * )LALCalloc( 1, sizeof( REAL8ListElement ) );
1521 List->val = 0;
1522 List->nextVal = NULL;
1523 if ( !( *head ) ) {
1524 *head = List;
1525 } else {
1526 ( *tail )->nextVal = ( REAL8ListElement * )List;
1527 }
1528 *tail = List;
1529
1531
1532 RETURN( status );
1533
1534
1535}
1536
1538 CrossCorrBeamFnListElement **beamHead,
1539 CrossCorrBeamFnListElement **beamTail )
1540{
1541
1543
1544 INITSTATUS( status );
1546
1547
1548 beamList = ( CrossCorrBeamFnListElement * )LALCalloc( 1, sizeof( CrossCorrBeamFnListElement ) );
1549 beamList->beamfn.a = 0;
1550 beamList->beamfn.b = 0;
1551 beamList->nextBeamfn = NULL;
1552 if ( !( *beamHead ) ) {
1553 *beamHead = beamList;
1554 } else {
1555 ( *beamTail )->nextBeamfn = ( CrossCorrBeamFnListElement * )beamList;
1556 }
1557 *beamTail = beamList;
1558
1560
1561 RETURN( status );
1562
1563}
1564
1566 SFTListElement **sftHead )
1567{
1568 SFTListElement *sftList;
1569 SFTtype *tmpSFT;
1570
1571 INITSTATUS( status );
1573
1574
1575 sftList = *sftHead;
1576 *sftHead = ( SFTListElement * )( *sftHead )->nextSFT;
1577
1578 tmpSFT = &( sftList->sft );
1579 XLALDestroySFT( tmpSFT );
1580 sftList = NULL;
1581
1583
1584 RETURN( status );
1585
1586}
1587
1589 PSDListElement **psdHead )
1590{
1591 PSDListElement *psdList;
1592 REAL8FrequencySeries *tmpPSD;
1593
1594 INITSTATUS( status );
1596
1597 psdList = *psdHead;
1598 *psdHead = ( PSDListElement * )( *psdHead )->nextPSD;
1599 tmpPSD = &( psdList->psd );
1601 psdList = NULL;
1602
1604
1605 RETURN( status );
1606
1607}
1608
1610 REAL8ListElement **head )
1611{
1612 REAL8ListElement *List;
1613 REAL8 *tmpVal;
1614
1615 INITSTATUS( status );
1617
1618
1619 List = *head;
1620 *head = ( REAL8ListElement * )( *head )->nextVal;
1621 tmpVal = &( List->val );
1622 LALFree( tmpVal );
1623 List = NULL;
1624
1626
1627 RETURN( status );
1628
1629}
1630
1632 CrossCorrBeamFnListElement **beamHead )
1633{
1635 CrossCorrBeamFn *beamfns;
1636
1637 INITSTATUS( status );
1639
1640 beamList = *beamHead;
1641 *beamHead = ( CrossCorrBeamFnListElement * )( *beamHead )->nextBeamfn;
1642 beamfns = &( beamList->beamfn );
1643 LALFree( beamfns );
1644 beamList = NULL;
1645
1647
1648 RETURN( status );
1649
1650}
1651
1653 REAL8Vector *fdots,
1654 REAL8 f0,
1655 REAL8 q1,
1656 REAL8 q2,
1657 REAL8 n )
1658{
1659 INITSTATUS( status );
1661
1663
1664
1665
1666 q1 = q1 / pow( uvar_fRef, 5 );
1667 q2 = q2 / pow( uvar_fRef, n );
1668
1669 /* hard code each derivative. symbolic differentiation too hard */
1670 fdots->data[0] = -( q1 * pow( f0, 5 ) ) - ( q2 * pow( f0, n ) );
1671
1672 fdots->data[1] = -( 5.0 * q1 * pow( f0, 4 ) * fdots->data[0] )
1673 - ( n * q2 * pow( f0, n - 1 ) * fdots->data[0] );
1674
1675 fdots->data[2] = -q1 * ( 20.0 * CUBE( f0 ) * SQUARE( fdots->data[0] ) + 5.0 * pow( f0, 4 ) * fdots->data[1] )
1676 - q2 * ( ( n - 1 ) * n * pow( f0, n - 2 ) * SQUARE( fdots->data[0] ) + n * pow( f0, n - 1 ) * fdots->data[1] );
1677
1678 fdots->data[3] = -q1 * ( 60.0 * SQUARE( f0 ) * CUBE( fdots->data[0] ) + 60.0 * CUBE( f0 ) * fdots->data[0] * fdots->data[1]
1679 + 5.0 * pow( f0, 4 ) * fdots->data[2] )
1680 - q2 * ( ( n - 2 ) * ( n - 1 ) * n * pow( f0, n - 3 ) * CUBE( fdots->data[0] )
1681 + 3 * n * ( n - 1 ) * pow( f0, n - 2 ) * fdots->data[0] * fdots->data[1] + n * pow( f0, n - 1 ) * fdots->data[2] );
1682
1683 fdots->data[4] = -q1 * ( 120.0 * f0 * pow( fdots->data[0], 4 ) + 360.0 * SQUARE( f0 ) * SQUARE( fdots->data[0] ) * fdots->data[1]
1684 + 60.0 * CUBE( f0 ) * SQUARE( fdots->data[1] ) + 80.0 * CUBE( f0 ) * fdots->data[0] * fdots->data[2]
1685 + 5.0 * SQUARE( f0 ) * SQUARE( f0 ) * fdots->data[3] )
1686 - q2 * ( ( n - 3 ) * ( n - 2 ) * ( n - 1 ) * n * pow( f0, n - 4 ) * pow( fdots->data[0], 4 )
1687 + 6.0 * ( n - 2 ) * ( n - 1 ) * n * pow( f0, n - 3 ) * SQUARE( fdots->data[0] ) * fdots->data[1]
1688 + 3.0 * ( n - 1 ) * n * pow( f0, n - 2 ) * SQUARE( fdots->data[1] )
1689 + 4.0 * ( n - 1 ) * n * pow( f0, n - 2 ) * fdots->data[0] * fdots->data[2] + n * pow( f0, n - 1 ) * fdots->data[3] );
1690
1691 fdots->data[5] = -q1 * ( 120.0 * pow( fdots->data[0], 5 ) + 1200.0 * f0 * CUBE( fdots->data[0] ) * fdots->data[1]
1692 + 900.0 * SQUARE( f0 ) * fdots->data[0] * SQUARE( fdots->data[1] ) + 600.0 * SQUARE( f0 ) * SQUARE( fdots->data[0] ) * fdots->data[2]
1693 + 200.0 * CUBE( f0 ) * fdots->data[1] * fdots->data[2] + 100.0 * CUBE( f0 ) * fdots->data[0] * fdots->data[3]
1694 + 5.0 * pow( f0, 4 ) * fdots->data[4] )
1695 - q2 * ( ( n - 4 ) * ( n - 3 ) * ( n - 2 ) * ( n - 1 ) * n * pow( f0, n - 5 ) * pow( fdots->data[0], 5 )
1696 + 10.0 * ( n - 3 ) * ( n - 2 ) * ( n - 1 ) * n * pow( f0, n - 4 ) * CUBE( fdots->data[0] ) * fdots->data[1]
1697 + 15.0 * ( n - 2 ) * ( n - 1 ) * n * pow( f0, n - 3 ) * fdots->data[0] * SQUARE( fdots->data[1] )
1698 + 10.0 * ( n - 2 ) * ( n - 1 ) * n * pow( f0, n - 3 ) * SQUARE( fdots->data[1] ) * fdots->data[2]
1699 + 10.0 * ( n - 1 ) * n * pow( f0, n - 2 ) * fdots->data[1] * fdots->data[2]
1700 + 5.0 * ( n - 1 ) * n * pow( f0, n - 2 ) * fdots->data[0] * fdots->data[3]
1701 + n * pow( f0, n - 1 ) * fdots->data[4] );
1702
1704
1705 RETURN( status );
1706}
1707
1709{
1710 INITSTATUS( status );
1712
1713
1714
1715 uvar_maxlag = 0;
1716
1723
1724 uvar_detChoice = 2;
1725 uvar_f0 = F0;
1726 uvar_fBand = FBAND;
1728 uvar_startTime = 0.0;
1730 uvar_fdot = 0.0;
1731 uvar_fdotBand = 0.0;
1732 uvar_fdotResolution = 0.0;
1733 uvar_fddot = 0.0;
1734 uvar_fddotBand = 0.0;
1736 uvar_dAlpha = 0.0;
1737 uvar_dDelta = 0.0;
1738 uvar_psi = 0.0;
1739 uvar_refTime = 0.0;
1740 uvar_cosi = 0.0;
1741 uvar_q1 = 1e-24;
1742 uvar_q1Band = 0.0;
1743 uvar_q1Resolution = 1e-25;
1744 uvar_q2 = 1e-20;
1745 uvar_q2Band = 0.0;
1746 uvar_q2Resolution = 1e-21;
1750 uvar_fRef = 1.0;
1751
1752 uvar_ephemEarth = XLALStringDuplicate( "earth00-40-DE405.dat.gz" );
1753 uvar_ephemSun = XLALStringDuplicate( "sun00-40-DE405.dat.gz" );
1754
1756 strcpy( uvar_dirnameOut, DIROUT );
1757
1759 strcpy( uvar_filenameOut, FILEOUT );
1760
1761 uvar_debugOut = ( CHAR * )LALCalloc( MAXFILENAMELENGTH, sizeof( CHAR ) );
1762 strcpy( uvar_debugOut, DEBUGOUT );
1763
1764 uvar_skyfile = ( CHAR * )LALCalloc( MAXFILENAMELENGTH, sizeof( CHAR ) );
1765 strcpy( uvar_skyfile, SKYFILE );
1766
1767 /* register user input variables */
1768 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_averagePsi, "averagePsi", BOOLEAN, 0, OPTIONAL, "Use average over psi" ) == XLAL_SUCCESS, XLAL_EFUNC );
1769 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_averageIota, "averageIota", BOOLEAN, 0, OPTIONAL, "Use average over iota" ) == XLAL_SUCCESS, XLAL_EFUNC );
1770 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_autoCorrelate, "autoCorrelate", BOOLEAN, 0, OPTIONAL, "Include autocorrelations" ) == XLAL_SUCCESS, XLAL_EFUNC );
1771 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_timingOn, "timingOn", BOOLEAN, 0, OPTIONAL, "Print code timing information" ) == XLAL_SUCCESS, XLAL_EFUNC );
1772 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_f0, "f0", REAL8, 'f', OPTIONAL, "Start search frequency" ) == XLAL_SUCCESS, XLAL_EFUNC );
1773 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_fdot, "fdot", REAL8, 0, OPTIONAL, "Start search frequency derivative" ) == XLAL_SUCCESS, XLAL_EFUNC );
1774 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_fBand, "fBand", REAL8, 'b', OPTIONAL, "Search frequency band" ) == XLAL_SUCCESS, XLAL_EFUNC );
1775 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_fResolution, "fRes", REAL8, 0, OPTIONAL, "Search frequency resolution. Default: 1/T" ) == XLAL_SUCCESS, XLAL_EFUNC );
1776 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_fdotBand, "fdotBand", REAL8, 0, OPTIONAL, "Search frequency derivative band" ) == XLAL_SUCCESS, XLAL_EFUNC );
1777 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_fdotResolution, "fdotRes", REAL8, 'r', OPTIONAL, "Search frequency derivative resolution. Default: 1/T^2" ) == XLAL_SUCCESS, XLAL_EFUNC );
1778 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_fRef, "fRef", REAL8, 0, OPTIONAL, "Reference frequency" ) == XLAL_SUCCESS, XLAL_EFUNC );
1779 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_fddot, "fddot", REAL8, 0, OPTIONAL, "Start frequency double derivative" ) == XLAL_SUCCESS, XLAL_EFUNC );
1780 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_fddotBand, "fddotBand", REAL8, 0, OPTIONAL, "Search frequency double derivative band" ) == XLAL_SUCCESS, XLAL_EFUNC );
1781 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_fddotResolution, "fddotRes", REAL8, 0, OPTIONAL, "Search frequency double derivative resolution. Default: 1/T^3" ) == XLAL_SUCCESS, XLAL_EFUNC );
1782 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_startTime, "startTime", REAL8, 0, OPTIONAL, "GPS start time of observation (SFT timestamps must be >= this)" ) == XLAL_SUCCESS, XLAL_EFUNC );
1783 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_endTime, "endTime", REAL8, 0, OPTIONAL, "GPS end time of observation (SFT timestamps must be < this)" ) == XLAL_SUCCESS, XLAL_EFUNC );
1784 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_skyRegion, "skyRegion", STRING, 0, OPTIONAL, "sky-region polygon (or 'allsky')" ) == XLAL_SUCCESS, XLAL_EFUNC );
1785 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_dAlpha, "dAlpha", REAL8, 0, OPTIONAL, "Sky resolution (flat/isotropic) (rad)" ) == XLAL_SUCCESS, XLAL_EFUNC );
1786 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_dDelta, "dDelta", REAL8, 0, OPTIONAL, "Sky resolution (flat/isotropic) (rad)" ) == XLAL_SUCCESS, XLAL_EFUNC );
1787 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_psi, "psi", REAL8, 0, OPTIONAL, "Polarisation angle (rad)" ) == XLAL_SUCCESS, XLAL_EFUNC );
1788 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_skyfile, "skyfile", STRING, 0, OPTIONAL, "Alternative: input skypatch file" ) == XLAL_SUCCESS, XLAL_EFUNC );
1789
1790 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_ephemEarth, "ephemEarth", STRING, 0, OPTIONAL, "Earth ephemeris file to use" ) == XLAL_SUCCESS, XLAL_EFUNC );
1791 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_ephemSun, "ephemSun", STRING, 0, OPTIONAL, "Sun ephemeris file to use" ) == XLAL_SUCCESS, XLAL_EFUNC );
1792
1793 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_sftDir, "sftDir", STRING, 'D', REQUIRED, "SFT filename pattern. Possibilities are:\n"
1794 " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'" ) == XLAL_SUCCESS, XLAL_EFUNC );
1795 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_maxlag, "maxlag", REAL8, 0, OPTIONAL, "Maximum time lag for correlating sfts" ) == XLAL_SUCCESS, XLAL_EFUNC );
1796 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_dirnameOut, "dirnameOut", STRING, 'o', OPTIONAL, "Output directory" ) == XLAL_SUCCESS, XLAL_EFUNC );
1797 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_filenameOut, "filenameOut", STRING, 0, OPTIONAL, "Output filename" ) == XLAL_SUCCESS, XLAL_EFUNC );
1798 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_debugOut, "debugOut", STRING, 0, OPTIONAL, "Debugging output filename" ) == XLAL_SUCCESS, XLAL_EFUNC );
1799 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_blocksRngMed, "blocksRngMed", INT4, 0, OPTIONAL, "Running Median block size" ) == XLAL_SUCCESS, XLAL_EFUNC );
1800 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_detChoice, "detChoice", INT4, 0, OPTIONAL, "0: Correlate SFTs from same IFOs only; 1: different IFOs only; 2: all IFOs" ) == XLAL_SUCCESS, XLAL_EFUNC );
1801 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_refTime, "refTime", REAL8, 0, OPTIONAL, "Pulsar reference time (SSB)" ) == XLAL_SUCCESS, XLAL_EFUNC );
1802 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_cosi, "cosi", REAL8, 0, OPTIONAL, "cos(iota) inclination angle" ) == XLAL_SUCCESS, XLAL_EFUNC );
1803 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_q1, "q1", REAL8, 0, OPTIONAL, "Starting Q1 value" ) == XLAL_SUCCESS, XLAL_EFUNC );
1804 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_q1Band, "q1Band", REAL8, 0, OPTIONAL, "Q1 search band" ) == XLAL_SUCCESS, XLAL_EFUNC );
1805 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_q1Resolution, "q1Res", REAL8, 0, OPTIONAL, "Pulsar ellipticity search resolution" ) == XLAL_SUCCESS, XLAL_EFUNC );
1806 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_q2, "q2", REAL8, 0, OPTIONAL, "Starting Q2 value" ) == XLAL_SUCCESS, XLAL_EFUNC );
1807 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_q2Band, "q2Band", REAL8, 0, OPTIONAL, "Q2 search band" ) == XLAL_SUCCESS, XLAL_EFUNC );
1808 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_q2Resolution, "q2Res", REAL8, 0, OPTIONAL, "Q2 search resolution" ) == XLAL_SUCCESS, XLAL_EFUNC );
1809 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_brakingindex, "braking", REAL8, 0, OPTIONAL, "Pulsar electromagnetic braking index" ) == XLAL_SUCCESS, XLAL_EFUNC );
1810 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_brakingindexBand, "brakingBand", REAL8, 0, OPTIONAL, "Pulsar electromagnetic braking index search band" ) == XLAL_SUCCESS, XLAL_EFUNC );
1811 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_brakingindexResolution, "brakingRes", REAL8, 0, OPTIONAL, "Pulsar electromagnetic braking index search resolution" ) == XLAL_SUCCESS, XLAL_EFUNC );
1812 XLAL_CHECK_LAL( status, XLALRegisterNamedUvar( &uvar_QCoeffs, "useQCoeffs", BOOLEAN, 0, OPTIONAL, "Search over pulsar spindown parameters instead of frequency derivatives" ) == XLAL_SUCCESS, XLAL_EFUNC );
1813
1814
1816 RETURN( status );
1817}
void InitDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan, const DopplerSkyScanInit *init)
Definition: DopplerScan.c:273
void FreeDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan)
Definition: DopplerScan.c:324
int XLALNextDopplerSkyPos(PulsarDopplerParams *pos, DopplerSkyScanState *skyScan)
NextDopplerSkyPos(): step through sky-grid return 0 = OK, -1 = ERROR.
Definition: DopplerScan.c:127
@ STATE_FINISHED
all templates have been read
Definition: DopplerScan.h:84
@ GRID_ISOTROPIC
approximately isotropic sky-grid
Definition: DopplerScan.h:92
#define LAL_INT4_MAX
void REPORTSTATUS(LALStatus *status)
lal_errhandler_t lal_errhandler
int LAL_ERR_EXIT(LALStatus *stat, const char *func, const char *file, const int line, volatile const char *id)
#define LAL_CALL(function, statusptr)
int k
void LALCheckMemoryLeaks(void)
#define LALRealloc(p, n)
#define LALCalloc(m, n)
#define LALFree(p)
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
#define ABORT(statusptr, code, mesg)
#define XLAL_CHECK_LAL(sp, assertion,...)
#define TRY(func, statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define ASSERT(assertion, statusptr, code, mesg)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define STRING(a)
#define fscanf
#define fprintf
double e
void XLALDestroyDetectorStateSeries(DetectorStateSeries *detStates)
Get rid of a DetectorStateSeries.
DetectorStateSeries * XLALGetDetectorStates(const LIGOTimeGPSVector *timestamps, const LALDetector *detector, const EphemerisData *edat, REAL8 tOffset)
Get the 'detector state' (ie detector-tensor, position, velocity, etc) for the given vector of timest...
int XLALExtrapolatePulsarSpinRange(PulsarSpinRange *range1, const PulsarSpinRange *range0, const REAL8 dtau)
General pulsar-spin extrapolation function: given a "spin-range" (ie spins + spin-bands) range0 at ti...
void XLALDestroyREAL8FrequencySeries(REAL8FrequencySeries *series)
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
void XLALDestroyAMCoeffs(AMCoeffs *amcoef)
Destroy a AMCoeffs structure.
Definition: LALComputeAM.c:497
AMCoeffs * XLALComputeAMCoeffs(const DetectorStateSeries *DetectorStates, SkyPosition skypos)
Compute the 'amplitude coefficients' , as defined in for a series of timestamps.
Definition: LALComputeAM.c:297
#define LAL_PI
unsigned char BOOLEAN
#define XLAL_INIT_MEM(x)
double REAL8
#define XLAL_INIT_DECL(var,...)
int64_t INT8
char CHAR
uint32_t UINT4
int32_t INT4
#define lalDebugLevel
void XLALFree(void *p)
#define LAL_INT8_FORMAT
char char * XLALStringDuplicate(const char *s)
#define VTOT
Total detector velocity/c TO BE CHANGED DEPENDING ON DETECTOR.
Definition: LUT.h:207
int XLALNormalizeSFT(REAL8FrequencySeries *rngmed, SFTtype *sft, UINT4 blockSize, const REAL8 assumeSqrtS)
Normalize an sft based on RngMed estimated PSD, and returns running-median.
@ LAL_PMETRIC_NONE
Definition: PtoleMetric.h:96
void LALCalculateEstimators(LALStatus *status, REAL8 *aplussq1, REAL8 *aplussq2, REAL8 *acrossq1, REAL8 *acrossq2, COMPLEX16Vector *yalpha, COMPLEX16Vector *gplus, COMPLEX16Vector *gcross, REAL8Vector *sigmaAlphasq)
void LALCalculateAveUalpha(LALStatus *status, COMPLEX16 *out, REAL8 phiI, REAL8 phiJ, REAL8 freqI, REAL8 freqJ, REAL8 deltaF, CrossCorrBeamFn beamfnsI, CrossCorrBeamFn beamfnsJ, REAL8 sigmasq)
Calculate pair weights (U_alpha) for an average over Psi and cos(iota)
void LALCalculateSigmaAlphaSq(LALStatus *status, REAL8 *out, UINT4 bin1, UINT4 bin2, REAL8FrequencySeries *psd1, REAL8FrequencySeries *psd2)
void LALGetSignalFrequencyInSFT(LALStatus *status, REAL8 *out, LIGOTimeGPS *epoch, PulsarDopplerParams *dopp, REAL8Vector *vel_c)
Calculate the frequency of the SFT at a given epoch.
void LALCorrelateSingleSFTPair(LALStatus *status, COMPLEX16 *out, COMPLEX8FrequencySeries *sft1, COMPLEX8FrequencySeries *sft2, REAL8FrequencySeries *psd1, REAL8FrequencySeries *psd2, UINT4 bin1, UINT4 bin2)
Correlate a single pair of SFT at a parameter space point.
void LALCalculateUalpha(LALStatus *status, COMPLEX16 *out, CrossCorrAmps amplitudes, REAL8 phiI, REAL8 phiJ, REAL8 freqI, REAL8 freqJ, REAL8 deltaF, CrossCorrBeamFn beamfnsI, CrossCorrBeamFn beamfnsJ, REAL8 sigmasq, REAL8 *psi, COMPLEX16 *gplus, COMPLEX16 *gcross)
Calculate pair weights (U_alpha) for the general case.
void LALGetSignalPhaseInSFT(LALStatus *status, REAL8 *out, LIGOTimeGPS *epoch, PulsarDopplerParams *dopp, REAL8Vector *r_c)
Get signal phase at a given epoch.
void LALCalculateCrossCorrPower(LALStatus *status, REAL8 *out, COMPLEX16Vector *yalpha, COMPLEX16Vector *ualpha)
void LALNormaliseCrossCorrPower(LALStatus *status, REAL8 *out, COMPLEX16Vector *ualpha, REAL8Vector *sigmaAlphasq)
DetChoice
@ ALL
#define PULSAR_MAX_SPINS
maximal number of spin-parameters (Freq + spindowns) we can handle
static const INT4 r
const LALDetector * XLALGetSiteInfo(const CHAR *name)
Find the site geometry-information 'LALDetector' for given a detector name (or prefix).
Definition: SFTnaming.c:218
void XLALDestroySFTVector(SFTVector *vect)
XLAL interface to destroy an SFTVector.
Definition: SFTtypes.c:300
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
Definition: SFTcatalog.c:329
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:47
int XLALCopySFT(SFTtype *dest, const SFTtype *src)
Copy an entire SFT-type into another.
Definition: SFTtypes.c:202
SFTVector * XLALLoadSFTs(const SFTCatalog *catalog, REAL8 fMin, REAL8 fMax)
Load the given frequency-band [fMin, fMax) (half-open) from the SFT-files listed in the SFT-'catalogu...
Definition: SFTfileIO.c:87
SFTCatalog * XLALSFTdataFind(const CHAR *file_pattern, const SFTConstraints *constraints)
Find the list of SFTs matching the file_pattern and satisfying the given constraints,...
Definition: SFTcatalog.c:71
void XLALDestroySFT(SFTtype *sft)
Destructor for one SFT.
Definition: SFTtypes.c:176
void XLALDestroyTimestampVector(LIGOTimeGPSVector *vect)
De-allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:69
REAL8Sequence * XLALCreateREAL8Sequence(size_t length)
COORDINATESYSTEM_EQUATORIAL
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
void XLALDestroyUserVars(void)
#define XLALRegisterNamedUvar(cvar, name, type, option, category,...)
int XLALUserVarWasSet(const void *cvar)
REAL8Vector * XLALResizeREAL8Vector(REAL8Vector *vector, UINT4 length)
COMPLEX16Vector * XLALResizeCOMPLEX16Vector(COMPLEX16Vector *vector, UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
void XLALDestroyCOMPLEX16Vector(COMPLEX16Vector *vector)
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_MAIN(assertion,...)
XLAL_SUCCESS
XLAL_EFUNC
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
LIGOTimeGPS * XLALGPSSet(LIGOTimeGPS *epoch, INT4 gpssec, INT8 gpsnan)
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
REAL8 XLALGPSDiff(const LIGOTimeGPS *t1, const LIGOTimeGPS *t0)
float data[BLOCKSIZE]
Definition: hwinject.c:360
head
n
out
desc
ts
INT4 uvar_detChoice
#define N_SPINDOWN_DERIVS
REAL8 uvar_cosi
REAL8 uvar_maxlag
CHAR * uvar_skyRegion
REAL8 uvar_fResolution
#define SKYFILE
REAL8 uvar_fdotResolution
int main(int argc, char *argv[])
void initUserVars(LALStatus *status)
void DeleteSFTHead(LALStatus *status, SFTListElement **sftHead)
REAL8 uvar_q1Resolution
REAL8 uvar_dAlpha
REAL8 uvar_q2Resolution
void AddBeamFntoList(LALStatus *status, CrossCorrBeamFnListElement **beamHead, CrossCorrBeamFnListElement **beamTail)
REAL8 uvar_refTime
BOOLEAN uvar_averagePsi
#define F0
REAL8 uvar_fdot
void CopySFTFromCatalog(LALStatus *status, SFTCatalog *catalog, SFTVector **sft, REAL8 fMin, REAL8 fMax, INT4 sftindex)
BOOLEAN uvar_QCoeffs
void AddPSDtoList(LALStatus *status, PSDListElement **psdHead, PSDListElement **psdTail, INT4 length)
#define DEBUGOUT
CHAR * uvar_sftDir
REAL8 uvar_brakingindexResolution
REAL8 uvar_q2Band
REAL8 uvar_brakingindex
REAL8 uvar_dDelta
void DeletePSDHead(LALStatus *status, PSDListElement **psdHead)
REAL8 uvar_fddotBand
REAL8 uvar_q2
REAL8 uvar_brakingindexBand
void CalculateFdots(LALStatus *status, REAL8Vector *fdots, REAL8 f0, REAL8 q1, REAL8 q2, REAL8 n)
void AddSFTtoList(LALStatus *status, SFTListElement **sftHead, SFTListElement **sftTail, SFTtype *sft)
REAL8 uvar_fddotResolution
#define DIROUT
void InitDoppParams(LALStatus *status, REAL8Vector *fdots, PulsarDopplerParams *thisPoint, LIGOTimeGPS refTime, REAL8 f_current, REAL8 q1_current, REAL8 q2_current, REAL8 n_current, REAL8 fdot_current, REAL8 fddot_current)
void DeleteREAL8Head(LALStatus *status, REAL8ListElement **head)
CHAR * uvar_ephemEarth
Earth ephemeris file to use.
REAL8 uvar_startTime
REAL8 uvar_fRef
BOOLEAN uvar_timingOn
REAL8 uvar_f0
BOOLEAN uvar_autoCorrelate
void DeleteBeamFnHead(LALStatus *status, CrossCorrBeamFnListElement **beamHead)
#define SQUARE(x)
#define TRUE
#define FALSE
CHAR * uvar_filenameOut
CHAR * uvar_skyfile
#define CUBE(x)
REAL8 uvar_fdotBand
CHAR * uvar_dirnameOut
#define FILEOUT
REAL8 uvar_q1Band
BOOLEAN uvar_averageIota
#define BLOCKSRNGMED
REAL8 uvar_fBand
void GetBeamInfo(LALStatus *status, CrossCorrBeamFnListElement *beamHead, SFTListElement *sftHead, REAL8ListElement *freqHead, REAL8ListElement *phaseHead, SkyPosition skypos, EphemerisData *edat, PulsarDopplerParams *thisPoint)
CHAR * uvar_debugOut
#define FBAND
REAL8 uvar_q1
INT4 uvar_blocksRngMed
void AddREAL8toList(LALStatus *status, REAL8ListElement **head, REAL8ListElement **tail)
#define MAXFILENAMELENGTH
REAL8 uvar_fddot
REAL8 uvar_endTime
REAL8 uvar_psi
CHAR * uvar_ephemSun
Sun ephemeris file to use.
void SetUpRadiometerSkyPatches(LALStatus *status, SkyPatchesInfo *out, CHAR *skyFileName, CHAR *skyRegion, REAL8 dAlpha, REAL8 dDelta)
Set up location of skypatch centers and sizes If user specified skyRegion then use DopplerScan functi...
Header for CW cross-correlation search.
#define PULSAR_CROSSCORR_ENULL
#define PULSAR_CROSSCORR_MSGEBAD
#define PULSAR_CROSSCORR_EBAD
#define PULSAR_CROSSCORR_MSGEMEM
#define PULSAR_CROSSCORR_EMEM
#define PULSAR_CROSSCORR_MSGENULL
#define PULSAR_CROSSCORR_MSGEVAL
#define PULSAR_CROSSCORR_EFILE
#define PULSAR_CROSSCORR_EVAL
#define PULSAR_CROSSCORR_MSGEFILE
This structure contains the per-SFT (weighted) antenna-pattern functions , with the SFT-index,...
Definition: LALComputeAM.h:63
REAL4Vector * b
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:65
REAL4Vector * a
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:64
COMPLEX16 * data
CHAR name[LALNameLength]
A vector of COMPLEX8FrequencySeries.
COMPLEX8FrequencySeries * data
Pointer to the data array.
struct tagCrossCorrBeamFnListElement * nextBeamfn
REAL8 vDetector[3]
Cart.
REAL8 rDetector[3]
Cartesian coords of detector position in ICRS J2000.
Timeseries of DetectorState's, representing the detector-info at different timestamps.
DetectorState * data
array of DetectorState entries
LALDetector detector
detector-info corresponding to this timeseries
initialization-structure passed to InitDopplerSkyScan()
Definition: DopplerScan.h:134
this structure reflects the current state of a DopplerSkyScan
Definition: DopplerScan.h:152
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
A vector of 'timestamps' of type LIGOTimeGPS.
Definition: SFTfileIO.h:188
REAL8FrequencySeries psd
struct tagPSDListElement * nextPSD
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
PulsarSpins fkdot
Intrinsic spins: [Freq, f1dot, f2dot, ... ] where fkdot = d^kFreq/dt^k.
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 .
PulsarSpins fkdot
Vector of spin-values .
LIGOTimeGPS refTime
SSB reference GPS-time at which spin-range is defined.
PulsarSpins fkdotBand
Vector of spin-bands , MUST be same length as fkdot.
REAL4 * data
REAL8Sequence * data
struct tagREAL8ListElement * nextVal
REAL8 * data
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238
SFTDescriptor * data
array of data-entries describing matched SFTs
Definition: SFTfileIO.h:243
UINT4 length
number of SFTs in catalog
Definition: SFTfileIO.h:242
'Constraints' for SFT-matching: which detector, within which time-stretch and which timestamps exactl...
Definition: SFTfileIO.h:212
CHAR * detector
2-char channel-prefix describing the detector (eg 'H1', 'H2', 'L1', 'G1' etc)
Definition: SFTfileIO.h:213
LIGOTimeGPSVector * timestamps
list of timestamps
Definition: SFTfileIO.h:216
LIGOTimeGPS * maxStartTime
only include SFTs whose epoch is < maxStartTime
Definition: SFTfileIO.h:215
LIGOTimeGPS * minStartTime
only include SFTs whose epoch is >= minStartTime
Definition: SFTfileIO.h:214
A 'descriptor' of an SFT: basically containing the header-info plus an opaque description of where ex...
Definition: SFTfileIO.h:226
SFTtype header
SFT-header info.
Definition: SFTfileIO.h:228
UINT4 version
SFT-specification version.
Definition: SFTfileIO.h:233
struct tagSFTLocator * locator
internal description of where to find this SFT [opaque!]
Definition: SFTfileIO.h:227
UINT4 numBins
number of frequency-bins in this SFT
Definition: SFTfileIO.h:232
UINT8 crc64
crc64 checksum
Definition: SFTfileIO.h:234
CHAR * comment
comment-entry in SFT-header
Definition: SFTfileIO.h:231
struct tagSFTListElement * nextSFT
struct holding info about skypoints
REAL8 longitude
REAL8 latitude
CoordinateSystem system