forked from scottransom/presto
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpresto.h
1550 lines (1373 loc) · 82.1 KB
/
presto.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#include <sys/times.h>
#include "chkio.h"
#include <math.h>
#include <string.h>
#include <time.h>
#include "clk_tck.h"
#include "ransomfft.h"
#include "meminfo.h"
#include "vectors.h"
#include "makeinf.h"
#include "orbint.h"
#include "misc_utils.h"
#include "cdflib.h"
#include "database.h"
#include "makedata.h"
#ifndef SQRT2
#define SQRT2 1.4142135623730950488016887242096980785696718753769
#endif
#ifndef DBLCORRECT
#define DBLCORRECT 1e-14
#endif
#ifndef PI
#define PI 3.1415926535897932384626433832795028841971693993751
#endif
#ifndef TWOPI
#define TWOPI 6.2831853071795864769252867665590057683943387987502
#endif
#ifndef DEGTORAD
#define DEGTORAD 0.017453292519943295769236907684886127134428718885417
#endif
#ifndef RADTODEG
#define RADTODEG 57.29577951308232087679815481410517033240547246656
#endif
#ifndef PIBYTWO
#define PIBYTWO 1.5707963267948966192313216916397514420985846996876
#endif
#ifndef SOL
#define SOL 299792458.0
#endif
#ifndef SECPERJULYR
#define SECPERJULYR 31557600.0
#endif
#ifndef SECPERDAY
#define SECPERDAY 86400.0
#endif
#ifndef ARCSEC2RAD
#define ARCSEC2RAD 4.8481368110953599358991410235794797595635330237270e-6
#endif
#ifndef SEC2RAD
#define SEC2RAD 7.2722052166430399038487115353692196393452995355905e-5
#endif
#ifndef __GNUC__
#define __inline__
#endif
/* Maximum number of input files to try and patch together */
#define MAXPATCHFILES 100
/* Blocksize to use when reading datafiles or subbands */
#define SUBSBLOCKLEN 1024
/* various function-like macros */
#ifndef SWAP
/* Swaps two variables of undetermined type */
#define SWAP(a,b) tempzz=(a);(a)=(b);(b)=tempzz;
#endif
#ifndef POWER
/* Returns unnormalized Fourier power */
/* Requires the following variables in calling function */
/* double powargr, powargi; */
#define POWER(r,i) (powargr=(r),powargi=(i),\
powargr*powargr+powargi*powargi)
#endif
#ifndef PHASE
/* Returns Fourier phase (degrees) */
/* Requires the following variables in calling function */
/* double phsargr, phsargi, phstmp; */
#define PHASE(r,i) (phsargr=(r),phsargi=(i),\
((phstmp=RADTODEG*atan2(phsargi,phsargr)) > 0.0) ? \
phstmp : phstmp+360.0)
#endif
#ifndef RADIAN_PHASE
/* Returns Fourier phase (radians) */
/* Requires the following variables in calling function */
/* double radargr, radargi, radtmp; */
#define RADIAN_PHASE(r,i) (radargr=(r),radargi=(i),\
((radtmp=atan2(radargi,radargr)) > 0.0) ? \
radtmp : radtmp+TWOPI)
#endif
#define GET_BIT(c, n) (*(c+(n>>3)) >> (7-(n&7)) & 1)
#define SET_BIT(c, n) (*(c+(n>>3)) |= 1 << (7-(n&7)))
#define UNSET_BIT(c, n) (*(c+(n>>3)) &= ~(1 << (7-(n&7))))
/* Number of bins (total) to average for local power: */
/* Must be an even number (1/2 on each side). */
#define NUMLOCPOWAVG 20
/* Number of bins next to freq in question to */
/* ignore (on each side) when determining local power. */
#define DELTAAVGBINS 5
/* Number of bins on each side of the central frequency */
/* to sum for Fourier interpolation (low accuracy) */
#define NUMFINTBINS 16
/* Used for raw-data handling */
typedef enum {
IF0, IF1, SUMIFS
} IFs;
/* Constants used in the interpolation routines */
typedef enum {
LOWACC, HIGHACC
} presto_interp_acc;
/* Constants used in the binary search routines */
typedef enum {
INTERBIN, INTERPOLATE
} presto_interptype;
typedef enum {
NO_CHECK_ALIASED, CHECK_ALIASED
} presto_checkaliased;
/* Constants used in the correlation/convolution routines */
typedef enum {
CONV, CORR, INPLACE_CONV, INPLACE_CORR
} presto_optype;
typedef enum {
FFTDK, FFTD, FFTK, NOFFTS
} presto_ffts;
typedef enum {
RAW, PREPPED, FFT, SAME
} presto_datainf;
typedef struct FFTCAND {
int nsum; /* Number of harmonics summed */
double p; /* Summed normalized power */
double r; /* Frequency of lowest harmonic */
double sig; /* Significance of the candidate taking N into account */
} fftcand;
typedef struct POSITION {
float pow; /* Power normalized with local power */
double p1; /* r if rzw search, r_startfft if bin search */
double p2; /* z if rzw search, r_freqmod if bin search */
double p3; /* w if rzw search, numfftbins if bin search */
} position;
typedef struct POWINDEX {
float pow; /* Power of Fourier bin */
long ind; /* Array index of that Fourier bin */
} powindex;
typedef struct RDERIVS {
double pow; /* Power normalized with local power */
double phs; /* Signal phase */
double dpow; /* 1st deriv of power wrt fourier freq */
double dphs; /* 1st deriv of phase wrt fourier freq */
double d2pow; /* 2nd deriv of power wrt fourier freq */
double d2phs; /* 2nd deriv of power wrt fourier freq */
double locpow; /* Local mean power level */
} rderivs;
typedef struct FOURIERPROPS {
double r; /* Fourier freq at start of run */
float rerr; /* Error in fourier freq */
double z; /* Fourier freq deriv at start of run = T^2*fdot */
float zerr; /* Error in fourier freq deriv */
double w; /* Fourier 2nd freq deriv = T^3*fdotdot */
float werr; /* Error in 2nd fourier freq deriv */
float pow; /* Power normalized with local power */
float powerr; /* Error in power */
float sig; /* Approx. Sigma level */
float rawpow; /* Raw signal power (unnormalized) */
float phs; /* Signal phase */
float phserr; /* Error in phase */
float cen; /* Centroid of signal pulsations over run */
float cenerr; /* Error in centroid */
float pur; /* Purity: rms duration of signal over run */
float purerr; /* Error in purity */
float locpow; /* Local mean power level */
} fourierprops;
typedef struct BINARYPROPS {
double ppsr; /* Pulsar period (s) */
double fpsr; /* Pulsar freq (hz) */
double rpsr; /* Pulsar Fourier freq (bins) */
double pbin; /* Binary period (s) */
double rbin; /* Binary Fourier freq (bins) */
double z; /* Amplitude of frequency modulation */
double asinic; /* Projected semi-major axis of orbit (lt-sec) */
double rdetect; /* Mini-FFT candidate detected freq (bins) */
long nfftbins; /* Number of FFT bins in mini-fft */
long lowbin; /* Start bin of mini-FFT in original FFT */
float ppsrerr; /* Error in pulsar period (s) */
float fpsrerr; /* Error in pulsar freq (hz) */
float rpsrerr; /* Error in pulsar Fourier freq (bins) */
float pbinerr; /* Error in binary period (s) */
float rbinerr; /* Error in pulsar Fourier freq (bins) */
float zerr; /* Error in freq modulation ampltude */
float asinicerr; /* Error in a*sin(i)/c in (lt-sec) */
float rdetecterr; /* Error in Mini-FFT detected freq (bins) */
float sig; /* Approx. Sigma level */
float phs; /* Signal phase */
float phserr; /* Error in phase */
float cen; /* Centroid of signal pulsations over run */
float cenerr; /* Error in centroid */
float pur; /* Purity: rms duration of signal over run */
float purerr; /* Error in purity */
float pow; /* Power normalized with Nph equivalent */
float powerr; /* Error in signal power */
} binaryprops;
typedef struct RAWBINCAND {
double full_N; /* Number of points in original time series */
double full_T; /* Length (s) of original time series */
double full_lo_r; /* Lowest Fourier bin that was miniFFTd */
double mini_N; /* Number of points in short (mini) FFT */
double mini_r; /* Candidate Fourier bin in miniFFT */
double mini_power; /* Candidate normalized power in miniFFT */
double mini_numsum; /* Number of powers summed to get candidate */
double mini_sigma; /* Equivalent candidate sigma (for sum pow) */
double psr_p; /* Approx PSR period (miniFFT center bin) */
double orb_p; /* Approx orbital period (s) */
} rawbincand;
typedef struct RZWERRS {
double p; /* Pulsation period at start of run (s) */
double pd; /* Period derivative at start of run (s/s) */
double pdd; /* Period 2nd derivative (s/s^2) */
double f; /* Pulsation frequency at start of run (hz) */
double fd; /* Frequency derivative at start of run (hz/s) */
double fdd; /* Frequency 2nd derivative (hz/s^2) */
float perr; /* Error in pulsation period (s) */
float pderr; /* Error in period derivative (s/s) */
float pdderr; /* Error in 2nd period derivative (s/s^2) */
float ferr; /* Error in frequency (hz) */
float fderr; /* Error in frequency derivative (hz/s) */
float fdderr; /* Error in 2nd frequency derivative (hh/s^2) */
} rzwerrs;
#ifndef ORBITPARAMS_TYPE
typedef struct orbitparams {
double p; /* Orbital period (s) */
double e; /* Orbital eccentricity */
double x; /* Projected semi-major axis (lt-sec) */
double w; /* Longitude of periapsis (deg) */
double t; /* Time since last periastron passage (s) */
double pd; /* Orbital period derivative (s/yr) */
double wd; /* Advance of longitude of periapsis (deg/yr) */
} orbitparams;
#define ORBITPARAMS_TYPE 1
#endif
typedef struct foldstats {
double numdata; /* Number of data bins folded */
double data_avg; /* Average level of the data bins */
double data_var; /* Variance of the data bins */
double numprof; /* Number of bins in the profile */
double prof_avg; /* Average level of the profile bins */
double prof_var; /* Variance of the profile bins */
double redchi; /* Reduced chi-squared of the profile */
} foldstats;
typedef struct bird{
double lobin;
double hibin;
} bird;
/***** Function Prototypes *****/
/* From swapendian.c: */
float swap_float(float var);
double swap_double(double var);
long double swap_longdouble(long double var);
long long swap_longlong(long long var);
int swap_int(int var);
unsigned int swap_uint(unsigned int var);
short swap_short(short var);
unsigned short swap_ushort(unsigned short var);
/* From responses.c: */
int r_resp_halfwidth(presto_interp_acc accuracy);
/* Return the approximate kernel half width in FFT bins required */
/* to achieve a fairly high accuracy correlation based correction */
/* or interpolation for a standard Fourier signal. */
/* Arguments: */
/* 'accuracy' is either LOWACC or HIGHACC. */
/* Notes: */
/* The result must be multiplied by 2*'numbetween' to get the */
/* length of the array required to hold such a kernel. */
int z_resp_halfwidth(double z, presto_interp_acc accuracy);
/* Return the approximate kernel half width in FFT bins required */
/* to achieve a fairly high accuracy correlation based correction */
/* or interpolation for a Fourier signal with constant f-dot. (i.e */
/* a constant frequency derivative) */
/* Arguments: */
/* 'z' is the Fourier Frequency derivative (# of bins the signal */
/* smears over during the observation). */
/* 'accuracy' is either LOWACC or HIGHACC. */
/* Notes: */
/* The result must be multiplied by 2*'numbetween' to get the */
/* length of the array required to hold such a kernel. */
int w_resp_halfwidth(double z, double w, presto_interp_acc accuracy);
/* Return the approximate kernel half width in FFT bins required */
/* to achieve a fairly high accuracy correlation based correction */
/* or interpolation for a Fourier signal with an f-dot that (i.e */
/* varies linearly in time -- a constant f-dotdot) */
/* Arguments: */
/* 'z' is the average Fourier Frequency derivative (# of bins */
/* the signal smears over during the observation). */
/* 'w' is the Fourier Frequency 2nd derivative (change in the */
/* Fourier f-dot during the observation). */
/* 'accuracy' is either LOWACC or HIGHACC. */
/* Notes: */
/* The result must be multiplied by 2*'numbetween' to get the */
/* length of the array required to hold such a kernel. */
void binary_velocity(double T, orbitparams * orbit,
double *minv, double *maxv);
/* Return the minimum and maximum orbital velocities of a pulsar */
/* during an observation as a fraction of the speed of light. */
/* Arguments: */
/* 'T' is the length of the observation in seconds. */
/* 'orbit' is a ptr to a orbitparams structure containing the */
/* Keplerian orbital parameters of the binary system. */
int bin_resp_halfwidth(double ppsr, double T, orbitparams * orbit);
/* Return the approximate kernel half width in FFT bins required */
/* to achieve a fairly high accuracy correlation based correction */
/* or interpolation for a pulsar in a binary orbit. */
/* Arguments: */
/* 'ppsr' is the period of the pusar in seconds. */
/* 'T' is the length of the observation in seconds. */
/* 'orbit' is a ptr to a orbitparams structure containing the */
/* Keplerian orbital parameters of the binary system. */
/* Notes: */
/* The result must be multiplied by 2 * 'numbetween' to get the */
/* length of the array required to hold such a kernel. */
fcomplex *gen_r_response(double roffset, int numbetween, int numkern);
/* Generate a complex response function for Fourier interpolation. */
/* Arguments: */
/* 'roffset' is the offset in Fourier bins for the full response */
/* (i.e. At this point, the response would equal 1.0) */
/* 'numbetween' is the number of points to interpolate between */
/* each standard FFT bin. (i.e. 'numbetween' = 1 = interbins) */
/* 'numkern' is the number of complex points that the kernel will */
/* contain. */
fcomplex *gen_z_response(double roffset, int numbetween, double z, \
int numkern);
/* Generate the response function for Fourier f-dot interpolation. */
/* Arguments: */
/* 'roffset' is the offset in Fourier bins for the full response */
/* (i.e. At this point, the response would equal 1.0) */
/* 'numbetween' is the number of points to interpolate between */
/* each standard FFT bin. (i.e. 'numbetween' = 1 = interbins) */
/* 'z' is the Fourier Frequency derivative (# of bins the signal */
/* smears over during the observation). */
/* 'numkern' is the number of complex points that the kernel will */
/* contain. */
fcomplex *gen_w_response(double roffset, int numbetween, double z, \
double w, int numkern);
/* Generate the response function for Fourier f-dot interpolation. */
/* Arguments: */
/* 'roffset' is the offset in Fourier bins for the full response */
/* (i.e. At this point, the response would equal 1.0) */
/* 'numbetween' is the number of points to interpolate between */
/* each standard FFT bin. (i.e. 'numbetween' = 1 = interbins) */
/* 'z' is the average Fourier Frequency derivative (# of bins */
/* the signal smears over during the observation). */
/* 'w' is the Fourier Frequency 2nd derivative (change in the */
/* Fourier f-dot during the observation). */
/* 'numkern' is the number of complex points that the kernel will */
/* contain. */
fcomplex *gen_w_response2(double roffset, int numbetween, double z, \
double w, int numkern);
/* Generate the response function for Fourier f-dot interpolation. */
/* Arguments: */
/* 'roffset' is the offset in Fourier bins for the full response */
/* (i.e. At this point, the response would equal 1.0) */
/* 'numbetween' is the number of points to interpolate between */
/* each standard FFT bin. (i.e. 'numbetween' = 1 = interbins) */
/* 'z' is the average Fourier Frequency derivative (# of bins */
/* the signal smears over during the observation). */
/* 'w' is the Fourier Frequency 2nd derivative (change in the */
/* Fourier f-dot during the observation). */
/* 'numkern' is the number of complex points that the kernel will */
/* contain. */
fcomplex *gen_bin_response(double roffset, int numbetween, double ppsr, \
double T, orbitparams * orbit, int numkern);
/* Generate the Fourier response function for a sinusoidal PSR */
/* signal from a binary orbit. */
/* Arguments: */
/* 'roffset' is the offset in Fourier bins for the full response */
/* (i.e. At this point, the response would equal 1.0) */
/* 'numbetween' is the number of points to interpolate between */
/* each standard FFT bin. (i.e. 'numbetween' = 1 = interbins) */
/* 'ppsr' is the period of the pusar in seconds. */
/* 'T' is the length of the observation in seconds. */
/* 'orbit' is a ptr to a orbitparams structure containing the */
/* Keplerian orbital parameters of the binary system. */
/* 'numkern' is the number of complex points that the kernel will */
/* contain. */
/* From characteristics.c: */
/* Routines to calculate characteristics of a candidate signal */
float get_numphotons(FILE * file);
/* Return the total number of photons in the FFT file */
/* i.e. it returns the value of the 0th frequency bin. */
/* Arguments: */
/* 'file' is a pointer to the file you want to access. */
double get_localpower(fcomplex *data, long numdata, double r);
/* Return the local power level at specific FFT frequency. */
/* Arguments: */
/* 'data' is a pointer to a complex FFT. */
/* 'numdata' is the number of complex points in 'data'. */
/* 'r' is the Fourier frequency in data that we want to */
/* interpolate. */
double get_localpower3d(fcomplex *data, long numdata, double r, \
double z, double w);
/* Return the local power level around a specific FFT */
/* frequency, f-dot, and f-dotdot. */
/* Arguments: */
/* 'data' is a pointer to a complex FFT. */
/* 'numdata' is the number of complex points in 'data'. */
/* 'r' is the Fourier frequency in data that we want to */
/* interpolate. */
/* 'z' is the Fourier Frequency derivative (# of bins the */
/* signal smears over during the observation). */
/* 'w' is the Fourier Frequency 2nd derivative (change in the */
/* Fourier f-dot during the observation). */
void get_derivs3d(fcomplex *data, long numdata, double r, \
double z, double w, double localpower, \
rderivs *result);
/* Return an rderives structure that contains the power, */
/* phase, and their first and second derivatives at a point */
/* in the F/F-dot/F-dortdot volume. */
/* Arguments: */
/* 'data' is a pointer to a complex FFT. */
/* 'numdata' is the number of complex points in 'data'. */
/* 'r' is the Fourier frequency in data that we want to */
/* interpolate. */
/* 'z' is the Fourier Frequency derivative (# of bins the */
/* signal smears over during the observation). */
/* 'w' is the Fourier Frequency 2nd derivative (change in */
/* the Fourier f-dot during the observation). */
/* 'localpower' is the local power level around the signal. */
/* 'result' is a pointer to an rderivs structure that will */
/* contain the results. */
void calc_props(rderivs data, double r, double z, double w, \
fourierprops * result);
/* Return a fourierprops structure that contains the various */
/* properties of a signal described by Middleditch, Deich, */
/* and Kulkarni in _Isolated_Pulsars_, 1993, p372. */
/* Arguments: */
/* 'data' is a pointer to an rderivs structure containing */
/* derivative information about the peak in question. */
/* 'r' is the Fourier frequency in data that we want to */
/* interpolate. */
/* 'z' is the Fourier Frequency derivative (# of bins the */
/* signal smears over during the observation). */
/* 'w' is the Fourier Frequency second derivative. */
/* 'result' is a pointer to an fourierprops structure that */
/* will contain the results. */
void calc_binprops(fourierprops * props, double T, int lowbin, \
int nfftbins, binaryprops * result);
/* Return a binaryprops structure that contains the various */
/* estimates of the binary pulsar system from a mini-FFT. */
/* Arguments: */
/* 'props' is a pointer to the candidate's fourierprops. */
/* 'T' is the total length (sec) of the original time series. */
/* 'lowbin' is the Fourier bin number from the original FFT */
/* the lowest bin in the mini-FFT. */
/* 'nfftbins' is the number of bins in the mini-FFT. */
/* 'absnorm' is the value of the power normalization */
/* constant for this mini-FFT. */
/* 'result' is the returned binaryprops structure. */
void calc_rzwerrs(fourierprops *props, double T, rzwerrs *result);
/* Calculate periods, frequencies, their derivatives */
/* and their errors. */
/* Arguments: */
/* 'props' is a pointer to a fourierprops structure. */
/* 'T' is the length of the data set in sec (i.e. N*dt). */
/* 'result' is a pointer to the returned rzwerrs struct. */
double equivalent_gaussian_sigma(double logp);
/* Return the approximate significance in Gaussian sigmas */
/* corresponding to a natural log probability logp */
double chi2_logp(double chi2, double dof);
/* Return the natural log probability corresponding to a chi^2 value */
/* of chi2 given dof degrees of freedom. */
double chi2_sigma(double chi2, double dof);
/* Return the approximate significance in Gaussian sigmas */
/* sigmas of a chi^2 value of chi2 given dof degrees of freedom. */
double candidate_sigma(double power, int numsum, double numtrials);
/* Return the approximate significance in Gaussian */
/* sigmas of a candidate of numsum summed powers, */
/* taking into account the number of independent trials. */
double power_for_sigma(double sigma, int numsum, double numtrials);
/* Return the approximate summed power level required */
/* to get a Gaussian significance of 'sigma', taking */
/* into account the number of independent trials. */
double chisqr(double *data, int numdata, double avg, double var);
/* Calculates the chi-square of the 'data' which has average */
/* 'avg', and variance 'var'. */
double z2n(double *data, int numdata, double var, int n);
/* Calculates the chi-square of the 'data' which has variance 'var'. */
void switch_f_and_p(double in, double ind, double indd,
double *out, double *outd, double *outdd);
/* Convert p, p-dot, and p-dotdot into f, f-dot, */
/* and f-dotdot or vise-versa. */
/* dispersion.c: */
/* Functions to de-disperse data */
double tree_max_dm(int numchan, double dt, double lofreq, double hifreq);
/* Return the maximum Dispersion Measure (dm) in cm-3 pc, the */
/* tree de-dispersion technique can correct for given a sample */
/* interval 'dt', the number of channels 'numchan', and the */
/* low and high observation frequencies in MHz. */
double smearing_from_bw(double dm, double center_freq, double bandwidth);
/* Return the smearing in seconds caused by dispersion, given */
/* a Dispersion Measure (dm) in cm-3 pc, the central frequency */
/* and the bandwith of the observation in MHz. */
double delay_from_dm(double dm, double freq_emitted);
/* Return the delay in seconds caused by dispersion, given */
/* a Dispersion Measure (dm) in cm-3 pc, and the emitted */
/* frequency (freq_emitted) of the pulsar in MHz. */
double dm_from_delay(double delay, double freq_emitted);
/* Return the Dispersion Measure in cm-3 pc, that would */
/* cause a pulse emitted at frequency 'freq_emitted' to be */
/* delayed by 'delay' seconds. */
double *dedisp_delays(int numchan, double dm, double lofreq,
double chanwidth, double voverc);
/* Return an array of delays (sec) for dedispersing 'numchan' */
/* channels at a DM of 'dm'. 'lofreq' is the center frequency */
/* in MHz of the lowest frequency channel. 'chanwidth' is the */
/* width in MHz of each channel. 'voverc' is the observatory's */
/* velocity towards or away from the source. This is to adjust */
/* the frequencies for doppler effects (for no correction use */
/* voverc=0). The returned array is allocated by this routine. */
void dedisp(unsigned char *data, unsigned char *lastdata, int numpts,
int numchan, double *dispdelays, float *result);
/* De-disperse a stretch of data with numpts * numchan points. */
/* The delays (in bins) are in dispdelays for each channel. */
/* The result is returned in result. The input data and */
/* dispdelays are always in ascending frequency order. */
/* Input data are ordered in time, with the channels stored */
/* together at each time point. */
double *subband_delays(int numchan, int numsubbands, double dm,
double lofreq, double chanwidth,
double voverc);
/* Return an array of delays (sec) for the highest frequency */
/* channels of each subband used in a subband de-dispersion. */
/* These are the delays described in the 'Note:' in the */
/* description of subband_search_delays(). See the comments */
/* for dedisp_delays() for more info. */
double *subband_search_delays(int numchan, int numsubbands, double dm,
double lofreq, double chanwidth,
double voverc);
/* Return an array of delays (sec) for a subband DM search. The */
/* delays are calculated normally for each of the 'numchan' channels */
/* using the appropriate frequencies at the 'dm'. Then the delay */
/* from the highest frequency channel of each of the 'numsubbands' */
/* subbands is subtracted from each subband. This gives the subbands */
/* the correct delays for each freq in the subband, but the subbands */
/* themselves are offset as if they had not been de-dispersed. This */
/* way, we can call float_dedisp() on the subbands if needed. */
/* 'lofreq' is the center frequency in MHz of the lowest frequency */
/* channel. 'chanwidth' is the width in MHz of each channel. The */
/* returned array is allocated by this routine. 'voverc' is used to */
/* correct the input frequencies for doppler effects. See the */
/* comments in dedisp_delays() for more info. */
/* Note: When performing a subband search, the delays for each */
/* subband must be calculated with the frequency of the highest */
/* channel in each subband, _not_ the center subband frequency. */
void dedisp_subbands(float *data, float *lastdata,
int numpts, int numchan,
int *delays, int numsubbands, float *result);
// De-disperse a stretch of data with numpts * numchan points into
// numsubbands subbands. Each time point for each subband is a float
// in the result array. The result array order is subbands of
// increasing frequency together at each time pt. The delays (in
// bins) are in delays for each channel. The input data and
// dispdelays are always in ascending frequency order. Input data are
// ordered in time, with the channels stored together at each time
// point.
void float_dedisp(float *data, float *lastdata,
int numpts, int numchan,
int *delays, float approx_mean, float *result);
// De-disperse a stretch of data with numpts * numchan points. The
// delays (in bins) are in delays for each channel. The result is
// returned in result. The input data and delays are always in
// ascending frequency order. Input data are ordered in time, with
// the channels stored together at each time point.
void combine_subbands(double *inprofs, foldstats *stats,
int numparts, int numsubbands, int proflen,
int *delays, double *outprofs,
foldstats *outprofstats);
/* Combine 'nparts' sets of 'numsubbands' profiles, each of length */
/* 'proflen' into a 'nparts' de-dispersed profiles. The de-dispersion */
/* uses the 'delays' (of which there are 'numsubbands' many) to */
/* show how many bins to shift each profile to the right. Only */
/* positive numbers may be used (left shifts may be accomplished using */
/* the shift modulo 'proflen'). The 'stats' about the profiles are */
/* combined as well and the combined stats are returned in */
/* 'outprofstats'. All arrays must be pre-allocated. */
/* output.c: */
/* Functions for text-based ouput of information */
int nice_output_1(char *output, double val, double err, int len);
/* Generates a string in "output" of length len with "val" rounded */
/* to the appropriate decimal place and the error in parenthesis */
/* as in scientific journals. The error has 1 decimal place. */
/* Note: len should be ~ 20 to show full double precision */
/* if the base 10 exponent of the error needs to be shown. */
/* If len == 0, left-justified minimum length string is returned. */
/* If len > 0, the string returned has is right justified. */
int nice_output_2(char *output, double val, double err, int len);
/* Generates a string in "output" of length len with "val" rounded */
/* to the appropriate decimal place and the error in parenthesis */
/* as in scientific journals. The error has 2 decimal places. */
/* Note: len should be ~ 20 to show full double precision */
/* if the base 10 exponent of the error needs to be shown. */
/* If len == 0, left-justified minimum length string is returned. */
/* If len > 0, the string returned has is right justified. */
void print_candidate(fourierprops * cand, double dt, long N, \
double nph, int numerrdigits);
/* Outputs a 2 column summary of all the properties or a fourier peak */
void print_bin_candidate(binaryprops * cand, int numerrdigits);
/* Outputs a 2 column summary of all the properties or a fourier peak */
void file_reg_candidates(fourierprops cand[], char *notes, int numcands, \
double dt, long N, double nph,
char name[], char longname[]);
/* Outputs a .ps file describing all the candidates from a search. */
void file_bin_candidates(binaryprops cand[], char *notes, \
int numcands, char name[]);
/* Outputs a .ps file describing all the binary candidates from a */
/* binary search.*/
/* get_candidates.c: */
/* Functions for manipulating candidate files */
int read_rzw_cand(FILE *file, fourierprops *cands);
/* Read the next rzw candidate from the file */
/* If successful, return 1, else 0 */
int read_bin_cand(FILE *file, binaryprops *cands);
/* Read the next binary candidate from the file */
/* If successful, return 1, else 0 */
int read_rawbin_cand(FILE *file, rawbincand *cands);
/* Read the next rawbin candidate from the file */
/* If successful, return 1, else 0 */
void get_rzw_cand(char *filenm, int candnum, fourierprops * cand);
/* Read the rzw candidate file 'filenm' and return a */
/* pointer to the fourierprops that describes it. */
void get_bin_cand(char *filenm, int candnum, binaryprops * cand);
/* Read the bin candidate file 'filenm' and return a */
/* pointer to the binaryprops that describes it. */
void get_rawbin_cand(char *filenm, int candnum, rawbincand * cand);
/* Read the rawbin candidate file 'filenm' and return a */
/* pointer to the rawbincand that describe it. */
/* read_fft.c: */
/* Functions for getting information from an FFT file */
fcomplex *read_fcomplex_file(FILE *file, long firstpt, long numpts);
/* Return an fcomplex vector with complex data taken from a file. */
/* Argumants: */
/* 'file' is a pointer to the file you want to access. */
/* 'firstpt' is the number of the first point to get. (0 = 1st */
/* point in the file). If < 0, the resulting array will */
/* be zero padded. */
/* 'numpts' is the number of points to get from the file. */
/* If the number of bins to read takes us past the end of */
/* file, the returned vector will be zero padded. */
float *read_float_file(FILE *file, long firstpt, long numpts);
/* Return a float vector with complex data taken from a file. */
/* Argumants: */
/* 'file' is a pointer to the file you want to access. */
/* 'firstpt' is the number of the first point to get. (0 = 1st */
/* point in the file). If < 0, the resulting array will */
/* be zero padded. */
/* 'numpts' is the number of points to get from the file. */
/* If the number of bins to read takes us past the end of */
/* file, the returned vector will be zero padded. */
/* The following routines are used by the routines above to do */
/* their thing.. */
/* In select.c */
int prune_powers(float *arr, int n, int numsumpows);
/* Sets powers that are more than approx PRUNELEV standard */
/* devs above the median value to NEWLEV times the median */
/* value. This helps 'clean' the spectrum of high power */
/* signals that probably have nothing to do with a phase */
/* modulation spectrum (i.e. they are RF noise or strong */
/* solitary pulsars. */
void hpselect(unsigned long m, unsigned long n, \
float arr[], powindex heap[]);
/* Selects the m largest values from the array arr */
/* and stores them and their indices in heap and heapindex. */
/* In median.c */
float median(float arr[], int n);
/* Finds the median (but messes up the array order) */
int comp_psr_to_cand(fourierprops * cand, infodata * idata, char *output, \
int full);
/* Compares a pulsar candidate defined by its properties found in */
/* *cand, and *idata with all of the pulsars in the pulsar */
/* database. It returns a string (verbose if full==1) describing */
/* the results of the search in *output. */
int comp_bin_to_cand(binaryprops * cand, infodata * idata, \
char *output, int full);
/* Compares a binary PSR candidate defined by its props found in */
/* *cand, and *idata with all of the pulsars in the pulsar */
/* database. It returns a string (verbose if full==1) describing */
/* the results of the search in *output. */
double dms2rad(int deg, int min, double sec);
/* Convert degrees, minutes and seconds into radians */
double hms2rad(int hour, int min, double sec);
/* Convert hours, minutes and seconds into radians */
double sphere_ang_diff(double ra1, double dec1, double ra2, double dec2);
/* Return the spherical angle (radians) between two RA and DECS */
/* In sorter.c */
float percolate(position * list, int nlist, int spot);
/* Pushes a position structure as far up a sorted list of positions */
/* as it needs to go to keep the list sorted. Returns the new low */
/* power in the list. */
float percolate_bin(binaryprops * list, int nlist);
/* Pushes a binaryprops structure as far up a sorted list of structs*/
/* as it needs to go to keep the list sorted. Returns the new low */
/* power in the list. */
/* From prep_corr.c */
int next_good_fftlen(int N);
/* Return one of the shortest, yet best performing, FFT lengths larger
* than N. This assumes FFTW. */
int fftlen_from_kernwidth(int kernwidth);
/* return the length of the optimal FFT to use for correlations with
* some kernel width kernwidth. This assumes FFTW. */
void spread_with_pad(fcomplex *data, int numdata, \
fcomplex *result, int numresult, \
int numbetween, int numpad);
/* Prepare the data array for correlation by spreading */
/* the input data array and padding it. */
/* Arguments: */
/* 'data' is the FFT array to be prepared */
/* 'numdata' is the number of complex points in 'data' */
/* 'result' is the prepped data array */
/* 'numresult' is the number of complex points in 'result' */
/* 'numbetween' is the number of interpolated pts per bin */
/* 'numpad' is the number of bins to use as zero padding */
void spread_no_pad(fcomplex *data, int numdata, \
fcomplex *result, int numresult, \
int numbetween);
/* Prepare the data array for correlation by spreading */
/* the input data array. */
/* Arguments: */
/* 'data' is the FFT array to be prepared */
/* 'numdata' is the number of complex points in 'data' */
/* 'result' is the prepped data array */
/* 'numresult' is the number of complex points in 'result' */
/* 'numbetween' is the number of interpolated pts per bin */
void paddata(fcomplex *data, int numdata, int numpad);
/* Pad the last 'numpad' bins of 'data' with zeros. */
/* Arguments: */
/* 'data' is the FFT array to be padded */
/* 'numdata' is the number of complex points in 'data' */
/* 'numpad' is the number of bins to use as zero padding */
void place_complex_kernel(fcomplex *kernel, int numkernel, \
fcomplex *result, int numresult);
/* This routine places the kernel in a zero filled array */
/* with half of the response at the beginning and half */
/* of the response at the end of the result array. See */
/* Numerical Recipes in C 2ed, p 541 for more info. */
/* Arguments: */
/* 'kernel' is a complex response function. Bin zero */
/* response is in bin numkernel/2. */
/* 'numkernel' is the number of points in the kernel. */
/* This should be an even number. */
/* 'result' is the result array. */
/* 'numresult' is the number of points in the result. */
void place_real_kernel(float *kernel, int numkernel, \
float *result, int numresult);
/* This routine places the kernel in a zero filled array */
/* with half of the response at the beginning and half */
/* of the response at the end of the result array. See */
/* Numerical Recipes in C 2ed, p 541 for more info. */
/* Arguments: */
/* 'kernel' is a real-valued response function. Bin */
/* zero response is in bin numkernel/2. */
/* 'numkernel' is the number of points in the kernel. */
/* This should be an even number. */
/* 'result' is the result array. */
/* 'numresult' is the number of points in the result. */
void chop_complex_ends(fcomplex *data, int numdata, \
fcomplex *result, int numresult, \
int chopbins);
/* Chop the contaminated ends off of an array that has */
/* been correlated/convolved. */
/* Arguments: */
/* 'data' is the array to chop. */
/* 'numdata' is the number of points in data. */
/* 'result' is the resultant array. */
/* 'numresult' is the number of points in the result. */
/* 'chopbins' is the number of bins to chop on each */
/* end of the data array. */
void chop_real_ends(float *data, int numdata, \
float *result, int numresult, \
int chopbins);
/* Chop the contaminated ends off of an array that has */
/* been correlated/convolved. */
/* Arguments: */
/* 'data' is the array to chop. */
/* 'numdata' is the number of points in data. */
/* 'result' is the resultant array. */
/* 'numresult' is the number of points in the result. */
/* 'chopbins' is the number of bins to chop on each */
/* end of the data array. */
/* In correlations.c */
fcomplex *complex_corr_conv(fcomplex *data, fcomplex *kernel, \
int numdata, \
presto_ffts ffts, presto_optype type);
/* Perform and return a complex correlation or convolution. */
/* Arguments: */
/* 'data' is the complex array to correlate/convolve. */
/* 'kernel' is the correlation/convolution kernel. */
/* 'numdata' is the length of 'data', 'kernel' and the result. */
/* 'ffts' describes how to perform the convolution/correlation. */
/* 'ffts' = FFTDK: FFT both the 'data' and the 'kernel'. */
/* 'ffts' = FFTD: FFT only the 'data' not the 'kernel'. */
/* 'ffts' = FFTK: FFT only the 'kernel' not the 'data'. */
/* 'ffts' = NOFFTS: Don't FFT the 'data' or the 'kernel'. */
/* 'type' is the type of operation to perform. */
/* 'type' = CONV: return a convolution in a new vector. */
/* 'type' = CORR: return a correlation in a new vector. */
/* 'type' = INPLACE_CONV: convolution over-writes 'data'. */
/* 'type' = INPLACE_CORR: correlation over-writes 'data'. */
float *real_corr_conv(float *data, float *kernel, int numdata, \
presto_ffts ffts, presto_optype type);
/* Perform and return a real-valued correlation or convolution. */
/* Arguments: */
/* 'data' is the complex array to correlate/convolve. */
/* 'kernel' is the correlation/convolution kernel. */
/* 'numdata' is the length of 'data', 'kernel' and the result. */
/* 'ffts' describes how to perform the convolution/correlation. */
/* 'ffts' = FFTDK: FFT both the 'data' and the 'kernel'. */
/* 'ffts' = FFTD: FFT only the 'data' not the 'kernel'. */
/* 'ffts' = FFTK: FFT only the 'kernel' not the 'data'. */
/* 'ffts' = NOFFTS: Don't FFT the 'data' or the 'kernel'. */
/* 'type' is the type of operation to perform. */
/* 'type' = CONV: return a convolution in a new vector. */
/* 'type' = CORR: return a correlation in a new vector. */
/* 'type' = INPLACE_CONV: convolution over-writes 'data'. */
/* 'type' = INPLACE_CORR: correlation over-writes 'data'. */
/* In corr_routines.c */
int corr_complex(fcomplex *data, int numdata, presto_datainf datainf, \
fcomplex *kern, int numkern, presto_datainf kerninf, \
fcomplex *result, int numresult, int lobin, \
int numbetween, int kern_half_width, presto_optype optype);
/* This routine is a general correlation or convolution routine */
/* for complex data. It can perform convolutions or correlations */
/* on raw complex data, data that is prepared for a convolution/ */
/* correlation but not FFTd, or already FFTd data. The kernel */
/* that it uses can also be raw, prepped, or FFTd. If you call */
/* the routine multiple times with either the same kernel or data */
/* array, it uses a saved version of the array from the previous */
/* call to cut down on many processing steps. The return value */
/* tells how many usable (i.e. non-contaminated) points were */
/* returned in the result array (the first value will be that of */
/* 'lobin'). This routine will _not_ perform in-place */
/* correlations or convolutions (i.e. it ignores those choices */
/* for 'optype'). */
/* Arguments: */
/* 'data' is a complex array of the data to be interpolated. */
/* 'numdata' is the number of complex points in 'data'. */
/* 'datainf' is one of the following that describes the data: */
/* RAW = Normal un-altered complex data. */
/* PREPPED = Data has been padded and spread based */
/* on 'kern_half_width' and 'numbetween' */
/* and is ready to be FFTd. */
/* FFT = Data has already been prepared and FFTd. */
/* SAME = Data is the same as the previous call. */
/* The routine uses its saved data. */
/* 'kern' is the correlation kernel. */
/* 'numkern' is the number of complex points in 'kern'. */
/* 'kerninf' is one of the same choices as 'datainf' above. */
/* 'result' is the resulting complex array (must already exist). */
/* 'numresult' is the number of complex points in 'result'. */
/* 'lobin' is the lowest fourier bin to convolve/correlate. */
/* 'numbetween' is the number of bins to spread the data points. */
/* 'kern_half_width' is half the width (bins) of the raw kernel. */
/* 'optype' is either CORR or CONV (correlation or convolution). */
/* Notes: */
/* If either 'datainf' or 'kerninf' are of type PREPPED or FFT, */
/* then the length of the FFTs used in the correlation/ */
/* convolution calculations will be of length 'numdata' or */
/* 'numkern'. If both 'datainf' and 'kerninf' are of type */
/* PREPPED or FFT then 'numdata' and 'numkern' must have the */
/* same value. In order for SAME values of 'datainf' and */
/* 'kerninf' to help out, the routine must be called with the */
/* same values for 'kern_half_width' and 'numbetween' as well. */
void stretch_fft(fcomplex *data, int numdata, \
fcomplex *result, int numresult);
/* This routine stretches and/or interpolates an FFT of length */
/* numdata. It zeros 'result' where end-effects have occurred. */
/* This routine is usually used to co-add stretched FFTs to */
/* increase the signal-to-noise ratios of a detection. */
/* Arguments: */
/* 'data' is a pointer to a complex FFT. */