diff --git a/clig/accelsearch.1 b/clig/accelsearch.1 index 1496b6367..32c9873f9 100644 --- a/clig/accelsearch.1 +++ b/clig/accelsearch.1 @@ -11,7 +11,7 @@ .\" it edited by clig, remove the respective pair of cligPart-lines. .\" .\" cligPart TITLE -.TH "accelsearch" 1 "07Dec16" "Clig-manuals" "Programmer's Manual" +.TH "accelsearch" 1 "19Sep17" "Clig-manuals" "Programmer's Manual" .\" cligPart TITLE end .\" cligPart NAME @@ -26,6 +26,7 @@ accelsearch \- Search an FFT or short time series for pulsars using a Fourier do [-lobin lobin] [-numharm numharm] [-zmax zmax] +[-wmax wmax] [-sigma sigma] [-rlo rlo] [-rhi rhi] @@ -69,6 +70,10 @@ The max (+ and -) Fourier freq deriv to search, 1 Int value between 0 and 1200. .br Default: `200' +.IP -wmax +The max (+ and -) Fourier freq double derivs to search, +.br +1 Int value between 0 and 4000. .IP -sigma Cutoff sigma for choosing candidates, .br diff --git a/clig/accelsearch_cmd.cli b/clig/accelsearch_cmd.cli index 962cb7348..51e12b0be 100644 --- a/clig/accelsearch_cmd.cli +++ b/clig/accelsearch_cmd.cli @@ -18,6 +18,8 @@ Int -numharm numharm {The number of harmonics to sum (power-of-two)}\ -r 1 32 -d 8 Int -zmax zmax {The max (+ and -) Fourier freq deriv to search} \ -r 0 1200 -d 200 +Int -wmax wmax {The max (+ and -) Fourier freq double derivs to search} \ + -r 0 4000 Float -sigma sigma {Cutoff sigma for choosing candidates}\ -r 1.0 30.0 -d 2.0 Double -rlo rlo {The lowest Fourier frequency (of the highest harmonic!) to search} \ diff --git a/docs/accelsearch.1 b/docs/accelsearch.1 index 1496b6367..32c9873f9 100644 --- a/docs/accelsearch.1 +++ b/docs/accelsearch.1 @@ -11,7 +11,7 @@ .\" it edited by clig, remove the respective pair of cligPart-lines. .\" .\" cligPart TITLE -.TH "accelsearch" 1 "07Dec16" "Clig-manuals" "Programmer's Manual" +.TH "accelsearch" 1 "19Sep17" "Clig-manuals" "Programmer's Manual" .\" cligPart TITLE end .\" cligPart NAME @@ -26,6 +26,7 @@ accelsearch \- Search an FFT or short time series for pulsars using a Fourier do [-lobin lobin] [-numharm numharm] [-zmax zmax] +[-wmax wmax] [-sigma sigma] [-rlo rlo] [-rhi rhi] @@ -69,6 +70,10 @@ The max (+ and -) Fourier freq deriv to search, 1 Int value between 0 and 1200. .br Default: `200' +.IP -wmax +The max (+ and -) Fourier freq double derivs to search, +.br +1 Int value between 0 and 4000. .IP -sigma Cutoff sigma for choosing candidates, .br diff --git a/include/accel.h b/include/accel.h index 3c4a9343d..b723c7868 100644 --- a/include/accel.h +++ b/include/accel.h @@ -25,94 +25,112 @@ #define ACCEL_DZ 2 /* Reciprocal of ACCEL_DZ */ #define ACCEL_RDZ 0.5 +/* Stepsize in Fourier F-dot-dot */ +#define ACCEL_DW 20 +/* Reciprocal of ACCEL_DW */ +#define ACCEL_RDW 0.05 /* Closest candidates we will accept as independent */ #define ACCEL_CLOSEST_R 15.0 /* Padding for .dat file reading so that we don't SEGFAULT */ #define ACCEL_PADDING 2000 typedef struct accelobs{ - long long N; /* Number of data points in observation */ - long long numbins; /* Number of spectral bins in the file */ - long long lobin; /* Lowest spectral bin present in the file */ - long long highestbin;/* Highest spectral bin present in the file */ - int fftlen; /* Length of short FFTs to us in search */ - int numharmstages; /* Number of stages of harmonic summing */ - int numz; /* Number of f-dots searched */ - int numbetween; /* Highest fourier freq resolution (2=interbin) */ - int numzap; /* Number of birdies to zap */ - int dat_input; /* The input file is a short time series */ - int mmap_file; /* The file number if using MMAP */ - int inmem; /* True if we want to keep the full f/fdot plan in RAM */ - int norm_type; /* 0 = old-style block median, 1 = local-means power norm */ - double dt; /* Data sample length (s) */ - double T; /* Total observation length */ - double rlo; /* Minimum fourier freq to search */ - double rhi; /* Maximum fourier freq to search */ - double dr; /* Stepsize in fourier freq (1/numbetween) */ - double zlo; /* Minimum fourier fdot to search */ - double zhi; /* Maximum fourier fdot to search */ - double dz; /* Stepsize in fourier fdot */ - double baryv; /* Average barycentric velocity during observation */ - float nph; /* Freq 0 level if requested, 0 otherwise */ - float sigma; /* Cutoff sigma to choose a candidate */ - float *powcut; /* Cutoff powers to choose a cand (per harmsummed) */ - float *ffdotplane; /* The full F-Fdot plane if working in memory */ - double *lobins; /* The low Fourier freq boundaries to zap (RFI) */ - double *hibins; /* The high Fourier freq boundaries to zap (RFI) */ - long long *numindep; /* Number of independent spectra (per harmsummed) */ - FILE *fftfile; /* The FFT file that we are analyzing */ - FILE *workfile; /* A text file with candidates as they are found */ - fcomplex *fft; /* A pointer to the FFT for MMAPing or input time series */ - char *rootfilenm; /* The root filename for associated files. */ - char *candnm; /* The fourierprop save file for the fundamentals */ - char *accelnm; /* The filename of the final candidates in text */ - char *workfilenm; /* The filename of the working candidates in text */ - int use_harmonic_polishing; /* Should we force harmonics to be related */ + long long N; /* Number of data points in observation */ + long long numbins; /* Number of spectral bins in the file */ + long long lobin; /* Lowest spectral bin present in the file */ + long long highestbin;/* Highest spectral bin present in the file */ + int maxkernlen; /* Maximum full width (in points, not Fourier bins) of corr kernels */ + int corr_uselen; /* Number of good data points we will get from high-harm correlations */ + int fftlen; /* Length of short FFTs to us in search */ + int numharmstages; /* Number of stages of harmonic summing */ + int numz; /* Number of f-dots searched */ + int numw; /* Number of f-dot-dots searched */ + int numbetween; /* Highest fourier freq resolution (2=interbin) */ + int numzap; /* Number of birdies to zap */ + int dat_input; /* The input file is a short time series */ + int mmap_file; /* The file number if using MMAP */ + int inmem; /* True if we want to keep the full f-fdot plane in RAM */ + int norm_type; /* 0 = old-style block median, 1 = local-means power norm */ + double dt; /* Data sample length (s) */ + double T; /* Total observation length */ + double rlo; /* Minimum fourier freq to search */ + double rhi; /* Maximum fourier freq to search */ + double dr; /* Stepsize in fourier freq (1/numbetween) */ + double zlo; /* Minimum fourier fdot to search */ + double zhi; /* Maximum fourier fdot to search */ + double dz; /* Stepsize in fourier fdot */ + double wlo; /* Minimum fourier f-dot-dot to search */ + double whi; /* Maximum fourier f-dot-dot to search */ + double dw; /* Stepsize in fourier f-dot-dot */ + double baryv; /* Average barycentric velocity during observation */ + float nph; /* Freq 0 level if requested, 0 otherwise */ + float sigma; /* Cutoff sigma to choose a candidate */ + float *powcut; /* Cutoff powers to choose a cand (per harmsummed) */ + float *ffdotplane; /* The full f-fdot-fdotdot plane if working in memory */ + double *lobins; /* The low Fourier freq boundaries to zap (RFI) */ + double *hibins; /* The high Fourier freq boundaries to zap (RFI) */ + long long *numindep; /* Number of independent spectra (per harmsummed) */ + FILE *fftfile; /* The FFT file that we are analyzing */ + FILE *workfile; /* A text file with candidates as they are found */ + fcomplex *fft; /* A pointer to the FFT for MMAPing or input time series */ + char *rootfilenm; /* The root filename for associated files. */ + char *candnm; /* The fourierprop save file for the fundamentals */ + char *accelnm; /* The filename of the final candidates in text */ + char *workfilenm; /* The filename of the working candidates in text */ + int use_harmonic_polishing; /* Should we force harmonics to be related */ } accelobs; typedef struct accelcand{ - float power; /* Summed power level (normalized) */ - float sigma; /* Equivalent sigma based on numindep (above) */ - int numharm; /* Number of harmonics summed */ - double r; /* Fourier freq of first harmonic */ - double z; /* Fourier f-dot of first harmonic */ - double *pows; /* Optimized powers for the harmonics */ - double *hirs; /* Optimized freqs for the harmonics */ - double *hizs; /* Optimized fdots for the harmonics */ - rderivs *derivs; /* An rderivs structure for each harmonic */ + float power; /* Summed power level (normalized) */ + float sigma; /* Equivalent sigma based on numindep (above) */ + int numharm; /* Number of harmonics summed */ + double r; /* Fourier freq of first harmonic */ + double z; /* Fourier f-dot of first harmonic */ + double w; /* Fourier f-dot-dot of first harmonic */ + double *pows; /* Optimized powers for the harmonics */ + double *hirs; /* Optimized freqs for the harmonics */ + double *hizs; /* Optimized f-dots for the harmonics */ + double *hiws; /* Optimized f-dot-dots for the harmonics */ + rderivs *derivs; /* An rderivs structure for each harmonic */ } accelcand; typedef struct kernel{ - int z; /* The fourier f-dot of the kernel */ - int fftlen; /* Number of complex points in the kernel */ - int numgoodbins; /* The number of good points you can get back */ - int numbetween; /* Fourier freq resolution (2=interbin) */ - int kern_half_width; /* Half width (bins) of the raw kernel. */ - fcomplex *data; /* The FFTd kernel itself */ + int z; /* The fourier f-dot of the kernel */ + int w; /* The fourier f-dot-dot of the kernel */ + int fftlen; /* Number of complex points in the kernel */ + int numgoodbins; /* The number of good points you can get back */ + int numbetween; /* Fourier freq resolution (2=interbin) */ + int kern_half_width; /* Half width (bins) of the raw kernel. */ + fcomplex *data; /* The FFTd kernel itself */ } kernel; typedef struct subharminfo{ - int numharm; /* The number of sub-harmonics */ - int harmnum; /* The sub-harmonic number (fundamental = numharm) */ - int zmax; /* The maximum Fourier f-dot for this harmonic */ - int numkern; /* Number of kernels in the vector */ - kernel *kern; /* The kernels themselves */ - unsigned short *rinds; /* Table of indices for Fourier Freqs */ + int numharm; /* The number of sub-harmonics */ + int harmnum; /* The sub-harmonic number (fundamental = numharm) */ + int zmax; /* The maximum Fourier f-dot for this harmonic */ + int wmax; /* The maximum Fourier f-dot-dot for this harmonic */ + int numkern_zdim; /* Number of kernels calculated in the z dimension */ + int numkern_wdim; /* Number of kernels calculated in the w dimension */ + int numkern; /* Total number of kernels in the vector */ + kernel **kern; /* A 2D array of the kernels themselves, with dimensions of z and w */ + unsigned short *rinds; /* Table of lookup indices for Fourier Freqs: subharmonic r values corresponding to "fundamental" r values */ + unsigned short *zinds; /* Table of lookup indices for Fourier F-dots */ } subharminfo; typedef struct ffdotpows{ - int numrs; /* Number of Fourier freqs present */ - int numzs; /* Number of Fourier f-dots present */ - int rlo; /* Lowest Fourier freq present */ - int zlo; /* Lowest Fourier f-dot present */ - float **powers; /* Matrix of the powers */ - unsigned short *rinds; /* Table of indices for Fourier Freqs */ + long long rlo; /* Lowest Fourier freq present */ + int zlo; /* Lowest Fourier f-dot present */ + int wlo; /* Lowest Fourier f-dot-dot present */ + int numrs; /* Number of Fourier freqs present */ + int numzs; /* Number of Fourier f-dots present */ + int numws; /* Number of Fourier f-dot-dots present */ + float ***powers; /* 3D Matrix of the powers */ + unsigned short *rinds; /* Table of lookup indices for Fourier Freqs */ + unsigned short *zinds; /* Table of lookup indices for Fourier f-dots */ } ffdotpows; /* accel_utils.c */ -/* accel_utils.c */ - subharminfo **create_subharminfos(accelobs *obs); void free_subharminfos(accelobs *obs, subharminfo **shis); void create_accelobs(accelobs *obs, infodata *idata, @@ -127,7 +145,7 @@ void output_harmonics(GSList *list, accelobs *obs, infodata *idata); void free_accelcand(gpointer data, gpointer user_data); void print_accelcand(gpointer data, gpointer user_data); fcomplex *get_fourier_amplitudes(long long lobin, int numbins, accelobs *obs); -ffdotpows *subharm_ffdot_plane(int numharm, int harmnum, +ffdotpows *subharm_fderivs_vol(int numharm, int harmnum, double fullrlo, double fullrhi, subharminfo *shi, accelobs *obs); ffdotpows *copy_ffdotpows(ffdotpows *orig); diff --git a/include/accelsearch_cmd.h b/include/accelsearch_cmd.h index 0897fc7a5..39ee9a964 100644 --- a/include/accelsearch_cmd.h +++ b/include/accelsearch_cmd.h @@ -25,6 +25,10 @@ typedef struct s_Cmdline { char zmaxP; int zmax; int zmaxC; + /***** -wmax: The max (+ and -) Fourier freq double derivs to search */ + char wmaxP; + int wmax; + int wmaxC; /***** -sigma: Cutoff sigma for choosing candidates */ char sigmaP; float sigma; diff --git a/include/presto.h b/include/presto.h index 7920d1dca..a8af5059d 100644 --- a/include/presto.h +++ b/include/presto.h @@ -429,7 +429,7 @@ float get_numphotons(FILE * file); /* Arguments: */ /* 'file' is a pointer to the file you want to access. */ -double get_localpower(fcomplex *data, int numdata, double r); +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. */ @@ -437,7 +437,7 @@ double get_localpower(fcomplex *data, int numdata, double r); /* 'r' is the Fourier frequency in data that we want to */ /* interpolate. */ -double get_localpower3d(fcomplex *data, int numdata, double r, \ +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. */ @@ -451,7 +451,7 @@ double get_localpower3d(fcomplex *data, int numdata, double r, \ /* 'w' is the Fourier Frequency 2nd derivative (change in the */ /* Fourier f-dot during the observation). */ -void get_derivs3d(fcomplex *data, int numdata, double r, \ +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, */ @@ -718,7 +718,7 @@ void get_rawbin_cand(char *filenm, int candnum, rawbincand * cand); /* read_fft.c: */ /* Functions for getting information from an FFT file */ -fcomplex *read_fcomplex_file(FILE *file, int firstpt, int numpts); +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. */ @@ -729,7 +729,7 @@ fcomplex *read_fcomplex_file(FILE *file, int firstpt, int numpts); /* 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, int firstpt, int numpts); +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. */ @@ -801,7 +801,15 @@ float percolate_bin(binaryprops * list, int nlist); /* From prep_corr.c */ -void spread_with_pad(fcomplex *data, int numdata, \ +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 */ @@ -1114,12 +1122,12 @@ void rzw_interp(fcomplex *data, int numdata, double r, double z, \ /* In maximize_r.c and maximize_rw.c */ -double max_r_arr(fcomplex *data, int numdata, double rin, +double max_r_arr(fcomplex *data, long numdata, double rin, double *rout, rderivs *derivs); /* Return the Fourier frequency that maximizes the power. */ -double max_rz_arr(fcomplex *data, int numdata, double rin, double zin, \ +double max_rz_arr(fcomplex *data, long numdata, double rin, double zin, \ double *rout, double *zout, rderivs * derivs); /* Return the Fourier frequency and Fourier f-dot that */ /* maximizes the power. */ @@ -1130,33 +1138,32 @@ double max_rz_file(FILE *fftfile, double rin, double zin, \ /* maximizes the power of the candidate in 'fftfile'. */ -void max_rz_arr_harmonics(fcomplex * data[], int num_harmonics, - int r_offset[], - int numdata, double rin, double zin, - double *rout, double *zout, rderivs derivs[], - double power[]); +void max_rz_arr_harmonics(fcomplex data[], long numdata, + int num_harmonics, + double rin, double zin, + double *rout, double *zout, + rderivs derivs[], double powers[]); /* Return the Fourier frequency and Fourier f-dot that */ /* maximizes the power. */ -void max_rz_file_harmonics(FILE * fftfile, int num_harmonics, - int lobin, - double rin, double zin, - double *rout, double *zout, rderivs derivs[], - double maxpow[]); -/* Return the Fourier frequency and Fourier f-dot that */ -/* maximizes the power of the candidate in 'fftfile'. */ - - -double max_rzw_arr(fcomplex *data, int numdata, double rin, double zin, \ - double win, double *rout, double *zout, \ - double *wout, rderivs * derivs); - /* Return the Fourier frequency, f-dot, and fdotdot that */ - /* maximizes the power. */ +void max_rzw_arr_harmonics(fcomplex data[], long numdata, + int num_harmonics, + double rin, double zin, double win, + double *rout, double *zout, double *wout, + rderivs derivs[], double powers[]); +/* Return the Fourier frequency, f-dot, and f-dotdot that */ +/* maximizes the *summed* power of the multi-harmonic candidate */ + +double max_rzw_arr(fcomplex *data, long numdata, double rin, double zin, \ + double win, double *rout, double *zout, \ + double *wout, rderivs * derivs); +/* Return the Fourier frequency, f-dot, and fdotdot that */ +/* maximizes the power. */ -double max_rz_file(FILE *fftfile, double rin, double zin, \ - double *rout, double *zout, rderivs * derivs); - /* Return the Fourier frequency and Fourier f-dot that */ - /* maximizes the power of the candidate in 'fftfile'. */ +double max_rzw_file(FILE * fftfile, double rin, double zin, double win, \ + double *rout, double *zout, double *wout, rderivs * derivs); +/* Return the Fourier frequency, f-dot, and fdotdot that */ +/* maximizes the power of the candidate in 'fftfile'. */ /* In fold.c */ diff --git a/python/presto_src/__init__.py b/python/presto_src/__init__.py index e953d3e58..5ff149abc 100644 --- a/python/presto_src/__init__.py +++ b/python/presto_src/__init__.py @@ -177,10 +177,7 @@ def maximize_r(data, r, norm = None): """ rd = rderivs() (rmax, maxpow) = max_r_arr(data, r, rd) - if not norm: - maxpow = maxpow / rd.locpow - else: - maxpow = maxpow / norm + maxpow = maxpow / rd.locpow if norm is None else maxpow / norm return [maxpow, rmax, rd] def maximize_rz(data, r, z, norm = None): @@ -193,12 +190,33 @@ def maximize_rz(data, r, z, norm = None): """ rd = rderivs() (rmax, zmax, maxpow) = max_rz_arr(data, r, z, rd) - if not norm: - maxpow = maxpow / rd.locpow - else: - maxpow = maxpow / norm + maxpow = maxpow / rd.locpow if norm is None else maxpow / norm return [maxpow, rmax, zmax, rd] +def maximize_rz_harmonics(data, r, z, numharm, norm = None): + """ + maximize_rz_harmonics(data, r, z, numharm, norm = None): + Optimize the detection of a signal at location 'r', 'z' in + the F-Fdot plane, including harmonic summing of the harmonics. + The routine returns a list containing the optimized values of + the maximum normalized power, rmax, zmax, and a list of + rderivs structures for the peak. + """ + rds = [rderivs() for ii in range(numharm)] + derivdata = np.zeros(7 * numharm, dtype=np.float64) + rmax, zmax = max_rz_arr_harmonics(data, r, z, derivdata) + maxpow = 0.0 + for ii in range(numharm): + rds[ii].pow = derivdata[ii*7+0] + rds[ii].phs = derivdata[ii*7+1] + rds[ii].dpow = derivdata[ii*7+2] + rds[ii].dphs = derivdata[ii*7+3] + rds[ii].d2pow = derivdata[ii*7+4] + rds[ii].d2phs = derivdata[ii*7+5] + rds[ii].locpow = derivdata[ii*7+6] + maxpow += rds[ii].pow / rds[ii].locpow if norm is None else rds[ii].pow / norm + return [maxpow, rmax, zmax, rds] + def maximize_rzw(data, r, z, w, norm = None): """ maximize_rzw(data, r, z, w, norm = None): @@ -209,12 +227,33 @@ def maximize_rzw(data, r, z, w, norm = None): """ rd = rderivs() (rmax, zmax, wmax, maxpow) = max_rzw_arr(data, r, z, w, rd) - if not norm: - maxpow = maxpow / rd.locpow - else: - maxpow = maxpow / norm + maxpow = maxpow / rd.locpow if norm is None else maxpow / norm return [maxpow, rmax, zmax, wmax, rd] +def maximize_rzw_harmonics(data, r, z, w, numharm, norm = None): + """ + maximize_rzw_harmonics(data, r, z, w, numharm, norm = None): + Optimize the detection of a signal at location 'r', 'z', 'w' in + the F-Fd-Fdd volume, including harmonic summing of the harmonics. + The routine returns a list containing the optimized values of + the maximum normalized power, rmax, zmax, wmax, and a list of + rderivs structures for the peak. + """ + rds = [rderivs() for ii in range(numharm)] + derivdata = np.zeros(7 * numharm, dtype=np.float64) + rmax, zmax, wmax = max_rzw_arr_harmonics(data, r, z, w, derivdata) + maxpow = 0.0 + for ii in range(numharm): + rds[ii].pow = derivdata[ii*7+0] + rds[ii].phs = derivdata[ii*7+1] + rds[ii].dpow = derivdata[ii*7+2] + rds[ii].dphs = derivdata[ii*7+3] + rds[ii].d2pow = derivdata[ii*7+4] + rds[ii].d2phs = derivdata[ii*7+5] + rds[ii].locpow = derivdata[ii*7+6] + maxpow += rds[ii].pow / rds[ii].locpow if norm is None else rds[ii].pow / norm + return [maxpow, rmax, zmax, wmax, rds] + def search_fft(data, numcands, norm='default'): """ search_fft(data, numcands): diff --git a/python/presto_src/presto_wrap.c b/python/presto_src/presto_wrap.c index 907054b90..6156e85ca 100644 --- a/python/presto_src/presto_wrap.c +++ b/python/presto_src/presto_wrap.c @@ -4219,6 +4219,55 @@ SWIG_AsVal_unsigned_SS_long (PyObject *obj, unsigned long *val) } + void wrap_max_rz_arr_harmonics(fcomplex *data, long numdata, + double rin, double zin, + double *derivdata, int len, + double *rout, double *zout){ + int ii, numharm = len / 7; + double *powers = gen_dvect(numharm); + rderivs *derivs = (rderivs *)malloc(sizeof(rderivs) * numharm); + + max_rz_arr_harmonics(data, numdata, numharm, rin, zin, rout, zout, derivs, powers); + vect_free(powers); + // Hack to effectively return a array of rderivs + for (ii = 0 ; ii < numharm ; ii++) { + derivdata[ii*7+0] = derivs[ii].pow; + derivdata[ii*7+1] = derivs[ii].phs; + derivdata[ii*7+2] = derivs[ii].dpow; + derivdata[ii*7+3] = derivs[ii].dphs; + derivdata[ii*7+4] = derivs[ii].d2pow; + derivdata[ii*7+5] = derivs[ii].d2phs; + derivdata[ii*7+6] = derivs[ii].locpow; + } + free(derivs); + } + + + void wrap_max_rzw_arr_harmonics(fcomplex *data, long numdata, + double rin, double zin, double win, + double *derivdata, int len, + double *rout, double *zout, double *wout){ + int ii, numharm = len / 7; + double *powers = gen_dvect(numharm); + rderivs *derivs = (rderivs *)malloc(sizeof(rderivs) * numharm); + + max_rzw_arr_harmonics(data, numdata, numharm, rin, zin, win, + rout, zout, wout, derivs, powers); + vect_free(powers); + // Hack to effectively return a array of rderivs + for (ii = 0 ; ii < numharm ; ii++) { + derivdata[ii*7+0] = derivs[ii].pow; + derivdata[ii*7+1] = derivs[ii].phs; + derivdata[ii*7+2] = derivs[ii].dpow; + derivdata[ii*7+3] = derivs[ii].dphs; + derivdata[ii*7+4] = derivs[ii].d2pow; + derivdata[ii*7+5] = derivs[ii].d2phs; + derivdata[ii*7+6] = derivs[ii].locpow; + } + free(derivs); + } + + void wrap_max_rzw_arr(fcomplex * data, long numdata, double rin, double zin, double win, rderivs * derivs, double *rout, double *zout, double *wout, double *powout){ double pow; @@ -7727,10 +7776,10 @@ SWIGINTERN PyObject *_wrap_fresnl(PyObject *SWIGUNUSEDPARM(self), PyObject *args SWIGINTERN PyObject *_wrap_rderivs_pow_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7741,11 +7790,11 @@ SWIGINTERN PyObject *_wrap_rderivs_pow_set(PyObject *SWIGUNUSEDPARM(self), PyObj SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_pow_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_pow_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_pow_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->pow = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7760,7 +7809,7 @@ SWIGINTERN PyObject *_wrap_rderivs_pow_get(PyObject *SWIGUNUSEDPARM(self), PyObj void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_pow_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7768,8 +7817,8 @@ SWIGINTERN PyObject *_wrap_rderivs_pow_get(PyObject *SWIGUNUSEDPARM(self), PyObj SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_pow_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->pow); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->pow); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7779,10 +7828,10 @@ SWIGINTERN PyObject *_wrap_rderivs_pow_get(PyObject *SWIGUNUSEDPARM(self), PyObj SWIGINTERN PyObject *_wrap_rderivs_phs_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7793,11 +7842,11 @@ SWIGINTERN PyObject *_wrap_rderivs_phs_set(PyObject *SWIGUNUSEDPARM(self), PyObj SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_phs_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_phs_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_phs_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->phs = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7812,7 +7861,7 @@ SWIGINTERN PyObject *_wrap_rderivs_phs_get(PyObject *SWIGUNUSEDPARM(self), PyObj void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_phs_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7820,8 +7869,8 @@ SWIGINTERN PyObject *_wrap_rderivs_phs_get(PyObject *SWIGUNUSEDPARM(self), PyObj SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_phs_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->phs); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->phs); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7831,10 +7880,10 @@ SWIGINTERN PyObject *_wrap_rderivs_phs_get(PyObject *SWIGUNUSEDPARM(self), PyObj SWIGINTERN PyObject *_wrap_rderivs_dpow_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7845,11 +7894,11 @@ SWIGINTERN PyObject *_wrap_rderivs_dpow_set(PyObject *SWIGUNUSEDPARM(self), PyOb SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_dpow_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_dpow_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_dpow_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->dpow = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7864,7 +7913,7 @@ SWIGINTERN PyObject *_wrap_rderivs_dpow_get(PyObject *SWIGUNUSEDPARM(self), PyOb void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_dpow_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7872,8 +7921,8 @@ SWIGINTERN PyObject *_wrap_rderivs_dpow_get(PyObject *SWIGUNUSEDPARM(self), PyOb SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_dpow_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->dpow); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->dpow); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7883,10 +7932,10 @@ SWIGINTERN PyObject *_wrap_rderivs_dpow_get(PyObject *SWIGUNUSEDPARM(self), PyOb SWIGINTERN PyObject *_wrap_rderivs_dphs_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7897,11 +7946,11 @@ SWIGINTERN PyObject *_wrap_rderivs_dphs_set(PyObject *SWIGUNUSEDPARM(self), PyOb SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_dphs_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_dphs_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_dphs_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->dphs = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7916,7 +7965,7 @@ SWIGINTERN PyObject *_wrap_rderivs_dphs_get(PyObject *SWIGUNUSEDPARM(self), PyOb void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_dphs_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7924,8 +7973,8 @@ SWIGINTERN PyObject *_wrap_rderivs_dphs_get(PyObject *SWIGUNUSEDPARM(self), PyOb SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_dphs_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->dphs); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->dphs); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7935,10 +7984,10 @@ SWIGINTERN PyObject *_wrap_rderivs_dphs_get(PyObject *SWIGUNUSEDPARM(self), PyOb SWIGINTERN PyObject *_wrap_rderivs_d2pow_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7949,11 +7998,11 @@ SWIGINTERN PyObject *_wrap_rderivs_d2pow_set(PyObject *SWIGUNUSEDPARM(self), PyO SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_d2pow_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_d2pow_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_d2pow_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->d2pow = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7968,7 +8017,7 @@ SWIGINTERN PyObject *_wrap_rderivs_d2pow_get(PyObject *SWIGUNUSEDPARM(self), PyO void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_d2pow_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7976,8 +8025,8 @@ SWIGINTERN PyObject *_wrap_rderivs_d2pow_get(PyObject *SWIGUNUSEDPARM(self), PyO SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_d2pow_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->d2pow); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->d2pow); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7987,10 +8036,10 @@ SWIGINTERN PyObject *_wrap_rderivs_d2pow_get(PyObject *SWIGUNUSEDPARM(self), PyO SWIGINTERN PyObject *_wrap_rderivs_d2phs_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -8001,11 +8050,11 @@ SWIGINTERN PyObject *_wrap_rderivs_d2phs_set(PyObject *SWIGUNUSEDPARM(self), PyO SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_d2phs_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_d2phs_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_d2phs_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->d2phs = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -8020,7 +8069,7 @@ SWIGINTERN PyObject *_wrap_rderivs_d2phs_get(PyObject *SWIGUNUSEDPARM(self), PyO void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_d2phs_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -8028,8 +8077,8 @@ SWIGINTERN PyObject *_wrap_rderivs_d2phs_get(PyObject *SWIGUNUSEDPARM(self), PyO SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_d2phs_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->d2phs); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->d2phs); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -8039,10 +8088,10 @@ SWIGINTERN PyObject *_wrap_rderivs_d2phs_get(PyObject *SWIGUNUSEDPARM(self), PyO SWIGINTERN PyObject *_wrap_rderivs_locpow_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -8053,11 +8102,11 @@ SWIGINTERN PyObject *_wrap_rderivs_locpow_set(PyObject *SWIGUNUSEDPARM(self), Py SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_locpow_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_locpow_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_locpow_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->locpow = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -8072,7 +8121,7 @@ SWIGINTERN PyObject *_wrap_rderivs_locpow_get(PyObject *SWIGUNUSEDPARM(self), Py void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_locpow_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -8080,8 +8129,8 @@ SWIGINTERN PyObject *_wrap_rderivs_locpow_get(PyObject *SWIGUNUSEDPARM(self), Py SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_locpow_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->locpow); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->locpow); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -11140,31 +11189,30 @@ SWIGINTERN PyObject *_wrap_gen_bin_response(PyObject *SWIGUNUSEDPARM(self), PyOb SWIGINTERN PyObject *_wrap_get_localpower(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; fcomplex *arg1 = (fcomplex *) 0 ; - int arg2 ; + long arg2 ; double arg3 ; - void *argp1 = 0 ; - int res1 = 0 ; - int val2 ; - int ecode2 = 0 ; + PyArrayObject *array1 = NULL ; + int is_new_object1 = 0 ; double val3 ; int ecode3 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; - PyObject * obj2 = 0 ; float result; - if (!PyArg_ParseTuple(args,(char *)"OOO:get_localpower",&obj0,&obj1,&obj2)) SWIG_fail; - res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_FCOMPLEX, 0 | 0 ); - if (!SWIG_IsOK(res1)) { - SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "get_localpower" "', argument " "1"" of type '" "fcomplex *""'"); + if (!PyArg_ParseTuple(args,(char *)"OO:get_localpower",&obj0,&obj1)) SWIG_fail; + { + npy_intp size[1] = { + -1 + }; + array1 = obj_to_array_contiguous_allow_conversion(obj0, + NPY_CFLOAT, + &is_new_object1); + if (!array1 || !require_dimensions(array1, 1) || + !require_size(array1, size, 1)) SWIG_fail; + arg1 = (fcomplex*) array_data(array1); + arg2 = (long) array_size(array1,0); } - arg1 = (fcomplex *)(argp1); - ecode2 = SWIG_AsVal_int(obj1, &val2); - if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "get_localpower" "', argument " "2"" of type '" "int""'"); - } - arg2 = (int)(val2); - ecode3 = SWIG_AsVal_double(obj2, &val3); + ecode3 = SWIG_AsVal_double(obj1, &val3); if (!SWIG_IsOK(ecode3)) { SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "get_localpower" "', argument " "3"" of type '" "double""'"); } @@ -11187,8 +11235,20 @@ SWIGINTERN PyObject *_wrap_get_localpower(PyObject *SWIGUNUSEDPARM(self), PyObje } } resultobj = SWIG_From_float((float)(result)); + { + if (is_new_object1 && array1) + { + Py_DECREF(array1); + } + } return resultobj; fail: + { + if (is_new_object1 && array1) + { + Py_DECREF(array1); + } + } return NULL; } @@ -11196,14 +11256,12 @@ SWIGINTERN PyObject *_wrap_get_localpower(PyObject *SWIGUNUSEDPARM(self), PyObje SWIGINTERN PyObject *_wrap_get_localpower3d(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; fcomplex *arg1 = (fcomplex *) 0 ; - int arg2 ; + long arg2 ; double arg3 ; double arg4 ; double arg5 ; - void *argp1 = 0 ; - int res1 = 0 ; - int val2 ; - int ecode2 = 0 ; + PyArrayObject *array1 = NULL ; + int is_new_object1 = 0 ; double val3 ; int ecode3 = 0 ; double val4 ; @@ -11214,31 +11272,32 @@ SWIGINTERN PyObject *_wrap_get_localpower3d(PyObject *SWIGUNUSEDPARM(self), PyOb PyObject * obj1 = 0 ; PyObject * obj2 = 0 ; PyObject * obj3 = 0 ; - PyObject * obj4 = 0 ; float result; - if (!PyArg_ParseTuple(args,(char *)"OOOOO:get_localpower3d",&obj0,&obj1,&obj2,&obj3,&obj4)) SWIG_fail; - res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_FCOMPLEX, 0 | 0 ); - if (!SWIG_IsOK(res1)) { - SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "get_localpower3d" "', argument " "1"" of type '" "fcomplex *""'"); + if (!PyArg_ParseTuple(args,(char *)"OOOO:get_localpower3d",&obj0,&obj1,&obj2,&obj3)) SWIG_fail; + { + npy_intp size[1] = { + -1 + }; + array1 = obj_to_array_contiguous_allow_conversion(obj0, + NPY_CFLOAT, + &is_new_object1); + if (!array1 || !require_dimensions(array1, 1) || + !require_size(array1, size, 1)) SWIG_fail; + arg1 = (fcomplex*) array_data(array1); + arg2 = (long) array_size(array1,0); } - arg1 = (fcomplex *)(argp1); - ecode2 = SWIG_AsVal_int(obj1, &val2); - if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "get_localpower3d" "', argument " "2"" of type '" "int""'"); - } - arg2 = (int)(val2); - ecode3 = SWIG_AsVal_double(obj2, &val3); + ecode3 = SWIG_AsVal_double(obj1, &val3); if (!SWIG_IsOK(ecode3)) { SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "get_localpower3d" "', argument " "3"" of type '" "double""'"); } arg3 = (double)(val3); - ecode4 = SWIG_AsVal_double(obj3, &val4); + ecode4 = SWIG_AsVal_double(obj2, &val4); if (!SWIG_IsOK(ecode4)) { SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "get_localpower3d" "', argument " "4"" of type '" "double""'"); } arg4 = (double)(val4); - ecode5 = SWIG_AsVal_double(obj4, &val5); + ecode5 = SWIG_AsVal_double(obj3, &val5); if (!SWIG_IsOK(ecode5)) { SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "get_localpower3d" "', argument " "5"" of type '" "double""'"); } @@ -11261,8 +11320,20 @@ SWIGINTERN PyObject *_wrap_get_localpower3d(PyObject *SWIGUNUSEDPARM(self), PyOb } } resultobj = SWIG_From_float((float)(result)); + { + if (is_new_object1 && array1) + { + Py_DECREF(array1); + } + } return resultobj; fail: + { + if (is_new_object1 && array1) + { + Py_DECREF(array1); + } + } return NULL; } @@ -11270,7 +11341,7 @@ SWIGINTERN PyObject *_wrap_get_localpower3d(PyObject *SWIGUNUSEDPARM(self), PyOb SWIGINTERN PyObject *_wrap_get_derivs3d(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; fcomplex *arg1 = (fcomplex *) 0 ; - int arg2 ; + long arg2 ; double arg3 ; double arg4 ; double arg5 ; @@ -11278,7 +11349,7 @@ SWIGINTERN PyObject *_wrap_get_derivs3d(PyObject *SWIGUNUSEDPARM(self), PyObject rderivs *arg7 = (rderivs *) 0 ; void *argp1 = 0 ; int res1 = 0 ; - int val2 ; + long val2 ; int ecode2 = 0 ; double val3 ; int ecode3 = 0 ; @@ -11304,11 +11375,11 @@ SWIGINTERN PyObject *_wrap_get_derivs3d(PyObject *SWIGUNUSEDPARM(self), PyObject SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "get_derivs3d" "', argument " "1"" of type '" "fcomplex *""'"); } arg1 = (fcomplex *)(argp1); - ecode2 = SWIG_AsVal_int(obj1, &val2); + ecode2 = SWIG_AsVal_long(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "get_derivs3d" "', argument " "2"" of type '" "int""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "get_derivs3d" "', argument " "2"" of type '" "long""'"); } - arg2 = (int)(val2); + arg2 = (long)(val2); ecode3 = SWIG_AsVal_double(obj2, &val3); if (!SWIG_IsOK(ecode3)) { SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "get_derivs3d" "', argument " "3"" of type '" "double""'"); @@ -13216,6 +13287,209 @@ SWIGINTERN PyObject *_wrap_max_rz_arr(PyObject *SWIGUNUSEDPARM(self), PyObject * } +SWIGINTERN PyObject *_wrap_max_rz_arr_harmonics(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { + PyObject *resultobj = 0; + fcomplex *arg1 = (fcomplex *) 0 ; + long arg2 ; + double arg3 ; + double arg4 ; + double *arg5 = (double *) 0 ; + int arg6 ; + double *arg7 = (double *) 0 ; + double *arg8 = (double *) 0 ; + PyArrayObject *array1 = NULL ; + int i1 = 1 ; + double val3 ; + int ecode3 = 0 ; + double val4 ; + int ecode4 = 0 ; + PyArrayObject *array5 = NULL ; + int i5 = 1 ; + double temp7 ; + int res7 = SWIG_TMPOBJ ; + double temp8 ; + int res8 = SWIG_TMPOBJ ; + PyObject * obj0 = 0 ; + PyObject * obj1 = 0 ; + PyObject * obj2 = 0 ; + PyObject * obj3 = 0 ; + + arg7 = &temp7; + arg8 = &temp8; + if (!PyArg_ParseTuple(args,(char *)"OOOO:max_rz_arr_harmonics",&obj0,&obj1,&obj2,&obj3)) SWIG_fail; + { + array1 = obj_to_array_no_conversion(obj0, NPY_CFLOAT); + if (!array1 || !require_dimensions(array1,1) || !require_contiguous(array1) + || !require_native(array1)) SWIG_fail; + arg1 = (fcomplex*) array_data(array1); + arg2 = 1; + for (i1=0; i1 < array_numdims(array1); ++i1) arg2 *= array_size(array1,i1); + } + ecode3 = SWIG_AsVal_double(obj1, &val3); + if (!SWIG_IsOK(ecode3)) { + SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "max_rz_arr_harmonics" "', argument " "3"" of type '" "double""'"); + } + arg3 = (double)(val3); + ecode4 = SWIG_AsVal_double(obj2, &val4); + if (!SWIG_IsOK(ecode4)) { + SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "max_rz_arr_harmonics" "', argument " "4"" of type '" "double""'"); + } + arg4 = (double)(val4); + { + array5 = obj_to_array_no_conversion(obj3, NPY_DOUBLE); + if (!array5 || !require_dimensions(array5,1) || !require_contiguous(array5) + || !require_native(array5)) SWIG_fail; + arg5 = (double*) array_data(array5); + arg6 = 1; + for (i5=0; i5 < array_numdims(array5); ++i5) arg6 *= array_size(array5,i5); + } + { + errno = 0; + wrap_max_rz_arr_harmonics(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8); + + if (errno != 0) + { + switch(errno) + { + case ENOMEM: + PyErr_Format(PyExc_MemoryError, "Failed malloc()"); + break; + default: + PyErr_Format(PyExc_Exception, "Unknown exception"); + } + SWIG_fail; + } + } + resultobj = SWIG_Py_Void(); + if (SWIG_IsTmpObj(res7)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg7))); + } else { + int new_flags = SWIG_IsNewObj(res7) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg7), SWIGTYPE_p_double, new_flags)); + } + if (SWIG_IsTmpObj(res8)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg8))); + } else { + int new_flags = SWIG_IsNewObj(res8) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg8), SWIGTYPE_p_double, new_flags)); + } + return resultobj; +fail: + return NULL; +} + + +SWIGINTERN PyObject *_wrap_max_rzw_arr_harmonics(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { + PyObject *resultobj = 0; + fcomplex *arg1 = (fcomplex *) 0 ; + long arg2 ; + double arg3 ; + double arg4 ; + double arg5 ; + double *arg6 = (double *) 0 ; + int arg7 ; + double *arg8 = (double *) 0 ; + double *arg9 = (double *) 0 ; + double *arg10 = (double *) 0 ; + PyArrayObject *array1 = NULL ; + int i1 = 1 ; + double val3 ; + int ecode3 = 0 ; + double val4 ; + int ecode4 = 0 ; + double val5 ; + int ecode5 = 0 ; + PyArrayObject *array6 = NULL ; + int i6 = 1 ; + double temp8 ; + int res8 = SWIG_TMPOBJ ; + double temp9 ; + int res9 = SWIG_TMPOBJ ; + double temp10 ; + int res10 = SWIG_TMPOBJ ; + PyObject * obj0 = 0 ; + PyObject * obj1 = 0 ; + PyObject * obj2 = 0 ; + PyObject * obj3 = 0 ; + PyObject * obj4 = 0 ; + + arg8 = &temp8; + arg9 = &temp9; + arg10 = &temp10; + if (!PyArg_ParseTuple(args,(char *)"OOOOO:max_rzw_arr_harmonics",&obj0,&obj1,&obj2,&obj3,&obj4)) SWIG_fail; + { + array1 = obj_to_array_no_conversion(obj0, NPY_CFLOAT); + if (!array1 || !require_dimensions(array1,1) || !require_contiguous(array1) + || !require_native(array1)) SWIG_fail; + arg1 = (fcomplex*) array_data(array1); + arg2 = 1; + for (i1=0; i1 < array_numdims(array1); ++i1) arg2 *= array_size(array1,i1); + } + ecode3 = SWIG_AsVal_double(obj1, &val3); + if (!SWIG_IsOK(ecode3)) { + SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "max_rzw_arr_harmonics" "', argument " "3"" of type '" "double""'"); + } + arg3 = (double)(val3); + ecode4 = SWIG_AsVal_double(obj2, &val4); + if (!SWIG_IsOK(ecode4)) { + SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "max_rzw_arr_harmonics" "', argument " "4"" of type '" "double""'"); + } + arg4 = (double)(val4); + ecode5 = SWIG_AsVal_double(obj3, &val5); + if (!SWIG_IsOK(ecode5)) { + SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "max_rzw_arr_harmonics" "', argument " "5"" of type '" "double""'"); + } + arg5 = (double)(val5); + { + array6 = obj_to_array_no_conversion(obj4, NPY_DOUBLE); + if (!array6 || !require_dimensions(array6,1) || !require_contiguous(array6) + || !require_native(array6)) SWIG_fail; + arg6 = (double*) array_data(array6); + arg7 = 1; + for (i6=0; i6 < array_numdims(array6); ++i6) arg7 *= array_size(array6,i6); + } + { + errno = 0; + wrap_max_rzw_arr_harmonics(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9,arg10); + + if (errno != 0) + { + switch(errno) + { + case ENOMEM: + PyErr_Format(PyExc_MemoryError, "Failed malloc()"); + break; + default: + PyErr_Format(PyExc_Exception, "Unknown exception"); + } + SWIG_fail; + } + } + resultobj = SWIG_Py_Void(); + if (SWIG_IsTmpObj(res8)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg8))); + } else { + int new_flags = SWIG_IsNewObj(res8) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg8), SWIGTYPE_p_double, new_flags)); + } + if (SWIG_IsTmpObj(res9)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg9))); + } else { + int new_flags = SWIG_IsNewObj(res9) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg9), SWIGTYPE_p_double, new_flags)); + } + if (SWIG_IsTmpObj(res10)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg10))); + } else { + int new_flags = SWIG_IsNewObj(res10) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg10), SWIGTYPE_p_double, new_flags)); + } + return resultobj; +fail: + return NULL; +} + + SWIGINTERN PyObject *_wrap_max_rzw_arr(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; fcomplex *arg1 = (fcomplex *) 0 ; @@ -13706,6 +13980,8 @@ static PyMethodDef SwigMethods[] = { { (char *)"corr_rzw_vol", _wrap_corr_rzw_vol, METH_VARARGS, NULL}, { (char *)"max_r_arr", _wrap_max_r_arr, METH_VARARGS, NULL}, { (char *)"max_rz_arr", _wrap_max_rz_arr, METH_VARARGS, NULL}, + { (char *)"max_rz_arr_harmonics", _wrap_max_rz_arr_harmonics, METH_VARARGS, NULL}, + { (char *)"max_rzw_arr_harmonics", _wrap_max_rzw_arr_harmonics, METH_VARARGS, NULL}, { (char *)"max_rzw_arr", _wrap_max_rzw_arr, METH_VARARGS, NULL}, { (char *)"barycenter", _wrap_barycenter, METH_VARARGS, NULL}, { NULL, NULL, 0, NULL } diff --git a/python/presto_src/prestoswig.py b/python/presto_src/prestoswig.py index d07144f5b..c7f027aa2 100644 --- a/python/presto_src/prestoswig.py +++ b/python/presto_src/prestoswig.py @@ -730,12 +730,12 @@ def gen_bin_response(roffset, numbetween, numkern, ppsr, T, orbit): return _presto.gen_bin_response(roffset, numbetween, numkern, ppsr, T, orbit) gen_bin_response = _presto.gen_bin_response -def get_localpower(data, numdata, r): - return _presto.get_localpower(data, numdata, r) +def get_localpower(data, r): + return _presto.get_localpower(data, r) get_localpower = _presto.get_localpower -def get_localpower3d(data, numdata, r, z, w): - return _presto.get_localpower3d(data, numdata, r, z, w) +def get_localpower3d(data, r, z, w): + return _presto.get_localpower3d(data, r, z, w) get_localpower3d = _presto.get_localpower3d def get_derivs3d(data, numdata, r, z, w, localpower, result): @@ -854,6 +854,14 @@ def max_rz_arr(data, rin, zin, derivs): return _presto.max_rz_arr(data, rin, zin, derivs) max_rz_arr = _presto.max_rz_arr +def max_rz_arr_harmonics(data, rin, zin, derivdata): + return _presto.max_rz_arr_harmonics(data, rin, zin, derivdata) +max_rz_arr_harmonics = _presto.max_rz_arr_harmonics + +def max_rzw_arr_harmonics(data, rin, zin, win, derivdata): + return _presto.max_rzw_arr_harmonics(data, rin, zin, win, derivdata) +max_rzw_arr_harmonics = _presto.max_rzw_arr_harmonics + def max_rzw_arr(data, rin, zin, win, derivs): return _presto.max_rzw_arr(data, rin, zin, win, derivs) max_rzw_arr = _presto.max_rzw_arr diff --git a/python/wrappers/presto.i b/python/wrappers/presto.i index 211bd7a21..2263d3126 100644 --- a/python/wrappers/presto.i +++ b/python/wrappers/presto.i @@ -255,13 +255,13 @@ int fresnl(double xxa, double *ssa, double *cca); // Return the Fresnel inegrals typedef struct RDERIVS { - float pow; /* Power normalized with local power */ - float phs; /* Signal phase */ - float dpow; /* 1st deriv of power wrt fourier freq */ - float dphs; /* 1st deriv of phase wrt fourier freq */ - float d2pow; /* 2nd deriv of power wrt fourier freq */ - float d2phs; /* 2nd deriv of power wrt fourier freq */ - float locpow; /* Local mean power level */ + 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 { @@ -554,8 +554,8 @@ void wrap_gen_bin_response(double roffset, int numbetween, int numkern, %} %clear (fcomplex **vect, long *nn); -%apply (fcomplex* IN_ARRAY1, int DIM1) {(fcomplex *data, int numdata)}; -float get_localpower(fcomplex *data, int numdata, double r); +%apply (fcomplex* IN_ARRAY1, long DIM1) {(fcomplex *data, long numdata)}; +float 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. */ @@ -563,7 +563,7 @@ float get_localpower(fcomplex *data, int numdata, double r); /* 'r' is the Fourier frequency in data that we want to */ /* interpolate. */ -float get_localpower3d(fcomplex *data, int numdata, double r, +float 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. */ @@ -576,8 +576,9 @@ float get_localpower3d(fcomplex *data, int numdata, double r, /* signal smears over during the observation). */ /* 'w' is the Fourier Frequency 2nd derivative (change in the */ /* Fourier f-dot during the observation). */ +%clear (fcomplex *data, long numdata); -void get_derivs3d(fcomplex *data, int numdata, double r, +void get_derivs3d(fcomplex *data, long numdata, double r, double z, double w, float localpower, rderivs *result); /* Return an rderives structure that contains the power, */ @@ -625,7 +626,6 @@ void calc_binprops(fourierprops * props, double T, int lowbin, /* 'absnorm' is the value of the power normalization */ /* constant for this mini-FFT. */ /* 'result' is the returned binaryprops structure. */ -%clear (fcomplex *data, int numdata); void calc_rzwerrs(fourierprops *props, double T, rzwerrs *result); /* Calculate periods, frequencies, their derivatives */ @@ -778,7 +778,6 @@ double sphere_ang_diff(double ra1, double dec1, double ra2, double dec2); %} %apply double *OUTPUT { double *rout, double *zout, double *powout }; -%apply (fcomplex* INPLACE_ARRAY1, long DIM1) {(fcomplex *data, long numdata)}; %rename (max_rz_arr) wrap_max_rz_arr; %inline %{ void wrap_max_rz_arr(fcomplex * data, long numdata, double rin, double zin, @@ -789,8 +788,64 @@ double sphere_ang_diff(double ra1, double dec1, double ra2, double dec2); } %} +%apply double *OUTPUT { double *rout, double *zout }; +%apply (double* INPLACE_ARRAY1, int DIM1) {(double *derivdata, int len)}; +%rename (max_rz_arr_harmonics) wrap_max_rz_arr_harmonics; +%inline %{ + void wrap_max_rz_arr_harmonics(fcomplex *data, long numdata, + double rin, double zin, + double *derivdata, int len, + double *rout, double *zout){ + int ii, numharm = len / 7; + double *powers = gen_dvect(numharm); + rderivs *derivs = (rderivs *)malloc(sizeof(rderivs) * numharm); + + max_rz_arr_harmonics(data, numdata, numharm, rin, zin, rout, zout, derivs, powers); + vect_free(powers); + // Hack to effectively return a array of rderivs + for (ii = 0 ; ii < numharm ; ii++) { + derivdata[ii*7+0] = derivs[ii].pow; + derivdata[ii*7+1] = derivs[ii].phs; + derivdata[ii*7+2] = derivs[ii].dpow; + derivdata[ii*7+3] = derivs[ii].dphs; + derivdata[ii*7+4] = derivs[ii].d2pow; + derivdata[ii*7+5] = derivs[ii].d2phs; + derivdata[ii*7+6] = derivs[ii].locpow; + } + free(derivs); + } +%} + +%apply double *OUTPUT { double *rout, double *zout, double *wout }; +%rename (max_rzw_arr_harmonics) wrap_max_rzw_arr_harmonics; +%inline %{ + void wrap_max_rzw_arr_harmonics(fcomplex *data, long numdata, + double rin, double zin, double win, + double *derivdata, int len, + double *rout, double *zout, double *wout){ + int ii, numharm = len / 7; + double *powers = gen_dvect(numharm); + rderivs *derivs = (rderivs *)malloc(sizeof(rderivs) * numharm); + + max_rzw_arr_harmonics(data, numdata, numharm, rin, zin, win, + rout, zout, wout, derivs, powers); + vect_free(powers); + // Hack to effectively return a array of rderivs + for (ii = 0 ; ii < numharm ; ii++) { + derivdata[ii*7+0] = derivs[ii].pow; + derivdata[ii*7+1] = derivs[ii].phs; + derivdata[ii*7+2] = derivs[ii].dpow; + derivdata[ii*7+3] = derivs[ii].dphs; + derivdata[ii*7+4] = derivs[ii].d2pow; + derivdata[ii*7+5] = derivs[ii].d2phs; + derivdata[ii*7+6] = derivs[ii].locpow; + } + free(derivs); + } +%} +%clear (double *derivdata, int len); + %apply double *OUTPUT { double *rout, double *zout, double *wout, double *powout }; -%apply (fcomplex* INPLACE_ARRAY1, long DIM1) {(fcomplex *data, long numdata)}; %rename (max_rzw_arr) wrap_max_rzw_arr; %inline %{ void wrap_max_rzw_arr(fcomplex * data, long numdata, double rin, double zin, double win, @@ -800,6 +855,7 @@ double sphere_ang_diff(double ra1, double dec1, double ra2, double dec2); *powout = pow; } %} +%clear (fcomplex *data, long numdata); %apply (double* INPLACE_ARRAY1, long DIM1) {(double *topotimes, long N1)}; %apply (double* INPLACE_ARRAY1, long DIM1) {(double *barytimes, long N2)}; diff --git a/python/wrappers/presto.py b/python/wrappers/presto.py index d07144f5b..c7f027aa2 100644 --- a/python/wrappers/presto.py +++ b/python/wrappers/presto.py @@ -730,12 +730,12 @@ def gen_bin_response(roffset, numbetween, numkern, ppsr, T, orbit): return _presto.gen_bin_response(roffset, numbetween, numkern, ppsr, T, orbit) gen_bin_response = _presto.gen_bin_response -def get_localpower(data, numdata, r): - return _presto.get_localpower(data, numdata, r) +def get_localpower(data, r): + return _presto.get_localpower(data, r) get_localpower = _presto.get_localpower -def get_localpower3d(data, numdata, r, z, w): - return _presto.get_localpower3d(data, numdata, r, z, w) +def get_localpower3d(data, r, z, w): + return _presto.get_localpower3d(data, r, z, w) get_localpower3d = _presto.get_localpower3d def get_derivs3d(data, numdata, r, z, w, localpower, result): @@ -854,6 +854,14 @@ def max_rz_arr(data, rin, zin, derivs): return _presto.max_rz_arr(data, rin, zin, derivs) max_rz_arr = _presto.max_rz_arr +def max_rz_arr_harmonics(data, rin, zin, derivdata): + return _presto.max_rz_arr_harmonics(data, rin, zin, derivdata) +max_rz_arr_harmonics = _presto.max_rz_arr_harmonics + +def max_rzw_arr_harmonics(data, rin, zin, win, derivdata): + return _presto.max_rzw_arr_harmonics(data, rin, zin, win, derivdata) +max_rzw_arr_harmonics = _presto.max_rzw_arr_harmonics + def max_rzw_arr(data, rin, zin, win, derivs): return _presto.max_rzw_arr(data, rin, zin, win, derivs) max_rzw_arr = _presto.max_rzw_arr diff --git a/python/wrappers/presto_wrap.c b/python/wrappers/presto_wrap.c index 907054b90..6156e85ca 100644 --- a/python/wrappers/presto_wrap.c +++ b/python/wrappers/presto_wrap.c @@ -4219,6 +4219,55 @@ SWIG_AsVal_unsigned_SS_long (PyObject *obj, unsigned long *val) } + void wrap_max_rz_arr_harmonics(fcomplex *data, long numdata, + double rin, double zin, + double *derivdata, int len, + double *rout, double *zout){ + int ii, numharm = len / 7; + double *powers = gen_dvect(numharm); + rderivs *derivs = (rderivs *)malloc(sizeof(rderivs) * numharm); + + max_rz_arr_harmonics(data, numdata, numharm, rin, zin, rout, zout, derivs, powers); + vect_free(powers); + // Hack to effectively return a array of rderivs + for (ii = 0 ; ii < numharm ; ii++) { + derivdata[ii*7+0] = derivs[ii].pow; + derivdata[ii*7+1] = derivs[ii].phs; + derivdata[ii*7+2] = derivs[ii].dpow; + derivdata[ii*7+3] = derivs[ii].dphs; + derivdata[ii*7+4] = derivs[ii].d2pow; + derivdata[ii*7+5] = derivs[ii].d2phs; + derivdata[ii*7+6] = derivs[ii].locpow; + } + free(derivs); + } + + + void wrap_max_rzw_arr_harmonics(fcomplex *data, long numdata, + double rin, double zin, double win, + double *derivdata, int len, + double *rout, double *zout, double *wout){ + int ii, numharm = len / 7; + double *powers = gen_dvect(numharm); + rderivs *derivs = (rderivs *)malloc(sizeof(rderivs) * numharm); + + max_rzw_arr_harmonics(data, numdata, numharm, rin, zin, win, + rout, zout, wout, derivs, powers); + vect_free(powers); + // Hack to effectively return a array of rderivs + for (ii = 0 ; ii < numharm ; ii++) { + derivdata[ii*7+0] = derivs[ii].pow; + derivdata[ii*7+1] = derivs[ii].phs; + derivdata[ii*7+2] = derivs[ii].dpow; + derivdata[ii*7+3] = derivs[ii].dphs; + derivdata[ii*7+4] = derivs[ii].d2pow; + derivdata[ii*7+5] = derivs[ii].d2phs; + derivdata[ii*7+6] = derivs[ii].locpow; + } + free(derivs); + } + + void wrap_max_rzw_arr(fcomplex * data, long numdata, double rin, double zin, double win, rderivs * derivs, double *rout, double *zout, double *wout, double *powout){ double pow; @@ -7727,10 +7776,10 @@ SWIGINTERN PyObject *_wrap_fresnl(PyObject *SWIGUNUSEDPARM(self), PyObject *args SWIGINTERN PyObject *_wrap_rderivs_pow_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7741,11 +7790,11 @@ SWIGINTERN PyObject *_wrap_rderivs_pow_set(PyObject *SWIGUNUSEDPARM(self), PyObj SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_pow_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_pow_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_pow_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->pow = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7760,7 +7809,7 @@ SWIGINTERN PyObject *_wrap_rderivs_pow_get(PyObject *SWIGUNUSEDPARM(self), PyObj void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_pow_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7768,8 +7817,8 @@ SWIGINTERN PyObject *_wrap_rderivs_pow_get(PyObject *SWIGUNUSEDPARM(self), PyObj SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_pow_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->pow); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->pow); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7779,10 +7828,10 @@ SWIGINTERN PyObject *_wrap_rderivs_pow_get(PyObject *SWIGUNUSEDPARM(self), PyObj SWIGINTERN PyObject *_wrap_rderivs_phs_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7793,11 +7842,11 @@ SWIGINTERN PyObject *_wrap_rderivs_phs_set(PyObject *SWIGUNUSEDPARM(self), PyObj SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_phs_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_phs_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_phs_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->phs = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7812,7 +7861,7 @@ SWIGINTERN PyObject *_wrap_rderivs_phs_get(PyObject *SWIGUNUSEDPARM(self), PyObj void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_phs_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7820,8 +7869,8 @@ SWIGINTERN PyObject *_wrap_rderivs_phs_get(PyObject *SWIGUNUSEDPARM(self), PyObj SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_phs_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->phs); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->phs); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7831,10 +7880,10 @@ SWIGINTERN PyObject *_wrap_rderivs_phs_get(PyObject *SWIGUNUSEDPARM(self), PyObj SWIGINTERN PyObject *_wrap_rderivs_dpow_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7845,11 +7894,11 @@ SWIGINTERN PyObject *_wrap_rderivs_dpow_set(PyObject *SWIGUNUSEDPARM(self), PyOb SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_dpow_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_dpow_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_dpow_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->dpow = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7864,7 +7913,7 @@ SWIGINTERN PyObject *_wrap_rderivs_dpow_get(PyObject *SWIGUNUSEDPARM(self), PyOb void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_dpow_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7872,8 +7921,8 @@ SWIGINTERN PyObject *_wrap_rderivs_dpow_get(PyObject *SWIGUNUSEDPARM(self), PyOb SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_dpow_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->dpow); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->dpow); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7883,10 +7932,10 @@ SWIGINTERN PyObject *_wrap_rderivs_dpow_get(PyObject *SWIGUNUSEDPARM(self), PyOb SWIGINTERN PyObject *_wrap_rderivs_dphs_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7897,11 +7946,11 @@ SWIGINTERN PyObject *_wrap_rderivs_dphs_set(PyObject *SWIGUNUSEDPARM(self), PyOb SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_dphs_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_dphs_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_dphs_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->dphs = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7916,7 +7965,7 @@ SWIGINTERN PyObject *_wrap_rderivs_dphs_get(PyObject *SWIGUNUSEDPARM(self), PyOb void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_dphs_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7924,8 +7973,8 @@ SWIGINTERN PyObject *_wrap_rderivs_dphs_get(PyObject *SWIGUNUSEDPARM(self), PyOb SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_dphs_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->dphs); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->dphs); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7935,10 +7984,10 @@ SWIGINTERN PyObject *_wrap_rderivs_dphs_get(PyObject *SWIGUNUSEDPARM(self), PyOb SWIGINTERN PyObject *_wrap_rderivs_d2pow_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -7949,11 +7998,11 @@ SWIGINTERN PyObject *_wrap_rderivs_d2pow_set(PyObject *SWIGUNUSEDPARM(self), PyO SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_d2pow_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_d2pow_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_d2pow_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->d2pow = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -7968,7 +8017,7 @@ SWIGINTERN PyObject *_wrap_rderivs_d2pow_get(PyObject *SWIGUNUSEDPARM(self), PyO void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_d2pow_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -7976,8 +8025,8 @@ SWIGINTERN PyObject *_wrap_rderivs_d2pow_get(PyObject *SWIGUNUSEDPARM(self), PyO SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_d2pow_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->d2pow); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->d2pow); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -7987,10 +8036,10 @@ SWIGINTERN PyObject *_wrap_rderivs_d2pow_get(PyObject *SWIGUNUSEDPARM(self), PyO SWIGINTERN PyObject *_wrap_rderivs_d2phs_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -8001,11 +8050,11 @@ SWIGINTERN PyObject *_wrap_rderivs_d2phs_set(PyObject *SWIGUNUSEDPARM(self), PyO SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_d2phs_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_d2phs_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_d2phs_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->d2phs = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -8020,7 +8069,7 @@ SWIGINTERN PyObject *_wrap_rderivs_d2phs_get(PyObject *SWIGUNUSEDPARM(self), PyO void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_d2phs_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -8028,8 +8077,8 @@ SWIGINTERN PyObject *_wrap_rderivs_d2phs_get(PyObject *SWIGUNUSEDPARM(self), PyO SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_d2phs_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->d2phs); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->d2phs); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -8039,10 +8088,10 @@ SWIGINTERN PyObject *_wrap_rderivs_d2phs_get(PyObject *SWIGUNUSEDPARM(self), PyO SWIGINTERN PyObject *_wrap_rderivs_locpow_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; struct RDERIVS *arg1 = (struct RDERIVS *) 0 ; - float arg2 ; + double arg2 ; void *argp1 = 0 ; int res1 = 0 ; - float val2 ; + double val2 ; int ecode2 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; @@ -8053,11 +8102,11 @@ SWIGINTERN PyObject *_wrap_rderivs_locpow_set(PyObject *SWIGUNUSEDPARM(self), Py SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_locpow_set" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - ecode2 = SWIG_AsVal_float(obj1, &val2); + ecode2 = SWIG_AsVal_double(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_locpow_set" "', argument " "2"" of type '" "float""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "rderivs_locpow_set" "', argument " "2"" of type '" "double""'"); } - arg2 = (float)(val2); + arg2 = (double)(val2); if (arg1) (arg1)->locpow = arg2; resultobj = SWIG_Py_Void(); return resultobj; @@ -8072,7 +8121,7 @@ SWIGINTERN PyObject *_wrap_rderivs_locpow_get(PyObject *SWIGUNUSEDPARM(self), Py void *argp1 = 0 ; int res1 = 0 ; PyObject * obj0 = 0 ; - float result; + double result; if (!PyArg_ParseTuple(args,(char *)"O:rderivs_locpow_get",&obj0)) SWIG_fail; res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_RDERIVS, 0 | 0 ); @@ -8080,8 +8129,8 @@ SWIGINTERN PyObject *_wrap_rderivs_locpow_get(PyObject *SWIGUNUSEDPARM(self), Py SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "rderivs_locpow_get" "', argument " "1"" of type '" "struct RDERIVS *""'"); } arg1 = (struct RDERIVS *)(argp1); - result = (float) ((arg1)->locpow); - resultobj = SWIG_From_float((float)(result)); + result = (double) ((arg1)->locpow); + resultobj = SWIG_From_double((double)(result)); return resultobj; fail: return NULL; @@ -11140,31 +11189,30 @@ SWIGINTERN PyObject *_wrap_gen_bin_response(PyObject *SWIGUNUSEDPARM(self), PyOb SWIGINTERN PyObject *_wrap_get_localpower(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; fcomplex *arg1 = (fcomplex *) 0 ; - int arg2 ; + long arg2 ; double arg3 ; - void *argp1 = 0 ; - int res1 = 0 ; - int val2 ; - int ecode2 = 0 ; + PyArrayObject *array1 = NULL ; + int is_new_object1 = 0 ; double val3 ; int ecode3 = 0 ; PyObject * obj0 = 0 ; PyObject * obj1 = 0 ; - PyObject * obj2 = 0 ; float result; - if (!PyArg_ParseTuple(args,(char *)"OOO:get_localpower",&obj0,&obj1,&obj2)) SWIG_fail; - res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_FCOMPLEX, 0 | 0 ); - if (!SWIG_IsOK(res1)) { - SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "get_localpower" "', argument " "1"" of type '" "fcomplex *""'"); + if (!PyArg_ParseTuple(args,(char *)"OO:get_localpower",&obj0,&obj1)) SWIG_fail; + { + npy_intp size[1] = { + -1 + }; + array1 = obj_to_array_contiguous_allow_conversion(obj0, + NPY_CFLOAT, + &is_new_object1); + if (!array1 || !require_dimensions(array1, 1) || + !require_size(array1, size, 1)) SWIG_fail; + arg1 = (fcomplex*) array_data(array1); + arg2 = (long) array_size(array1,0); } - arg1 = (fcomplex *)(argp1); - ecode2 = SWIG_AsVal_int(obj1, &val2); - if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "get_localpower" "', argument " "2"" of type '" "int""'"); - } - arg2 = (int)(val2); - ecode3 = SWIG_AsVal_double(obj2, &val3); + ecode3 = SWIG_AsVal_double(obj1, &val3); if (!SWIG_IsOK(ecode3)) { SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "get_localpower" "', argument " "3"" of type '" "double""'"); } @@ -11187,8 +11235,20 @@ SWIGINTERN PyObject *_wrap_get_localpower(PyObject *SWIGUNUSEDPARM(self), PyObje } } resultobj = SWIG_From_float((float)(result)); + { + if (is_new_object1 && array1) + { + Py_DECREF(array1); + } + } return resultobj; fail: + { + if (is_new_object1 && array1) + { + Py_DECREF(array1); + } + } return NULL; } @@ -11196,14 +11256,12 @@ SWIGINTERN PyObject *_wrap_get_localpower(PyObject *SWIGUNUSEDPARM(self), PyObje SWIGINTERN PyObject *_wrap_get_localpower3d(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; fcomplex *arg1 = (fcomplex *) 0 ; - int arg2 ; + long arg2 ; double arg3 ; double arg4 ; double arg5 ; - void *argp1 = 0 ; - int res1 = 0 ; - int val2 ; - int ecode2 = 0 ; + PyArrayObject *array1 = NULL ; + int is_new_object1 = 0 ; double val3 ; int ecode3 = 0 ; double val4 ; @@ -11214,31 +11272,32 @@ SWIGINTERN PyObject *_wrap_get_localpower3d(PyObject *SWIGUNUSEDPARM(self), PyOb PyObject * obj1 = 0 ; PyObject * obj2 = 0 ; PyObject * obj3 = 0 ; - PyObject * obj4 = 0 ; float result; - if (!PyArg_ParseTuple(args,(char *)"OOOOO:get_localpower3d",&obj0,&obj1,&obj2,&obj3,&obj4)) SWIG_fail; - res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_FCOMPLEX, 0 | 0 ); - if (!SWIG_IsOK(res1)) { - SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "get_localpower3d" "', argument " "1"" of type '" "fcomplex *""'"); + if (!PyArg_ParseTuple(args,(char *)"OOOO:get_localpower3d",&obj0,&obj1,&obj2,&obj3)) SWIG_fail; + { + npy_intp size[1] = { + -1 + }; + array1 = obj_to_array_contiguous_allow_conversion(obj0, + NPY_CFLOAT, + &is_new_object1); + if (!array1 || !require_dimensions(array1, 1) || + !require_size(array1, size, 1)) SWIG_fail; + arg1 = (fcomplex*) array_data(array1); + arg2 = (long) array_size(array1,0); } - arg1 = (fcomplex *)(argp1); - ecode2 = SWIG_AsVal_int(obj1, &val2); - if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "get_localpower3d" "', argument " "2"" of type '" "int""'"); - } - arg2 = (int)(val2); - ecode3 = SWIG_AsVal_double(obj2, &val3); + ecode3 = SWIG_AsVal_double(obj1, &val3); if (!SWIG_IsOK(ecode3)) { SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "get_localpower3d" "', argument " "3"" of type '" "double""'"); } arg3 = (double)(val3); - ecode4 = SWIG_AsVal_double(obj3, &val4); + ecode4 = SWIG_AsVal_double(obj2, &val4); if (!SWIG_IsOK(ecode4)) { SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "get_localpower3d" "', argument " "4"" of type '" "double""'"); } arg4 = (double)(val4); - ecode5 = SWIG_AsVal_double(obj4, &val5); + ecode5 = SWIG_AsVal_double(obj3, &val5); if (!SWIG_IsOK(ecode5)) { SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "get_localpower3d" "', argument " "5"" of type '" "double""'"); } @@ -11261,8 +11320,20 @@ SWIGINTERN PyObject *_wrap_get_localpower3d(PyObject *SWIGUNUSEDPARM(self), PyOb } } resultobj = SWIG_From_float((float)(result)); + { + if (is_new_object1 && array1) + { + Py_DECREF(array1); + } + } return resultobj; fail: + { + if (is_new_object1 && array1) + { + Py_DECREF(array1); + } + } return NULL; } @@ -11270,7 +11341,7 @@ SWIGINTERN PyObject *_wrap_get_localpower3d(PyObject *SWIGUNUSEDPARM(self), PyOb SWIGINTERN PyObject *_wrap_get_derivs3d(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; fcomplex *arg1 = (fcomplex *) 0 ; - int arg2 ; + long arg2 ; double arg3 ; double arg4 ; double arg5 ; @@ -11278,7 +11349,7 @@ SWIGINTERN PyObject *_wrap_get_derivs3d(PyObject *SWIGUNUSEDPARM(self), PyObject rderivs *arg7 = (rderivs *) 0 ; void *argp1 = 0 ; int res1 = 0 ; - int val2 ; + long val2 ; int ecode2 = 0 ; double val3 ; int ecode3 = 0 ; @@ -11304,11 +11375,11 @@ SWIGINTERN PyObject *_wrap_get_derivs3d(PyObject *SWIGUNUSEDPARM(self), PyObject SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "get_derivs3d" "', argument " "1"" of type '" "fcomplex *""'"); } arg1 = (fcomplex *)(argp1); - ecode2 = SWIG_AsVal_int(obj1, &val2); + ecode2 = SWIG_AsVal_long(obj1, &val2); if (!SWIG_IsOK(ecode2)) { - SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "get_derivs3d" "', argument " "2"" of type '" "int""'"); + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "get_derivs3d" "', argument " "2"" of type '" "long""'"); } - arg2 = (int)(val2); + arg2 = (long)(val2); ecode3 = SWIG_AsVal_double(obj2, &val3); if (!SWIG_IsOK(ecode3)) { SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "get_derivs3d" "', argument " "3"" of type '" "double""'"); @@ -13216,6 +13287,209 @@ SWIGINTERN PyObject *_wrap_max_rz_arr(PyObject *SWIGUNUSEDPARM(self), PyObject * } +SWIGINTERN PyObject *_wrap_max_rz_arr_harmonics(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { + PyObject *resultobj = 0; + fcomplex *arg1 = (fcomplex *) 0 ; + long arg2 ; + double arg3 ; + double arg4 ; + double *arg5 = (double *) 0 ; + int arg6 ; + double *arg7 = (double *) 0 ; + double *arg8 = (double *) 0 ; + PyArrayObject *array1 = NULL ; + int i1 = 1 ; + double val3 ; + int ecode3 = 0 ; + double val4 ; + int ecode4 = 0 ; + PyArrayObject *array5 = NULL ; + int i5 = 1 ; + double temp7 ; + int res7 = SWIG_TMPOBJ ; + double temp8 ; + int res8 = SWIG_TMPOBJ ; + PyObject * obj0 = 0 ; + PyObject * obj1 = 0 ; + PyObject * obj2 = 0 ; + PyObject * obj3 = 0 ; + + arg7 = &temp7; + arg8 = &temp8; + if (!PyArg_ParseTuple(args,(char *)"OOOO:max_rz_arr_harmonics",&obj0,&obj1,&obj2,&obj3)) SWIG_fail; + { + array1 = obj_to_array_no_conversion(obj0, NPY_CFLOAT); + if (!array1 || !require_dimensions(array1,1) || !require_contiguous(array1) + || !require_native(array1)) SWIG_fail; + arg1 = (fcomplex*) array_data(array1); + arg2 = 1; + for (i1=0; i1 < array_numdims(array1); ++i1) arg2 *= array_size(array1,i1); + } + ecode3 = SWIG_AsVal_double(obj1, &val3); + if (!SWIG_IsOK(ecode3)) { + SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "max_rz_arr_harmonics" "', argument " "3"" of type '" "double""'"); + } + arg3 = (double)(val3); + ecode4 = SWIG_AsVal_double(obj2, &val4); + if (!SWIG_IsOK(ecode4)) { + SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "max_rz_arr_harmonics" "', argument " "4"" of type '" "double""'"); + } + arg4 = (double)(val4); + { + array5 = obj_to_array_no_conversion(obj3, NPY_DOUBLE); + if (!array5 || !require_dimensions(array5,1) || !require_contiguous(array5) + || !require_native(array5)) SWIG_fail; + arg5 = (double*) array_data(array5); + arg6 = 1; + for (i5=0; i5 < array_numdims(array5); ++i5) arg6 *= array_size(array5,i5); + } + { + errno = 0; + wrap_max_rz_arr_harmonics(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8); + + if (errno != 0) + { + switch(errno) + { + case ENOMEM: + PyErr_Format(PyExc_MemoryError, "Failed malloc()"); + break; + default: + PyErr_Format(PyExc_Exception, "Unknown exception"); + } + SWIG_fail; + } + } + resultobj = SWIG_Py_Void(); + if (SWIG_IsTmpObj(res7)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg7))); + } else { + int new_flags = SWIG_IsNewObj(res7) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg7), SWIGTYPE_p_double, new_flags)); + } + if (SWIG_IsTmpObj(res8)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg8))); + } else { + int new_flags = SWIG_IsNewObj(res8) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg8), SWIGTYPE_p_double, new_flags)); + } + return resultobj; +fail: + return NULL; +} + + +SWIGINTERN PyObject *_wrap_max_rzw_arr_harmonics(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { + PyObject *resultobj = 0; + fcomplex *arg1 = (fcomplex *) 0 ; + long arg2 ; + double arg3 ; + double arg4 ; + double arg5 ; + double *arg6 = (double *) 0 ; + int arg7 ; + double *arg8 = (double *) 0 ; + double *arg9 = (double *) 0 ; + double *arg10 = (double *) 0 ; + PyArrayObject *array1 = NULL ; + int i1 = 1 ; + double val3 ; + int ecode3 = 0 ; + double val4 ; + int ecode4 = 0 ; + double val5 ; + int ecode5 = 0 ; + PyArrayObject *array6 = NULL ; + int i6 = 1 ; + double temp8 ; + int res8 = SWIG_TMPOBJ ; + double temp9 ; + int res9 = SWIG_TMPOBJ ; + double temp10 ; + int res10 = SWIG_TMPOBJ ; + PyObject * obj0 = 0 ; + PyObject * obj1 = 0 ; + PyObject * obj2 = 0 ; + PyObject * obj3 = 0 ; + PyObject * obj4 = 0 ; + + arg8 = &temp8; + arg9 = &temp9; + arg10 = &temp10; + if (!PyArg_ParseTuple(args,(char *)"OOOOO:max_rzw_arr_harmonics",&obj0,&obj1,&obj2,&obj3,&obj4)) SWIG_fail; + { + array1 = obj_to_array_no_conversion(obj0, NPY_CFLOAT); + if (!array1 || !require_dimensions(array1,1) || !require_contiguous(array1) + || !require_native(array1)) SWIG_fail; + arg1 = (fcomplex*) array_data(array1); + arg2 = 1; + for (i1=0; i1 < array_numdims(array1); ++i1) arg2 *= array_size(array1,i1); + } + ecode3 = SWIG_AsVal_double(obj1, &val3); + if (!SWIG_IsOK(ecode3)) { + SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "max_rzw_arr_harmonics" "', argument " "3"" of type '" "double""'"); + } + arg3 = (double)(val3); + ecode4 = SWIG_AsVal_double(obj2, &val4); + if (!SWIG_IsOK(ecode4)) { + SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "max_rzw_arr_harmonics" "', argument " "4"" of type '" "double""'"); + } + arg4 = (double)(val4); + ecode5 = SWIG_AsVal_double(obj3, &val5); + if (!SWIG_IsOK(ecode5)) { + SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "max_rzw_arr_harmonics" "', argument " "5"" of type '" "double""'"); + } + arg5 = (double)(val5); + { + array6 = obj_to_array_no_conversion(obj4, NPY_DOUBLE); + if (!array6 || !require_dimensions(array6,1) || !require_contiguous(array6) + || !require_native(array6)) SWIG_fail; + arg6 = (double*) array_data(array6); + arg7 = 1; + for (i6=0; i6 < array_numdims(array6); ++i6) arg7 *= array_size(array6,i6); + } + { + errno = 0; + wrap_max_rzw_arr_harmonics(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9,arg10); + + if (errno != 0) + { + switch(errno) + { + case ENOMEM: + PyErr_Format(PyExc_MemoryError, "Failed malloc()"); + break; + default: + PyErr_Format(PyExc_Exception, "Unknown exception"); + } + SWIG_fail; + } + } + resultobj = SWIG_Py_Void(); + if (SWIG_IsTmpObj(res8)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg8))); + } else { + int new_flags = SWIG_IsNewObj(res8) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg8), SWIGTYPE_p_double, new_flags)); + } + if (SWIG_IsTmpObj(res9)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg9))); + } else { + int new_flags = SWIG_IsNewObj(res9) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg9), SWIGTYPE_p_double, new_flags)); + } + if (SWIG_IsTmpObj(res10)) { + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_From_double((*arg10))); + } else { + int new_flags = SWIG_IsNewObj(res10) ? (SWIG_POINTER_OWN | 0 ) : 0 ; + resultobj = SWIG_Python_AppendOutput(resultobj, SWIG_NewPointerObj((void*)(arg10), SWIGTYPE_p_double, new_flags)); + } + return resultobj; +fail: + return NULL; +} + + SWIGINTERN PyObject *_wrap_max_rzw_arr(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; fcomplex *arg1 = (fcomplex *) 0 ; @@ -13706,6 +13980,8 @@ static PyMethodDef SwigMethods[] = { { (char *)"corr_rzw_vol", _wrap_corr_rzw_vol, METH_VARARGS, NULL}, { (char *)"max_r_arr", _wrap_max_r_arr, METH_VARARGS, NULL}, { (char *)"max_rz_arr", _wrap_max_rz_arr, METH_VARARGS, NULL}, + { (char *)"max_rz_arr_harmonics", _wrap_max_rz_arr_harmonics, METH_VARARGS, NULL}, + { (char *)"max_rzw_arr_harmonics", _wrap_max_rzw_arr_harmonics, METH_VARARGS, NULL}, { (char *)"max_rzw_arr", _wrap_max_rzw_arr, METH_VARARGS, NULL}, { (char *)"barycenter", _wrap_barycenter, METH_VARARGS, NULL}, { NULL, NULL, 0, NULL } diff --git a/src/accel_utils.c b/src/accel_utils.c index 87159b678..300cb1712 100644 --- a/src/accel_utils.c +++ b/src/accel_utils.c @@ -50,6 +50,16 @@ static inline int twon_to_index(int n) return x; } + +static inline double calc_required_r(double harm_fract, double rfull) +/* Calculate the 'r' you need for subharmonic */ +/* harm_fract = harmnum / numharm if the */ +/* 'r' at the fundamental harmonic is 'rfull'. */ +{ + return rint(ACCEL_RDR * rfull * harm_fract) * ACCEL_DR; +} + + static inline int calc_required_z(double harm_fract, double zfull) /* Calculate the 'z' you need for subharmonic */ /* harm_fract = harmnum / numharm if the */ @@ -59,12 +69,12 @@ static inline int calc_required_z(double harm_fract, double zfull) } -static inline double calc_required_r(double harm_fract, double rfull) -/* Calculate the 'r' you need for subharmonic */ +static inline int calc_required_w(double harm_fract, double wfull) +/* Calculate the maximum 'w' needed for the given subharmonic */ /* harm_fract = harmnum / numharm if the */ -/* 'r' at the fundamental harmonic is 'rfull'. */ +/* 'w' at the fundamental harmonic is 'wfull'. */ { - return (int) (ACCEL_RDR * rfull * harm_fract + 0.5) * ACCEL_DR; + return NEAREST_INT(ACCEL_RDW * wfull * harm_fract) * ACCEL_DW; } @@ -78,12 +88,20 @@ static inline int index_from_r(double r, double lor) static inline int index_from_z(double z, double loz) /* Return an index for a Fourier Fdot given an array that */ -/* has stepsize ACCEL_DZ and low freq 'lor'. */ +/* has stepsize ACCEL_DZ and low freq dot 'loz'. */ { return (int) ((z - loz) * ACCEL_RDZ + DBLCORRECT); } +static inline int index_from_w(double w, double low) +/* Return an index for a Fourier Fdotdot given an array that */ +/* has stepsize ACCEL_DW and low freq dotdot 'low'. */ +{ + return (int) ((w - low) * ACCEL_RDW + DBLCORRECT); +} + + static void compare_rzw_cands(fourierprops * list, int nlist, char *notes) { int ii, jj, kk; @@ -116,35 +134,35 @@ static void compare_rzw_cands(fourierprops * list, int nlist, char *notes) } -static int calc_fftlen(int numharm, int harmnum, int max_zfull) +static int calc_fftlen(int numharm, int harmnum, int max_zfull, int max_wfull, accelobs * obs) /* The fft length needed to properly process a subharmonic */ { int bins_needed, end_effects; double harm_fract; harm_fract = (double) harmnum / (double) numharm; - bins_needed = (ACCEL_USELEN * harmnum) / numharm + 2; + bins_needed = (int) ceil(obs->corr_uselen * harm_fract) + 2; end_effects = 2 * ACCEL_NUMBETWEEN * - z_resp_halfwidth(calc_required_z(harm_fract, max_zfull), LOWACC); - //printf("bins_needed = %d end_effects = %d FFTlen = %lld\n", - // bins_needed, end_effects, next2_to_n(bins_needed + end_effects)); - return next2_to_n(bins_needed + end_effects); + w_resp_halfwidth(calc_required_z(harm_fract, max_zfull), + calc_required_w(harm_fract, max_wfull), LOWACC); + return next_good_fftlen(bins_needed + end_effects); } -static void init_kernel(int z, int fftlen, kernel * kern) +static void init_kernel(int z, int w, int fftlen, kernel * kern) { int numkern; fcomplex *tempkern; kern->z = z; + kern->w = w; kern->fftlen = fftlen; kern->numbetween = ACCEL_NUMBETWEEN; - kern->kern_half_width = z_resp_halfwidth((double) z, LOWACC); + kern->kern_half_width = w_resp_halfwidth((double) z, (double) w, LOWACC); numkern = 2 * kern->numbetween * kern->kern_half_width; kern->numgoodbins = kern->fftlen - numkern; kern->data = gen_cvect(kern->fftlen); - tempkern = gen_z_response(0.0, kern->numbetween, kern->z, numkern); + tempkern = gen_w_response(0.0, kern->numbetween, kern->z, kern->w, numkern); place_complex_kernel(tempkern, numkern, kern->data, kern->fftlen); vect_free(tempkern); COMPLEXFFT(kern->data, kern->fftlen, -1); @@ -157,67 +175,122 @@ static void free_kernel(kernel * kern) } -static void init_subharminfo(int numharm, int harmnum, int zmax, subharminfo * shi) -/* Note: 'zmax' is the overall maximum 'z' in the search */ +kernel **gen_kernmatrix(int numz, int numw) { + int ii; + kernel **kerns; + + kerns = (kernel **) malloc((size_t) numw * sizeof(kernel *)); + if (!kerns) { + perror("\nError in 1st malloc() in gen_kernmatrix()"); + printf("\n"); + exit(-1); + } + kerns[0] = (kernel *) malloc((size_t) ((numz * numw) * sizeof(kernel))); + if (!kerns[0]) { + perror("\nError in 2nd malloc() in gen_kernmatrix()"); + printf("\n"); + exit(-1); + } + for (ii = 1; ii < numw; ii++) + kerns[ii] = kerns[ii - 1] + numz; + return kerns; +} + + +static void init_subharminfo(int numharm, int harmnum, int zmax, int wmax, subharminfo * shi, accelobs * obs) +/* Note: 'zmax' is the overall maximum 'z' in the search while + 'wmax' is the overall maximum 'w' in the search */ { - int ii, fftlen; + int ii, jj, fftlen; double harm_fract; harm_fract = (double) harmnum / (double) numharm; shi->numharm = numharm; shi->harmnum = harmnum; shi->zmax = calc_required_z(harm_fract, zmax); - if (numharm > 1) - shi->rinds = - (unsigned short *) malloc(ACCEL_USELEN * sizeof(unsigned short)); - fftlen = calc_fftlen(numharm, harmnum, zmax); - shi->numkern = (shi->zmax / ACCEL_DZ) * 2 + 1; - shi->kern = (kernel *) malloc(shi->numkern * sizeof(kernel)); - for (ii = 0; ii < shi->numkern; ii++) - init_kernel(-shi->zmax + ii * ACCEL_DZ, fftlen, &shi->kern[ii]); + shi->wmax = calc_required_w(harm_fract, wmax); + if (numharm > 1) { + shi->rinds = (unsigned short *) malloc(obs->corr_uselen * sizeof(unsigned short)); + shi->zinds = (unsigned short *) malloc(obs->corr_uselen * sizeof(unsigned short)); + } + if (numharm==1 && harmnum==1) + fftlen = obs->fftlen; + else + fftlen = calc_fftlen(numharm, harmnum, zmax, wmax, obs); + shi->numkern_zdim = (shi->zmax / ACCEL_DZ) * 2 + 1; + shi->numkern_wdim = (shi->wmax / ACCEL_DW) * 2 + 1; + shi->numkern = shi->numkern_zdim * shi->numkern_wdim; + /* Allocate 2D array of kernels, with dimensions being z and w */ + shi->kern = gen_kernmatrix(shi->numkern_zdim, shi->numkern_wdim); + /* Actually append kernels to each array element */ + for (ii = 0; ii < shi->numkern_wdim; ii++) { + for (jj = 0; jj < shi->numkern_zdim; jj++) { + init_kernel(-shi->zmax + jj * ACCEL_DZ, + -shi->wmax + ii * ACCEL_DW, fftlen, &shi->kern[ii][jj]); + } + } } subharminfo **create_subharminfos(accelobs * obs) { - int ii, jj, harmtosum; + double kern_ram_use=0; + int ii, jj, harmtosum, fftlen; subharminfo **shis; - + shis = (subharminfo **) malloc(obs->numharmstages * sizeof(subharminfo *)); /* Prep the fundamental (actually, the highest harmonic) */ shis[0] = (subharminfo *) malloc(2 * sizeof(subharminfo)); - init_subharminfo(1, 1, (int) obs->zhi, &shis[0][0]); - printf - (" Harmonic 1/1 has %3d kernel(s) from z = %4d to %4d, FFT length = %d\n", - shis[0][0].numkern, -shis[0][0].zmax, shis[0][0].zmax, - calc_fftlen(1, 1, (int) obs->zhi)); + init_subharminfo(1, 1, (int) obs->zhi, (int) obs->whi, &shis[0][0], obs); + fftlen = obs->fftlen; + kern_ram_use += shis[0][0].numkern * fftlen * sizeof(fcomplex); // in Bytes + if (obs->numw) + printf(" Harm 1/1 : %5d kernels, %4d < z < %-4d and %5d < w < %-5d (%5d pt FFTs)\n", + shis[0][0].numkern, -shis[0][0].zmax, shis[0][0].zmax, + -shis[0][0].wmax, shis[0][0].wmax, fftlen); + else + printf(" Harm 1/1 : %5d kernels, %4d < z < %-4d (%d pt FFTs)\n", + shis[0][0].numkern, -shis[0][0].zmax, shis[0][0].zmax, fftlen); /* Prep the sub-harmonics if needed */ if (!obs->inmem) { for (ii = 1; ii < obs->numharmstages; ii++) { harmtosum = index_to_twon(ii); shis[ii] = (subharminfo *) malloc(harmtosum * sizeof(subharminfo)); for (jj = 1; jj < harmtosum; jj += 2) { - init_subharminfo(harmtosum, jj, (int) obs->zhi, &shis[ii][jj - 1]); - printf - (" Harmonic %2d/%-2d has %3d kernel(s) from z = %4d to %4d, FFT length = %d\n", - jj, harmtosum, shis[ii][jj - 1].numkern, -shis[ii][jj - 1].zmax, - shis[ii][jj - 1].zmax, - calc_fftlen(harmtosum, jj, (int) obs->zhi)); + init_subharminfo(harmtosum, jj, (int) obs->zhi, + (int) obs->whi, &shis[ii][jj - 1], obs); + fftlen = calc_fftlen(harmtosum, jj, (int) obs->zhi, (int) obs->whi, obs); + kern_ram_use += shis[ii][jj - 1].numkern * fftlen * sizeof(fcomplex); // in Bytes + if (obs->numw) + printf(" Harm %2d/%-2d: %5d kernels, %4d < z < %-4d and %5d < w < %-5d (%5d pt FFTs)\n", + jj, harmtosum, shis[ii][jj - 1].numkern, + -shis[ii][jj - 1].zmax, shis[ii][jj - 1].zmax, + -shis[ii][jj - 1].wmax, shis[ii][jj - 1].wmax, fftlen); + else + printf(" Harm %2d/%-2d: %5d kernels, %4d < z < %-4d (%d pt FFTs)\n", + jj, harmtosum, shis[ii][jj - 1].numkern, + -shis[ii][jj - 1].zmax, shis[ii][jj - 1].zmax, fftlen); } } } + printf("Total RAM used by correlation kernels: %.3f GB\n", kern_ram_use / (1 << 30)); return shis; } static void free_subharminfo(subharminfo * shi) { - int ii; - - for (ii = 0; ii < shi->numkern; ii++) - free_kernel(&shi->kern[ii]); - if (shi->numharm > 1) + int ii, jj; + + for (ii = 0; ii < shi->numkern_wdim; ii++) { + for (jj = 0; jj < shi->numkern_zdim; jj++) { + free_kernel(&shi->kern[ii][jj]); + } + } + if (shi->numharm > 1) { free(shi->rinds); + free(shi->zinds); + } free(shi->kern); } @@ -225,7 +298,7 @@ static void free_subharminfo(subharminfo * shi) void free_subharminfos(accelobs * obs, subharminfo ** shis) { int ii, jj, harmtosum; - + /* Free the sub-harmonics */ if (!obs->inmem) { for (ii = 1; ii < obs->numharmstages; ii++) { @@ -245,19 +318,21 @@ void free_subharminfos(accelobs * obs, subharminfo ** shis) static accelcand *create_accelcand(float power, float sigma, - int numharm, double r, double z) + int numharm, double r, double z, double w) { accelcand *obj; - + obj = (accelcand *) malloc(sizeof(accelcand)); obj->power = power; obj->sigma = sigma; obj->numharm = numharm; obj->r = r; obj->z = z; + obj->w = w; obj->pows = NULL; obj->hirs = NULL; obj->hizs = NULL; + obj->hiws = NULL; obj->derivs = NULL; return obj; } @@ -269,6 +344,7 @@ void free_accelcand(gpointer data, gpointer user_data) vect_free(((accelcand *) data)->pows); vect_free(((accelcand *) data)->hirs); vect_free(((accelcand *) data)->hizs); + vect_free(((accelcand *) data)->hiws); free(((accelcand *) data)->derivs); } free((accelcand *) data); @@ -299,7 +375,7 @@ GSList *sort_accelcands(GSList * list) static GSList *insert_new_accelcand(GSList * list, float power, float sigma, - int numharm, double rr, double zz, int *added) + int numharm, double rr, double zz, double ww, int *added) /* Checks the current list to see if there is already */ /* a candidate within ACCEL_CLOSEST_R bins. If not, */ /* it adds it to the list in increasing freq order. */ @@ -311,7 +387,7 @@ static GSList *insert_new_accelcand(GSList * list, float power, float sigma, if (!list) { new_list = g_slist_alloc(); new_list->data = - (gpointer *) create_accelcand(power, sigma, numharm, rr, zz); + (gpointer *) create_accelcand(power, sigma, numharm, rr, zz, ww); *added = 1; return new_list; } @@ -333,7 +409,7 @@ static GSList *insert_new_accelcand(GSList * list, float power, float sigma, if (((accelcand *) (prev_list->data))->sigma < sigma) { free_accelcand(prev_list->data, NULL); prev_list->data = (gpointer *) create_accelcand(power, sigma, - numharm, rr, zz); + numharm, rr, zz, ww); *added = 1; } if (next_diff_r < ACCEL_CLOSEST_R) { @@ -346,7 +422,7 @@ static GSList *insert_new_accelcand(GSList * list, float power, float sigma, } else { /* Overwrite the next cand */ tmp_list->data = (gpointer *) create_accelcand(power, sigma, - numharm, rr, zz); + numharm, rr, zz, ww); *added = 1; } } @@ -356,13 +432,13 @@ static GSList *insert_new_accelcand(GSList * list, float power, float sigma, if (((accelcand *) (tmp_list->data))->sigma < sigma) { free_accelcand(tmp_list->data, NULL); tmp_list->data = (gpointer *) create_accelcand(power, sigma, - numharm, rr, zz); + numharm, rr, zz, ww); *added = 1; } } else { /* This is a new candidate */ new_list = g_slist_alloc(); new_list->data = - (gpointer *) create_accelcand(power, sigma, numharm, rr, zz); + (gpointer *) create_accelcand(power, sigma, numharm, rr, zz, ww); *added = 1; if (!tmp_list->next && (((accelcand *) (tmp_list->data))->r < (rr - ACCEL_CLOSEST_R))) { @@ -461,65 +537,78 @@ GSList *eliminate_harmonics(GSList * cands, int *numcands) } -// FIXME: this shouldn't be a #define, or it shouldn't be here void optimize_accelcand(accelcand * cand, accelobs * obs) { int ii; - int *r_offset; - fcomplex **data; - double r, z; + double r, z, w; cand->pows = gen_dvect(cand->numharm); cand->hirs = gen_dvect(cand->numharm); cand->hizs = gen_dvect(cand->numharm); - r_offset = (int *) malloc(sizeof(int) * cand->numharm); - data = (fcomplex **) malloc(sizeof(fcomplex *) * cand->numharm); + cand->hiws = gen_dvect(cand->numharm); cand->derivs = (rderivs *) malloc(sizeof(rderivs) * cand->numharm); - - if (obs->use_harmonic_polishing) { - if (obs->mmap_file || obs->dat_input) { - for (ii = 0; ii < cand->numharm; ii++) { - r_offset[ii] = obs->lobin; - data[ii] = obs->fft; - } - max_rz_arr_harmonics(data, - cand->numharm, - r_offset, - obs->numbins, - cand->r - obs->lobin, - cand->z, &r, &z, cand->derivs, cand->pows); - } else { - max_rz_file_harmonics(obs->fftfile, + + if (obs->use_harmonic_polishing && + (obs->mmap_file || obs->dat_input)) { + if (obs->numw) { + max_rzw_arr_harmonics(obs->fft, obs->numbins, cand->numharm, - obs->lobin, cand->r - obs->lobin, - cand->z, &r, &z, cand->derivs, cand->pows); + cand->z, cand->w, &r, &z, &w, + cand->derivs, cand->pows); + + } else { + max_rz_arr_harmonics(obs->fft, obs->numbins, + cand->numharm, + cand->r - obs->lobin, + cand->z, &r, &z, + cand->derivs, cand->pows); } for (ii = 0; ii < cand->numharm; ii++) { cand->hirs[ii] = (r + obs->lobin) * (ii + 1); cand->hizs[ii] = z * (ii + 1); + cand->hiws[ii] = obs->numw ? w * (ii + 1) : 0.0; } } else { for (ii = 0; ii < cand->numharm; ii++) { - if (obs->mmap_file || obs->dat_input) - cand->pows[ii] = max_rz_arr(obs->fft, - obs->numbins, - cand->r * (ii + 1) - obs->lobin, - cand->z * (ii + 1), - &(cand->hirs[ii]), - &(cand->hizs[ii]), &(cand->derivs[ii])); - else - cand->pows[ii] = max_rz_file(obs->fftfile, - cand->r * (ii + 1) - obs->lobin, - cand->z * (ii + 1), - &(cand->hirs[ii]), - &(cand->hizs[ii]), &(cand->derivs[ii])); + if (obs->mmap_file || obs->dat_input) { + if (obs->numw) + cand->pows[ii] = max_rzw_arr(obs->fft, + obs->numbins, + cand->r * (ii + 1) - obs->lobin, + cand->z * (ii + 1), + cand->w * (ii + 1), + &(cand->hirs[ii]), + &(cand->hizs[ii]), + &(cand->hiws[ii]), + &(cand->derivs[ii])); + else + cand->pows[ii] = max_rz_arr(obs->fft, + obs->numbins, + cand->r * (ii + 1) - obs->lobin, + cand->z * (ii + 1), + &(cand->hirs[ii]), + &(cand->hizs[ii]), &(cand->derivs[ii])); + } else { + if (obs->numw) + cand->pows[ii] = max_rzw_file(obs->fftfile, + cand->r * (ii + 1) - obs->lobin, + cand->z * (ii + 1), + cand->w * (ii + 1), + &(cand->hirs[ii]), + &(cand->hizs[ii]), + &(cand->hiws[ii]), + &(cand->derivs[ii])); + else + cand->pows[ii] = max_rz_file(obs->fftfile, + cand->r * (ii + 1) - obs->lobin, + cand->z * (ii + 1), + &(cand->hirs[ii]), + &(cand->hizs[ii]), &(cand->derivs[ii])); + } cand->hirs[ii] += obs->lobin; } } - free(r_offset); - free(data); - cand->sigma = candidate_sigma(cand->power, cand->numharm, obs->numindep[twon_to_index(cand->numharm)]); } @@ -566,20 +655,20 @@ void output_fundamentals(fourierprops * props, GSList * list, accelobs * obs, infodata * idata) { double accel = 0.0, accelerr = 0.0, coherent_pow; - int ii, jj, numcols = 12, numcands, *width, *error; - int widths[12] = { 4, 5, 6, 8, 4, 16, 15, 15, 15, 11, 15, 20 }; - int errors[12] = { 0, 0, 0, 0, 0, 1, 1, 2, 1, 2, 2, 0 }; + int ii, jj, numcols = 13, numcands, *width, *error; + int widths[13] = { 4, 5, 6, 8, 4, 16, 15, 15, 15, 11, 11, 15, 20 }; + int errors[13] = { 0, 0, 0, 0, 0, 1, 1, 2, 1, 2, 2, 2, 0 }; char tmpstr[30], ctrstr[30], *notes; accelcand *cand; GSList *listptr; rzwerrs errs; static char **title; static char *titles1[] = { "", "", "Summed", "Coherent", "Num", "Period", - "Frequency", "FFT 'r'", "Freq Deriv", "FFT 'z'", + "Frequency", "FFT 'r'", "Freq Deriv", "FFT 'z'", "FFT 'w'", "Accel", "" }; static char *titles2[] = { "Cand", "Sigma", "Power", "Power", "Harm", "(ms)", - "(Hz)", "(bin)", "(Hz/s)", "(bins)", + "(Hz)", "(bin)", "(Hz/s)", "(bins)", "(bins)", "(m/s^2)", "Notes" }; @@ -617,8 +706,14 @@ void output_fundamentals(fourierprops * props, GSList * list, width = widths; title = titles1; for (ii = 0; ii < numcols - 1; ii++) { - center_string(ctrstr, *title++, *width++); - fprintf(obs->workfile, "%s ", ctrstr); + if (obs->numw==0 && ii==10) { // Skip jerk parts + title++; + width++; + continue; + } else { + center_string(ctrstr, *title++, *width++); + fprintf(obs->workfile, "%s ", ctrstr); + } } center_string(ctrstr, *title++, *width++); fprintf(obs->workfile, "%s\n", ctrstr); @@ -626,17 +721,28 @@ void output_fundamentals(fourierprops * props, GSList * list, width = widths; title = titles2; for (ii = 0; ii < numcols - 1; ii++) { - center_string(ctrstr, *title++, *width++); - fprintf(obs->workfile, "%s ", ctrstr); + if (obs->numw==0 && ii==10) { // Skip jerk parts + title++; + width++; + continue; + } else { + center_string(ctrstr, *title++, *width++); + fprintf(obs->workfile, "%s ", ctrstr); + } } center_string(ctrstr, *title++, *width++); fprintf(obs->workfile, "%s\n", ctrstr); width = widths; for (ii = 0; ii < numcols - 1; ii++) { - memset(tmpstr, '-', *width); - tmpstr[*width++] = '\0'; - fprintf(obs->workfile, "%s--", tmpstr); + if (obs->numw==0 && ii==10) { // Skip jerk parts + width++; + continue; + } else { + memset(tmpstr, '-', *width); + tmpstr[*width++] = '\0'; + fprintf(obs->workfile, "%s--", tmpstr); + } } memset(tmpstr, '-', *width++); tmpstr[widths[ii]] = '\0'; @@ -704,6 +810,13 @@ void output_fundamentals(fourierprops * props, GSList * list, write_val_with_err(obs->workfile, errs.fd, errs.fderr, *error++, *width++); write_val_with_err(obs->workfile, props[ii].z, props[ii].zerr, *error++, *width++); + if (obs->numw) { + write_val_with_err(obs->workfile, props[ii].w, props[ii].werr, + *error++, *width++); + } else { + error++; + width++; + } accel = props[ii].z * SOL / (obs->T * obs->T * errs.f); accelerr = props[ii].zerr * SOL / (obs->T * obs->T * errs.f); write_val_with_err(obs->workfile, accel, accelerr, *error++, *width++); @@ -718,20 +831,20 @@ void output_fundamentals(fourierprops * props, GSList * list, void output_harmonics(GSList * list, accelobs * obs, infodata * idata) { - int ii, jj, numcols = 13, numcands; - int widths[13] = { 5, 4, 5, 15, 11, 18, 13, 12, 9, 12, 10, 10, 20 }; - int errors[13] = { 0, 0, 0, 2, 0, 2, 0, 2, 0, 2, 2, 2, 0 }; + int ii, jj, numcols = 15, numcands; + int widths[15] = { 5, 4, 5, 15, 11, 18, 13, 12, 9, 12, 9, 12, 10, 10, 20 }; + int errors[15] = { 0, 0, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 2, 2, 0 }; char tmpstr[30], ctrstr[30], notes[21], *command; accelcand *cand; GSList *listptr; fourierprops props; rzwerrs errs; static char *titles1[] = { "", "", "", "Power /", "Raw", - "FFT 'r'", "Pred 'r'", "FFT 'z'", "Pred 'z'", + "FFT 'r'", "Pred 'r'", "FFT 'z'", "Pred 'z'", "FFT 'w'", "Pred 'w'", "Phase", "Centroid", "Purity", "" }; static char *titles2[] = { "Cand", "Harm", "Sigma", "Loc Pow", "Power", - "(bin)", "(bin)", "(bins)", "(bins)", + "(bin)", "(bin)", "(bins)", "(bins)", "(bins)", "(bins)", "(rad)", "(0-1)", "

= 1", "Notes" }; @@ -741,30 +854,42 @@ void output_harmonics(GSList * list, accelobs * obs, infodata * idata) /* Print the header */ for (ii = 0; ii < numcols - 1; ii++) { - center_string(ctrstr, titles1[ii], widths[ii]); - fprintf(obs->workfile, "%s ", ctrstr); + if (obs->numw==0 && (ii==9 || ii==10)) { // Skip jerk parts + continue; + } else { + center_string(ctrstr, titles1[ii], widths[ii]); + fprintf(obs->workfile, "%s ", ctrstr); + } } center_string(ctrstr, titles1[ii], widths[ii]); fprintf(obs->workfile, "%s\n", ctrstr); for (ii = 0; ii < numcols - 1; ii++) { - if (obs->nph > 0.0 && ii == 3) /* HAAACK!!! */ - center_string(ctrstr, "NumPhot", widths[ii]); - else - center_string(ctrstr, titles2[ii], widths[ii]); - fprintf(obs->workfile, "%s ", ctrstr); + if (obs->numw==0 && (ii==9 || ii==10)) { // Skip jerk parts + continue; + } else { + if (obs->nph > 0.0 && ii == 3) /* HAAACK!!! */ + center_string(ctrstr, "NumPhot", widths[ii]); + else + center_string(ctrstr, titles2[ii], widths[ii]); + fprintf(obs->workfile, "%s ", ctrstr); + } } center_string(ctrstr, titles2[ii], widths[ii]); fprintf(obs->workfile, "%s\n", ctrstr); for (ii = 0; ii < numcols - 1; ii++) { - memset(tmpstr, '-', widths[ii]); - tmpstr[widths[ii]] = '\0'; - fprintf(obs->workfile, "%s--", tmpstr); + if (obs->numw==0 && (ii==9 || ii==10)) { // Skip jerk parts + continue; + } else { + memset(tmpstr, '-', widths[ii]); + tmpstr[widths[ii]] = '\0'; + fprintf(obs->workfile, "%s--", tmpstr); + } } memset(tmpstr, '-', widths[ii]); tmpstr[widths[ii]] = '\0'; fprintf(obs->workfile, "%s\n", tmpstr); - /* Print the fundamentals */ + /* Print the harmonics */ for (ii = 0; ii < numcands; ii++) { cand = (accelcand *) (listptr->data); @@ -775,11 +900,11 @@ void output_harmonics(GSList * list, accelobs * obs, infodata * idata) tmp_locpow = cand->derivs[jj].locpow; cand->derivs[jj].locpow = obs->nph; calc_props(cand->derivs[jj], cand->hirs[jj], - cand->hizs[jj], 0.0, &props); + cand->hizs[jj], cand->hiws[jj], &props); cand->derivs[jj].locpow = tmp_locpow; } else { calc_props(cand->derivs[jj], cand->hirs[jj], - cand->hizs[jj], 0.0, &props); + cand->hizs[jj], cand->hiws[jj], &props); } calc_rzwerrs(&props, obs->T, &errs); if (strncmp(idata->telescope, "None", 4) != 0) { @@ -815,12 +940,19 @@ void output_harmonics(GSList * list, accelobs * obs, infodata * idata) sprintf(tmpstr, "%.2f", cand->z * (jj + 1)); center_string(ctrstr, tmpstr, widths[8]); fprintf(obs->workfile, "%s ", ctrstr); + if (obs->numw) { + write_val_with_err(obs->workfile, props.w, props.werr, + errors[9], widths[9]); + sprintf(tmpstr, "%.2f", cand->w * (jj + 1)); + center_string(ctrstr, tmpstr, widths[10]); + fprintf(obs->workfile, "%s ", ctrstr); + } write_val_with_err(obs->workfile, props.phs, props.phserr, - errors[9], widths[9]); + errors[11], widths[11]); write_val_with_err(obs->workfile, props.cen, props.cenerr, - errors[10], widths[10]); + errors[12], widths[12]); write_val_with_err(obs->workfile, props.pur, props.purerr, - errors[11], widths[11]); + errors[13], widths[13]); fprintf(obs->workfile, " %.20s\n", notes); fflush(obs->workfile); } @@ -840,8 +972,8 @@ void print_accelcand(gpointer data, gpointer user_data) accelcand *obj = (accelcand *) data; user_data = NULL; - printf("sigma: %-7.4f pow: %-7.2f harm: %-2d r: %-14.4f z: %-10.4f\n", - obj->sigma, obj->power, obj->numharm, obj->r, obj->z); + printf("sigma: %-7.4f pow: %-7.2f harm: %-2d r: %-14.4f z: %-10.4f w: %-10.2f\n", + obj->sigma, obj->power, obj->numharm, obj->r, obj->z, obj->w); } @@ -876,65 +1008,63 @@ fcomplex *get_fourier_amplitudes(long long lobin, int numbins, accelobs * obs) } } -ffdotpows *subharm_ffdot_plane(int numharm, int harmnum, +ffdotpows *subharm_fderivs_vol(int numharm, int harmnum, double fullrlo, double fullrhi, subharminfo * shi, accelobs * obs) { - int ii, lobin, hibin, numdata, nice_numdata, fftlen, binoffset; + int ii, numdata, fftlen, binoffset; + long long lobin; float powargr, powargi; - static int numrs_full = 0; double drlo, drhi, harm_fract; - ffdotpows *ffdot; fcomplex *data, *pdata; fftwf_plan invplan; - - if (numrs_full == 0) { - if (numharm == 1 && harmnum == 1) { - numrs_full = ACCEL_USELEN; - } else { - printf("You must call subharm_ffdot_plane() with numharm=1 and\n"); - printf("harnum=1 before you use other values! Exiting.\n\n"); - exit(0); - } - } - ffdot = (ffdotpows *) malloc(sizeof(ffdotpows)); + ffdotpows *ffdot = (ffdotpows *) malloc(sizeof(ffdotpows)); /* Calculate and get the required amplitudes */ harm_fract = (double) harmnum / (double) numharm; drlo = calc_required_r(harm_fract, fullrlo); drhi = calc_required_r(harm_fract, fullrhi); - ffdot->rlo = (int) floor(drlo); + ffdot->rlo = (long long) floor(drlo); ffdot->zlo = calc_required_z(harm_fract, obs->zlo); + ffdot->wlo = calc_required_w(harm_fract, obs->wlo); /* Initialize the lookup indices */ if (numharm > 1 && !obs->inmem) { double rr, subr; - for (ii = 0; ii < numrs_full; ii++) { + for (ii = 0; ii < obs->corr_uselen; ii++) { rr = fullrlo + ii * ACCEL_DR; subr = calc_required_r(harm_fract, rr); shi->rinds[ii] = index_from_r(subr, ffdot->rlo); } + double zz, subz; + for (ii = 0; ii < obs->numz; ii++) { + zz = obs->zlo + ii * ACCEL_DZ; + subz = calc_required_z(harm_fract, zz); + shi->zinds[ii] = index_from_z(subz, ffdot->zlo); + } } ffdot->rinds = shi->rinds; + // The +1 below is important! ffdot->numrs = (int) ((ceil(drhi) - floor(drlo)) * ACCEL_RDR + DBLCORRECT) + 1; if (numharm == 1 && harmnum == 1) { - ffdot->numrs = ACCEL_USELEN; + ffdot->numrs = obs->corr_uselen; } else { - if (ffdot->numrs % ACCEL_RDR) { + if (ffdot->numrs % ACCEL_RDR) ffdot->numrs = (ffdot->numrs / ACCEL_RDR + 1) * ACCEL_RDR; - } } - ffdot->numzs = shi->numkern; - binoffset = shi->kern[0].kern_half_width; - fftlen = shi->kern[0].fftlen; + ffdot->zinds = shi->zinds; + ffdot->numzs = shi->numkern_zdim; + ffdot->numws = shi->numkern_wdim; + + /* Determine the largest kernel halfwidth needed to analyze the current subharmonic */ + /* Verified numerically that, as long as we have symmetric z's and w's, */ + /* shi->kern[0][0].kern_half_width is the maximal halfwidth over the range of w's and z's */ + binoffset = shi->kern[0][0].kern_half_width; + fftlen = shi->kern[0][0].fftlen; lobin = ffdot->rlo - binoffset; - hibin = (int) ceil(drhi) + binoffset; - numdata = hibin - lobin + 1; - nice_numdata = next2_to_n(numdata); // for FFTs - if (nice_numdata != fftlen / ACCEL_NUMBETWEEN) - printf("WARNING!!: nice_numdata != fftlen/2 in subharm_ffdot_plane()!\n"); - data = get_fourier_amplitudes(lobin, nice_numdata, obs); + numdata = fftlen / ACCEL_NUMBETWEEN; + data = get_fourier_amplitudes(lobin, numdata, obs); if (!obs->mmap_file && !obs->dat_input && 0) printf("This is newly malloc'd!\n"); @@ -947,7 +1077,7 @@ ffdotpows *subharm_ffdot_plane(int numharm, int harmnum, data[ii].i *= norm; } } else if (obs->norm_type == 0) { - // old-style block median normalization + // default block median normalization float *powers; double norm; powers = gen_fvect(numdata); @@ -960,13 +1090,13 @@ ffdotpows *subharm_ffdot_plane(int numharm, int harmnum, data[ii].i *= norm; } } else { - // new-style running double-tophat local-power normalization + // optional running double-tophat local-power normalization float *powers, *loc_powers; - powers = gen_fvect(nice_numdata); - for (ii = 0; ii < nice_numdata; ii++) { + powers = gen_fvect(numdata); + for (ii = 0; ii < numdata; ii++) { powers[ii] = POWER(data[ii].r, data[ii].i); } - loc_powers = corr_loc_pow(powers, nice_numdata); + loc_powers = corr_loc_pow(powers, numdata); for (ii = 0; ii < numdata; ii++) { float norm = invsqrtf(loc_powers[ii]); data[ii].r *= norm; @@ -983,7 +1113,7 @@ ffdotpows *subharm_ffdot_plane(int numharm, int harmnum, COMPLEXFFT(pdata, fftlen, -1); // Create the output power array - ffdot->powers = gen_fmatrix(ffdot->numzs, ffdot->numrs); + ffdot->powers = gen_f3Darr(ffdot->numws, ffdot->numzs, ffdot->numrs); // Create a plan with temp arrays. We will reuse the plan // with the new-array FFTW execute functions @@ -1009,41 +1139,46 @@ ffdotpows *subharm_ffdot_plane(int numharm, int harmnum, // tmpdat gets overwritten during the correlation fcomplex *tmpdat = gen_cvect(fftlen); fcomplex *tmpout = gen_cvect(fftlen); + int jj; #ifdef _OPENMP +// #pragma omp for collapse(2) Do we want this somehow? #pragma omp for #endif - for (ii = 0; ii < ffdot->numzs; ii++) { - int jj; - float *fkern = (float *) shi->kern[ii].data; - float *fpdata = (float *) pdata; - float *fdata = (float *) tmpdat; - float *outpows = ffdot->powers[ii]; - // multiply data and kernel - // (using floats for better vectorization) -#if (defined(__GNUC__) || defined(__GNUG__)) && \ + /* Check, should we add the collapse to parallelize numws and numzs loops? */ + for (ii = 0; ii < ffdot->numws; ii++) { + for (jj = 0; jj < ffdot->numzs; jj++) { + int kk; + float *fkern = (float *) shi->kern[ii][jj].data; + float *fpdata = (float *) pdata; + float *fdata = (float *) tmpdat; + float *outpows = ffdot->powers[ii][jj]; + // multiply data and kernel + // (using floats for better vectorization) +#if (defined(__GNUC__) || defined(__GNUG__)) && \ !(defined(__clang__) || defined(__INTEL_COMPILER)) #pragma GCC ivdep #endif - for (jj = 0; jj < fftlen * 2; jj += 2) { - const float dr = fpdata[jj], di = fpdata[jj + 1]; - const float kr = fkern[jj], ki = fkern[jj + 1]; - fdata[jj] = dr * kr + di * ki; - fdata[jj + 1] = di * kr - dr * ki; - } - // Do the inverse FFT (tmpdat -> tmpout) - fftwf_execute_dft(invplan, (fftwf_complex *) tmpdat, - (fftwf_complex *) tmpout); - // Turn the good parts of the result into powers and store - // them in the output matrix - fdata = (float *) tmpout; -#if (defined(__GNUC__) || defined(__GNUG__)) && \ + for (kk = 0; kk < fftlen * 2; kk += 2) { + const float dr = fpdata[kk], di = fpdata[kk + 1]; + const float kr = fkern[kk], ki = fkern[kk + 1]; + fdata[kk] = dr * kr + di * ki; + fdata[kk + 1] = di * kr - dr * ki; + } + // Do the inverse FFT (tmpdat -> tmpout) + fftwf_execute_dft(invplan, (fftwf_complex *) tmpdat, + (fftwf_complex *) tmpout); + // Turn the good parts of the result into powers and store + // them in the output matrix + fdata = (float *) tmpout; +#if (defined(__GNUC__) || defined(__GNUG__)) && \ !(defined(__clang__) || defined(__INTEL_COMPILER)) #pragma GCC ivdep #endif - for (jj = 0; jj < ffdot->numrs; jj++) { - const int ind = 2 * (jj + offset); - outpows[jj] = (fdata[ind] * fdata[ind] + - fdata[ind + 1] * fdata[ind + 1]) * norm; + for (kk = 0; kk < ffdot->numrs; kk++) { + const int ind = 2 * (kk + offset); + outpows[kk] = (fdata[ind] * fdata[ind] + + fdata[ind + 1] * fdata[ind + 1]) * norm; + } } } vect_free(tmpdat); @@ -1060,15 +1195,17 @@ ffdotpows *copy_ffdotpows(ffdotpows * orig) { int ii; ffdotpows *copy; - + copy = (ffdotpows *) malloc(sizeof(ffdotpows)); copy->numrs = orig->numrs; copy->numzs = orig->numzs; + copy->numws = orig->numws; copy->rlo = orig->rlo; copy->zlo = orig->zlo; - copy->powers = gen_fmatrix(orig->numzs, orig->numrs); - for (ii = 0; ii < (orig->numzs * orig->numrs); ii++) - copy->powers[0][ii] = orig->powers[0][ii]; + copy->wlo = orig->wlo; + copy->powers = gen_f3Darr(orig->numws, orig->numzs, orig->numrs); + for (ii = 0; ii < (orig->numws * orig->numzs * orig->numrs); ii++) + copy->powers[0][0][ii] = orig->powers[0][0][ii]; return copy; } @@ -1078,14 +1215,14 @@ void fund_to_ffdotplane(ffdotpows * ffd, accelobs * obs) // This moves the fundamental's ffdot plane powers // into the one for the full array int ii; - long long rlen = (obs->highestbin + ACCEL_USELEN) * ACCEL_RDR; + long long rlen = (obs->highestbin + obs->corr_uselen) * ACCEL_RDR; long long offset; float *outpow; for (ii = 0; ii < ffd->numzs; ii++) { offset = ii * rlen; outpow = obs->ffdotplane + offset + ffd->rlo * ACCEL_RDR; - memcpy(outpow, ffd->powers[ii], ffd->numrs * sizeof(float)); + memcpy(outpow, ffd->powers[0][ii], ffd->numrs * sizeof(float)); } } @@ -1098,10 +1235,8 @@ void fund_to_ffdotplane_trans(ffdotpows * ffd, accelobs * obs) // memory local (since numz << numr) int ii, jj; float *outpow = obs->ffdotplane + ffd->rlo * ACCEL_RDR * ffd->numzs; - //printf("rlo = %d, numrs = %d, nextrlo = %d\n", - // ffd->rlo, ffd->numrs, (int)(ffd->rlo+ffd->numrs*ACCEL_DR)); for (ii = 0; ii < ffd->numrs; ii++) { - float *inpow = ffd->powers[0] + ii; + float *inpow = ffd->powers[0][0] + ii; for (jj = 0; jj < ffd->numzs; jj++, inpow += ffd->numrs) { *outpow++ = *inpow; } @@ -1111,6 +1246,7 @@ void fund_to_ffdotplane_trans(ffdotpows * ffd, accelobs * obs) void free_ffdotpows(ffdotpows * ffd) { + vect_free(ffd->powers[0][0]); vect_free(ffd->powers[0]); vect_free(ffd->powers); free(ffd); @@ -1119,16 +1255,19 @@ void free_ffdotpows(ffdotpows * ffd) void add_ffdotpows(ffdotpows * fundamental, ffdotpows * subharmonic, int numharm, int harmnum) { - int ii, jj, zz, rind, zind, subz; + int ii, jj, kk, ww, rind, zind, wind, subw; const double harm_fract = (double) harmnum / (double) numharm; - - for (ii = 0; ii < fundamental->numzs; ii++) { - zz = fundamental->zlo + ii * ACCEL_DZ; - subz = calc_required_z(harm_fract, zz); - zind = index_from_z(subz, subharmonic->zlo); - for (jj = 0; jj < fundamental->numrs; jj++) { - rind = subharmonic->rinds[jj]; - fundamental->powers[ii][jj] += subharmonic->powers[zind][rind]; + + for (ii = 0; ii < fundamental->numws; ii++) { + ww = fundamental->wlo + ii * ACCEL_DW; + subw = calc_required_w(harm_fract, ww); + wind = index_from_w(subw, subharmonic->wlo); + for (jj = 0; jj < fundamental->numzs; jj++) { + zind = subharmonic->zinds[jj]; + for (kk = 0; kk < fundamental->numrs; kk++) { + rind = subharmonic->rinds[kk]; + fundamental->powers[ii][jj][kk] += subharmonic->powers[wind][zind][rind]; + } } } } @@ -1136,23 +1275,27 @@ void add_ffdotpows(ffdotpows * fundamental, void add_ffdotpows_ptrs(ffdotpows * fundamental, ffdotpows * subharmonic, int numharm, int harmnum) { - int ii, jj, zz, zind, subz; - const int zlo = fundamental->zlo; + int ii, jj, kk, ww, wind, subw; + const int wlo = fundamental->wlo; const int numrs = fundamental->numrs; const int numzs = fundamental->numzs; + const int numws = fundamental->numws; const double harm_fract = (double) harmnum / (double) numharm; float *outpows, *inpows; - unsigned short *indsptr; - - for (ii = 0; ii < numzs; ii++) { - zz = zlo + ii * ACCEL_DZ; - subz = calc_required_z(harm_fract, zz); - zind = index_from_z(subz, subharmonic->zlo); - inpows = subharmonic->powers[zind]; - outpows = fundamental->powers[ii]; - indsptr = subharmonic->rinds; - for (jj = 0; jj < numrs; jj++) - *outpows++ += inpows[*indsptr++]; + unsigned short *rindsptr, *zindsptr; + + for (ii = 0; ii < numws; ii++) { + ww = wlo + ii * ACCEL_DW; + subw = calc_required_w(harm_fract, ww); + wind = index_from_w(subw, subharmonic->wlo); + zindsptr = subharmonic->zinds; + for (jj = 0; jj < numzs; jj++) { + inpows = subharmonic->powers[wind][*zindsptr++]; + outpows = fundamental->powers[ii][jj]; + rindsptr = subharmonic->rinds; + for (kk = 0; kk < numrs; kk++) + *outpows++ += inpows[*rindsptr++]; + } } } @@ -1164,24 +1307,24 @@ void inmem_add_ffdotpows(ffdotpows * fundamental, accelobs * obs, const int numrs = fundamental->numrs; const int numzs = fundamental->numzs; const double harm_fract = (double) harmnum / (double) numharm; - int *indices; + int *rinds; // Pre-compute the frequency lookup table - indices = gen_ivect(numrs); + rinds = gen_ivect(numrs); { int ii, rrint; for (ii = 0, rrint = ACCEL_RDR * rlo; ii < numrs; ii++, rrint++) - indices[ii] = (int) (rrint * harm_fract + 0.5); + rinds[ii] = (int) (rrint * harm_fract + 0.5); } // Now add all the powers #ifdef _OPENMP -#pragma omp parallel default(none) shared(indices,fundamental,obs) +#pragma omp parallel default(none) shared(rinds,fundamental,obs) #endif { const int zlo = fundamental->zlo; - const long long rlen = (obs->highestbin + ACCEL_USELEN) * ACCEL_RDR; - float *powptr = fundamental->powers[0]; + const long long rlen = (obs->highestbin + obs->corr_uselen) * ACCEL_RDR; + float *powptr = fundamental->powers[0][0]; float *fdp = obs->ffdotplane; int ii, jj, zz, zind, subz; float *inpows, *outpows; @@ -1201,10 +1344,10 @@ void inmem_add_ffdotpows(ffdotpows * fundamental, accelobs * obs, #pragma GCC ivdep #endif for (jj = 0; jj < numrs; jj++) - outpows[jj] += inpows[indices[jj]]; + outpows[jj] += inpows[rinds[jj]]; } } - vect_free(indices); + vect_free(rinds); } @@ -1215,23 +1358,23 @@ void inmem_add_ffdotpows_trans(ffdotpows * fundamental, accelobs * obs, const int numrs = fundamental->numrs; const int numzs = fundamental->numzs; const double harm_fract = (double) harmnum / (double) numharm; - long *indices; + long *rinds; // Pre-compute the frequency lookup table - indices = gen_lvect(numrs); + rinds = gen_lvect(numrs); { int ii, rrint; for (ii = 0, rrint = ACCEL_RDR * rlo; ii < numrs; ii++, rrint++) - indices[ii] = (long) (rrint * harm_fract + 0.5) * numzs; + rinds[ii] = (long) (rrint * harm_fract + 0.5) * numzs; } // Now add all the powers #ifdef _OPENMP -#pragma omp parallel default(none) shared(indices,fundamental,obs) +#pragma omp parallel default(none) shared(rinds,fundamental,obs) #endif { const int zlo = fundamental->zlo; - float *powptr = fundamental->powers[0]; + float *powptr = fundamental->powers[0][0]; float *fdp = obs->ffdotplane; int ii, jj, zz, zind, subz; float *inpows, *outpows; @@ -1249,10 +1392,10 @@ void inmem_add_ffdotpows_trans(ffdotpows * fundamental, accelobs * obs, #pragma GCC ivdep #endif for (jj = 0; jj < numrs; jj++) - outpows[jj] += inpows[indices[jj]]; + outpows[jj] += inpows[rinds[jj]]; } } - vect_free(indices); + vect_free(rinds); } @@ -1262,36 +1405,40 @@ GSList *search_ffdotpows(ffdotpows * ffdot, int numharm, int ii; float powcut; long long numindep; - + powcut = obs->powcut[twon_to_index(numharm)]; numindep = obs->numindep[twon_to_index(numharm)]; - + #ifdef _OPENMP #pragma omp parallel for shared(ffdot,powcut,obs,numharm,numindep) #endif - for (ii = 0; ii < ffdot->numzs; ii++) { + for (ii = 0; ii < ffdot->numws; ii++) { int jj; - for (jj = 0; jj < ffdot->numrs; jj++) { - if (ffdot->powers[ii][jj] > powcut) { - float pow, sig; - double rr, zz; - int added = 0; - - pow = ffdot->powers[ii][jj]; - sig = candidate_sigma(pow, numharm, numindep); - rr = (ffdot->rlo + jj * (double) ACCEL_DR) / (double) numharm; - zz = (ffdot->zlo + ii * (double) ACCEL_DZ) / (double) numharm; + for (jj = 0; jj < ffdot->numzs; jj++) { + int kk; + for (kk = 0; kk < ffdot->numrs; kk++) { + if (ffdot->powers[ii][jj][kk] > powcut) { + float pow, sig; + double rr, zz, ww; + int added = 0; + + pow = ffdot->powers[ii][jj][kk]; + sig = candidate_sigma(pow, numharm, numindep); + rr = (ffdot->rlo + kk * (double) ACCEL_DR) / (double) numharm; + zz = (ffdot->zlo + jj * (double) ACCEL_DZ) / (double) numharm; + ww = (ffdot->wlo + ii * (double) ACCEL_DW) / (double) numharm; #ifdef _OPENMP #pragma omp critical #endif - { - cands = insert_new_accelcand(cands, pow, sig, numharm, - rr, zz, &added); + { + cands = insert_new_accelcand(cands, pow, sig, numharm, + rr, zz, ww, &added); + } + if (added && !obs->dat_input) + fprintf(obs->workfile, + "%-7.2f %-7.4f %-2d %-14.4f %-14.9f %-10.4f %-10.4f\n", + pow, sig, numharm, rr, rr / obs->T, zz, ww); } - if (added && !obs->dat_input) - fprintf(obs->workfile, - "%-7.2f %-7.4f %-2d %-14.4f %-14.9f %-10.4f\n", - pow, sig, numharm, rr, rr / obs->T, zz); } } } @@ -1483,16 +1630,6 @@ void create_accelobs(accelobs * obs, infodata * idata, Cmdline * cmd, int usemma printf("done.\n\n"); } - /* Determine the output filenames */ - - rootlen = strlen(obs->rootfilenm) + 25; - obs->candnm = (char *) calloc(rootlen, 1); - obs->accelnm = (char *) calloc(rootlen, 1); - obs->workfilenm = (char *) calloc(rootlen, 1); - sprintf(obs->candnm, "%s_ACCEL_%d.cand", obs->rootfilenm, cmd->zmax); - sprintf(obs->accelnm, "%s_ACCEL_%d", obs->rootfilenm, cmd->zmax); - sprintf(obs->workfilenm, "%s_ACCEL_%d.txtcand", obs->rootfilenm, cmd->zmax); - /* Open the FFT file if it exists appropriately */ if (!obs->dat_input) { obs->fftfile = chkfopen(cmd->argv[0], "rb"); @@ -1525,8 +1662,6 @@ void create_accelobs(accelobs * obs, infodata * idata, Cmdline * cmd, int usemma if (cmd->zmax % ACCEL_DZ) cmd->zmax = (cmd->zmax / ACCEL_DZ + 1) * ACCEL_DZ; - if (!obs->dat_input) - obs->workfile = chkfopen(obs->workfilenm, "w"); obs->N = (long long) idata->N; if (cmd->photonP) { if (obs->mmap_file || obs->dat_input) { @@ -1567,14 +1702,53 @@ void create_accelobs(accelobs * obs, infodata * idata, Cmdline * cmd, int usemma if (cmd->numharm != 1 && cmd->numharm != 2 && cmd->numharm != 4 && - cmd->numharm != 8 && cmd->numharm != 16 && cmd->numharm != 32) { + cmd->numharm != 8 && + cmd->numharm != 16 && + cmd->numharm != 32) { printf("\n'numharm' = %d must be a power-of-two! Exiting\n\n", cmd->numharm); exit(1); } obs->numharmstages = twon_to_index(cmd->numharm) + 1; + obs->dz = ACCEL_DZ; obs->numz = (cmd->zmax / ACCEL_DZ) * 2 + 1; + + /* Setting extra parameters for jerk search */ + if (cmd->wmaxP) { + if (cmd->wmax % ACCEL_DW) + cmd->wmax = (cmd->wmax / ACCEL_DW + 1) * ACCEL_DW; + obs->whi = cmd->wmax; + obs->wlo = -cmd->wmax; + obs->dw = ACCEL_DW; + obs->numw = (cmd->wmax / ACCEL_DW) * 2 + 1; + if (cmd->wmax==0.0) + obs->numw = 0; + printf("Jerk search enabled with maximum fdotdot wmax = %d\n", cmd->wmax); + } else { + obs->whi = 0.0; + obs->wlo = 0.0; + obs->dw = 0.0; + obs->numw = 0; + } + + /* Determine the output filenames */ + rootlen = strlen(obs->rootfilenm) + 45; + obs->candnm = (char *) calloc(rootlen, 1); + obs->accelnm = (char *) calloc(rootlen, 1); + obs->workfilenm = (char *) calloc(rootlen, 1); + if (obs->numw) { + sprintf(obs->candnm, "%s_ACCEL_%d_JERK_%d.cand", obs->rootfilenm, cmd->zmax, cmd->wmax); + sprintf(obs->accelnm, "%s_ACCEL_%d_JERK_%d", obs->rootfilenm, cmd->zmax, cmd->wmax); + sprintf(obs->workfilenm, "%s_ACCEL_%d_JERK_%d.txtcand", obs->rootfilenm, cmd->zmax, cmd->wmax); + } else { + sprintf(obs->candnm, "%s_ACCEL_%d.cand", obs->rootfilenm, cmd->zmax); + sprintf(obs->accelnm, "%s_ACCEL_%d", obs->rootfilenm, cmd->zmax); + sprintf(obs->workfilenm, "%s_ACCEL_%d.txtcand", obs->rootfilenm, cmd->zmax); + } + if (!obs->dat_input) + obs->workfile = chkfopen(obs->workfilenm, "w"); + obs->numbetween = ACCEL_NUMBETWEEN; obs->dt = idata->dt; obs->T = idata->dt * idata->N; @@ -1629,13 +1803,19 @@ void create_accelobs(accelobs * obs, infodata * idata, Cmdline * cmd, int usemma obs->powcut = (float *) malloc(obs->numharmstages * sizeof(float)); obs->numindep = (long long *) malloc(obs->numharmstages * sizeof(long long)); for (ii = 0; ii < obs->numharmstages; ii++) { - if (obs->numz == 1) + if (obs->numz == 1 && obs->numw == 0) obs->numindep[ii] = (obs->rhi - obs->rlo) / index_to_twon(ii); - else + else if (obs->numz > 1 && obs->numw == 0) /* The numz+1 takes care of the small amount of */ /* search we get above zmax and below zmin. */ obs->numindep[ii] = (obs->rhi - obs->rlo) * (obs->numz + 1) * (obs->dz / 6.95) / index_to_twon(ii); + else + /* The numw+1 takes care of the small amount of */ + /* search we get above wmax and below wmin. */ + obs->numindep[ii] = (obs->rhi - obs->rlo) * \ + (obs->numz + 1) * (obs->dz / 6.95) * \ + (obs->numw + 1) * (obs->dw / 44.2) / index_to_twon(ii); obs->powcut[ii] = power_for_sigma(obs->sigma, index_to_twon(ii), obs->numindep[ii]); } @@ -1648,17 +1828,35 @@ void create_accelobs(accelobs * obs, infodata * idata, Cmdline * cmd, int usemma obs->numzap = 0; */ + /* Determine corr_uselen from zmax and wmax */ + obs->maxkernlen = 2 * ACCEL_NUMBETWEEN * w_resp_halfwidth(obs->zhi, obs->whi, LOWACC); + obs->fftlen = fftlen_from_kernwidth(obs->maxkernlen); + if (obs->fftlen < 2048) + obs->fftlen = 2048; // This gives slightly better speed empirically + obs->corr_uselen = obs->fftlen - obs->maxkernlen; + // Make sure that obs->corr_uselen is an integer number of + // full (i.e. un-interpolated) Fourier bins + if (obs->corr_uselen % ACCEL_RDR) + obs->corr_uselen = obs->corr_uselen / ACCEL_RDR * ACCEL_RDR; + /* Can we perform the search in-core memory? */ { long long memuse; double gb = (double) (1L << 30); - // This is the size of powers covering the full f-dot plane to search - // Need the extra ACCEL_USELEN since we generate the plane in blocks - memuse = sizeof(float) * (obs->highestbin + ACCEL_USELEN) - * obs->numbetween * obs->numz; - printf("Full f-fdot plane would need %.2f GB: ", (float) memuse / gb); - if (memuse < MAXRAMUSE || cmd->inmemP) { + // This is the size of powers covering the full f-dot-dot plane to search + // Need the extra obs->corr_uselen since we generate the plane in blocks + if (cmd->wmaxP) { + memuse = sizeof(float) * (obs->highestbin + obs->corr_uselen) + * obs->numbetween * obs->numz * obs->numw; + printf("Full f-dot-dot volume would need %.2f GB: ", (float) memuse / gb); + } else { + memuse = sizeof(float) * (obs->highestbin + obs->corr_uselen) + * obs->numbetween * obs->numz; + printf("Full f-fdot plane would need %.2f GB: ", (float) memuse / gb); + } + + if (!cmd->wmaxP && (memuse < MAXRAMUSE || cmd->inmemP)) { printf("using in-memory accelsearch.\n\n"); obs->inmem = 1; obs->ffdotplane = gen_fvect(memuse / sizeof(float)); diff --git a/src/accelsearch.c b/src/accelsearch.c index 54679baae..d5692cb7d 100644 --- a/src/accelsearch.c +++ b/src/accelsearch.c @@ -42,7 +42,7 @@ static void print_percent_complete(int current, int number, char *what, int rese int main(int argc, char *argv[]) { - int ii; + int ii, rstep; double ttim, utim, stim, tott; struct tms runtimes; subharminfo **subharminfs; @@ -73,12 +73,15 @@ int main(int argc, char *argv[]) #endif printf("\n\n"); - printf(" Fourier-Domain Acceleration Search Routine\n"); - printf(" by Scott M. Ransom\n\n"); + printf(" Fourier-Domain Acceleration and Jerk Search Routine\n"); + printf(" by Scott M. Ransom\n\n"); /* Create the accelobs structure */ create_accelobs(&obs, &idata, cmd, 1); + /* The step-size of blocks to walk through the input data */ + rstep = obs.corr_uselen * ACCEL_DR; + /* Zap birdies if requested and if in memory */ if (cmd->zaplistP && !obs.mmap_file && obs.fft) { int numbirds; @@ -108,11 +111,13 @@ int main(int argc, char *argv[]) 1 << (obs.numharmstages - 1)); printf(" f = %.1f to %.1f Hz\n", obs.rlo / obs.T, obs.rhi / obs.T); printf(" r = %.1f to %.1f Fourier bins\n", obs.rlo, obs.rhi); - printf(" z = %.1f to %.1f Fourier bins drifted\n\n", obs.zlo, obs.zhi); + printf(" z = %.1f to %.1f Fourier bins drifted\n", obs.zlo, obs.zhi); + if (obs.numw) + printf(" w = %.1f to %.1f Fourier-derivative bins drifted\n", obs.wlo, obs.whi); /* Generate the correlation kernels */ - printf("Generating correlation kernels:\n"); + printf("\nGenerating correlation kernels:\n"); subharminfs = create_subharminfos(&obs); printf("Done generating kernels.\n\n"); if (cmd->ncpus > 1) { @@ -131,47 +136,61 @@ int main(int argc, char *argv[]) obs.workfilenm); } - /* Start the main search loop */ + /* Function pointers to make code a bit cleaner */ + void (*fund_to_ffdot)() = NULL; + void (*add_subharm)() = NULL; + void (*inmem_add_subharm)() = NULL; + if (obs.inmem) { + if (cmd->otheroptP) { + fund_to_ffdot = &fund_to_ffdotplane_trans; + inmem_add_subharm = &inmem_add_ffdotpows_trans; + } else { + fund_to_ffdot = &fund_to_ffdotplane; + inmem_add_subharm = &inmem_add_ffdotpows; + } + } else { + if (cmd->otheroptP) { + add_subharm = &add_ffdotpows_ptrs; + } else { + add_subharm = &add_ffdotpows; + } + } + /* Start the main search loop */ { double startr, lastr, nextr; ffdotpows *fundamental; /* Populate the saved F-Fdot plane at low freqs for in-memory * searches of harmonics that are below obs.rlo */ - if (obs.inmem) { startr = 8; // Choose a very low Fourier bin lastr = 0; nextr = 0; while (startr < obs.rlo) { - nextr = startr + ACCEL_USELEN * ACCEL_DR; + nextr = startr + rstep; lastr = nextr - ACCEL_DR; // Compute the F-Fdot plane - fundamental = subharm_ffdot_plane(1, 1, startr, lastr, + fundamental = subharm_fderivs_vol(1, 1, startr, lastr, &subharminfs[0][0], &obs); // Copy it into the full in-core one - if (cmd->otheroptP) - fund_to_ffdotplane_trans(fundamental, &obs); - else - fund_to_ffdotplane(fundamental, &obs); + fund_to_ffdot(fundamental, &obs); free_ffdotpows(fundamental); startr = nextr; } } /* Reset indices if needed and search for real */ - startr = obs.rlo; lastr = 0; nextr = 0; - while (startr + ACCEL_USELEN * ACCEL_DR < obs.highestbin) { + while (startr + rstep < obs.highestbin) { /* Search the fundamental */ print_percent_complete(startr - obs.rlo, obs.highestbin - obs.rlo, "search", 0); - nextr = startr + ACCEL_USELEN * ACCEL_DR; + nextr = startr + rstep; lastr = nextr - ACCEL_DR; - fundamental = subharm_ffdot_plane(1, 1, startr, lastr, + fundamental = subharm_fderivs_vol(1, 1, startr, lastr, &subharminfs[0][0], &obs); cands = search_ffdotpows(fundamental, 1, &obs, cands); @@ -180,33 +199,19 @@ int main(int argc, char *argv[]) ffdotpows *subharmonic; // Copy the fundamental's ffdot plane to the full in-core one - if (obs.inmem) { - if (cmd->otheroptP) - fund_to_ffdotplane_trans(fundamental, &obs); - else - fund_to_ffdotplane(fundamental, &obs); - } + if (obs.inmem) + fund_to_ffdot(fundamental, &obs); for (stage = 1; stage < obs.numharmstages; stage++) { harmtosum = 1 << stage; for (harm = 1; harm < harmtosum; harm += 2) { if (obs.inmem) { - if (cmd->otheroptP) - inmem_add_ffdotpows_trans(fundamental, &obs, - harmtosum, harm); - else - inmem_add_ffdotpows(fundamental, &obs, harmtosum, - harm); + inmem_add_subharm(fundamental, &obs, harmtosum, harm); } else { subharmonic = - subharm_ffdot_plane(harmtosum, harm, startr, lastr, + subharm_fderivs_vol(harmtosum, harm, startr, lastr, &subharminfs[stage][harm - 1], &obs); - if (cmd->otheroptP) - add_ffdotpows_ptrs(fundamental, subharmonic, - harmtosum, harm); - else - add_ffdotpows(fundamental, subharmonic, harmtosum, - harm); + add_subharm(fundamental, subharmonic, harmtosum, harm); free_ffdotpows(subharmonic); } } @@ -224,18 +229,14 @@ int main(int argc, char *argv[]) free_subharminfos(&obs, subharminfs); { /* Candidate list trimming and optimization */ - int numcands; + int numcands = g_slist_length(cands); GSList *listptr; accelcand *cand; fourierprops *props; - - numcands = g_slist_length(cands); - if (numcands) { /* Sort the candidates according to the optimized sigmas */ - cands = sort_accelcands(cands); /* Eliminate (most of) the harmonically related candidates */ @@ -243,7 +244,6 @@ int main(int argc, char *argv[]) eliminate_harmonics(cands, &numcands); /* Now optimize each candidate and its harmonics */ - print_percent_complete(0, 0, NULL, 1); listptr = cands; for (ii = 0; ii < numcands; ii++) { @@ -255,7 +255,6 @@ int main(int argc, char *argv[]) print_percent_complete(ii, numcands, "optimization", 0); /* Calculate the properties of the fundamentals */ - props = (fourierprops *) malloc(sizeof(fourierprops) * numcands); listptr = cands; for (ii = 0; ii < numcands; ii++) { @@ -264,23 +263,21 @@ int main(int argc, char *argv[]) /* send the originally determined r and z from the */ /* harmonic sum in the search. Note that the derivs are */ /* not used for the computations with the fundamental. */ - calc_props(cand->derivs[0], cand->r, cand->z, 0.0, props + ii); + calc_props(cand->derivs[0], cand->r, cand->z, cand->w, props + ii); /* Override the error estimates based on power */ props[ii].rerr = (float) (ACCEL_DR) / cand->numharm; props[ii].zerr = (float) (ACCEL_DZ) / cand->numharm; + props[ii].werr = (float) (ACCEL_DW) / cand->numharm; listptr = listptr->next; } /* Write the fundamentals to the output text file */ - output_fundamentals(props, cands, &obs, &idata); /* Write the harmonics to the output text file */ - output_harmonics(cands, &obs, &idata); - + /* Write the fundamental fourierprops to the cand file */ - obs.workfile = chkfopen(obs.candnm, "wb"); chkfwrite(props, sizeof(fourierprops), numcands, obs.workfile); fclose(obs.workfile); diff --git a/src/accelsearch_cmd.c b/src/accelsearch_cmd.c index 0aa62b3ff..9869f1b1a 100644 --- a/src/accelsearch_cmd.c +++ b/src/accelsearch_cmd.c @@ -38,6 +38,10 @@ static Cmdline cmd = { /* zmaxP = */ 1, /* zmax = */ 200, /* zmaxC = */ 1, + /***** -wmax: The max (+ and -) Fourier freq double derivs to search */ + /* wmaxP = */ 0, + /* wmax = */ (int)0, + /* wmaxC = */ 0, /***** -sigma: Cutoff sigma for choosing candidates */ /* sigmaP = */ 1, /* sigma = */ 2.0, @@ -831,6 +835,18 @@ showOptionValues(void) } } + /***** -wmax: The max (+ and -) Fourier freq double derivs to search */ + if( !cmd.wmaxP ) { + printf("-wmax not found.\n"); + } else { + printf("-wmax found:\n"); + if( !cmd.wmaxC ) { + printf(" no values\n"); + } else { + printf(" value = `%d'\n", cmd.wmax); + } + } + /***** -sigma: Cutoff sigma for choosing candidates */ if( !cmd.sigmaP ) { printf("-sigma not found.\n"); @@ -978,7 +994,7 @@ showOptionValues(void) void usage(void) { - fprintf(stderr,"%s"," [-ncpus ncpus] [-lobin lobin] [-numharm numharm] [-zmax zmax] [-sigma sigma] [-rlo rlo] [-rhi rhi] [-flo flo] [-fhi fhi] [-inmem] [-photon] [-median] [-locpow] [-zaplist zaplist] [-baryv baryv] [-otheropt] [-noharmpolish] [-noharmremove] [--] infile\n"); + fprintf(stderr,"%s"," [-ncpus ncpus] [-lobin lobin] [-numharm numharm] [-zmax zmax] [-wmax wmax] [-sigma sigma] [-rlo rlo] [-rhi rhi] [-flo flo] [-fhi fhi] [-inmem] [-photon] [-median] [-locpow] [-zaplist zaplist] [-baryv baryv] [-otheropt] [-noharmpolish] [-noharmremove] [--] infile\n"); fprintf(stderr,"%s"," Search an FFT or short time series for pulsars using a Fourier domain acceleration search with harmonic summing.\n"); fprintf(stderr,"%s"," -ncpus: Number of processors to use with OpenMP\n"); fprintf(stderr,"%s"," 1 int value between 1 and oo\n"); @@ -992,6 +1008,8 @@ usage(void) fprintf(stderr,"%s"," -zmax: The max (+ and -) Fourier freq deriv to search\n"); fprintf(stderr,"%s"," 1 int value between 0 and 1200\n"); fprintf(stderr,"%s"," default: `200'\n"); + fprintf(stderr,"%s"," -wmax: The max (+ and -) Fourier freq double derivs to search\n"); + fprintf(stderr,"%s"," 1 int value between 0 and 4000\n"); fprintf(stderr,"%s"," -sigma: Cutoff sigma for choosing candidates\n"); fprintf(stderr,"%s"," 1 float value between 1.0 and 30.0\n"); fprintf(stderr,"%s"," default: `2.0'\n"); @@ -1019,7 +1037,7 @@ usage(void) fprintf(stderr,"%s"," -noharmremove: Do not remove harmonically related candidates (never removed for numharm = 1)\n"); fprintf(stderr,"%s"," infile: Input file name of the floating point .fft or .[s]dat file. A '.inf' file of the same name must also exist\n"); fprintf(stderr,"%s"," 1 value\n"); - fprintf(stderr,"%s"," version: 07Dec16\n"); + fprintf(stderr,"%s"," version: 19Sep17\n"); fprintf(stderr,"%s"," "); exit(EXIT_FAILURE); } @@ -1075,6 +1093,16 @@ parseCmdline(int argc, char **argv) continue; } + if( 0==strcmp("-wmax", argv[i]) ) { + int keep = i; + cmd.wmaxP = 1; + i = getIntOpt(argc, argv, i, &cmd.wmax, 1); + cmd.wmaxC = i-keep; + checkIntLower("-wmax", &cmd.wmax, cmd.wmaxC, 4000); + checkIntHigher("-wmax", &cmd.wmax, cmd.wmaxC, 0); + continue; + } + if( 0==strcmp("-sigma", argv[i]) ) { int keep = i; cmd.sigmaP = 1; diff --git a/src/amoeba.c b/src/amoeba.c index eb5f0d4fe..4cad30aff 100644 --- a/src/amoeba.c +++ b/src/amoeba.c @@ -6,6 +6,7 @@ #include #include +#include "ransomfft.h" #ifndef SWAP /* Swaps two variables of undetermined type */ @@ -13,10 +14,13 @@ #endif static double amotry(double p[3][2], double *y, double *psum, - double (*funk) (double[]), int ihi, double fac); + double (*funk) (double[], fcomplex[], long[], float[], int[], int[]), + int ihi, double fac, fcomplex data[], long *numdata, + float *locpows, int *numharm, int *kernhw); void amoeba(double p[3][2], double *y, double ftol, - double (*funk) (double[]), int *nfunk) + double (*funk) (double[], fcomplex[], long[], float[], int[], int[]), + int *nfunk, fcomplex data[], long *numdata, float *locpows, int *numharm, int *kernhw) { int ii, ihi, ilo, inhi; double rtol, ysave, ytry, psum[2], tempzz; @@ -51,18 +55,22 @@ void amoeba(double p[3][2], double *y, double ftol, return; } *nfunk += 2; - ytry = amotry(p, y, psum, funk, ihi, -1.0); + ytry = amotry(p, y, psum, funk, ihi, -1.0, data, + numdata, locpows, numharm, kernhw); if (ytry <= y[ilo]) - ytry = amotry(p, y, psum, funk, ihi, 2.0); + ytry = amotry(p, y, psum, funk, ihi, 2.0, data, + numdata, locpows, numharm, kernhw); else if (ytry >= y[inhi]) { ysave = y[ihi]; - ytry = amotry(p, y, psum, funk, ihi, 0.5); + ytry = amotry(p, y, psum, funk, ihi, 0.5, data, + numdata, locpows, numharm, kernhw); if (ytry >= ysave) { for (ii = 0; ii <= 2; ii++) { if (ii != ilo) { p[ii][0] = psum[0] = 0.5 * (p[ii][0] + p[ilo][0]); p[ii][1] = psum[1] = 0.5 * (p[ii][1] + p[ilo][1]); - y[ii] = (*funk) (psum); + y[ii] = (*funk) (psum, data, numdata, locpows, + numharm, kernhw); } } *nfunk += 2; @@ -76,7 +84,9 @@ void amoeba(double p[3][2], double *y, double ftol, static double amotry(double p[3][2], double *y, double *psum, - double (*funk) (double[]), int ihi, double fac) + double (*funk) (double[], fcomplex[], long[], float[], int[], int[]), + int ihi, double fac, fcomplex data[], long *numdata, float *locpows, + int *numharm, int *kernhw) { double fac1, fac2, ytry, ptry[2]; @@ -84,7 +94,7 @@ static double amotry(double p[3][2], double *y, double *psum, fac2 = fac1 - fac; ptry[0] = psum[0] * fac1 - p[ihi][0] * fac2; ptry[1] = psum[1] * fac1 - p[ihi][1] * fac2; - ytry = (*funk) (ptry); + ytry = (*funk) (ptry, data, numdata, locpows, numharm, kernhw); if (ytry < y[ihi]) { y[ihi] = ytry; psum[0] += ptry[0] - p[ihi][0]; diff --git a/src/characteristics.c b/src/characteristics.c index d9b620324..4125447df 100644 --- a/src/characteristics.c +++ b/src/characteristics.c @@ -23,7 +23,7 @@ float get_numphotons(FILE * file) } -double get_localpower(fcomplex * data, int numdata, double r) +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. */ @@ -32,7 +32,7 @@ double get_localpower(fcomplex * data, int numdata, double r) /* interpolate. */ { double powargr, powargi, sum = 0.0; - int ii, binsperside, lo1, lo2, hi1, hi2, intfreq; + long ii, binsperside, lo1, lo2, hi1, hi2, intfreq; intfreq = (long) floor(r); binsperside = NUMLOCPOWAVG / 2; @@ -74,7 +74,7 @@ double get_localpower(fcomplex * data, int numdata, double r) } -double get_localpower3d(fcomplex * data, int numdata, double r, double z, double w) +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: */ @@ -136,7 +136,7 @@ double get_localpower3d(fcomplex * data, int numdata, double r, double z, double } -void get_derivs3d(fcomplex * data, int numdata, double r, +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 */ diff --git a/src/corr_prep.c b/src/corr_prep.c index 4c9328f85..d72f58b79 100644 --- a/src/corr_prep.c +++ b/src/corr_prep.c @@ -1,5 +1,39 @@ #include "presto.h" +int next_good_fftlen(int N) +/* Return one of the shortest, yet best performing, FFT lengths larger + * than N. This assumes FFTW. */ +{ + int fftlens[17] = {128, 192, 256, 384, 512, 768, 1024, 1280, 2048, 4096, + 5120, 7680, 10240, 12288, 15360, 16384, 25600}; + int ii = 0; + if (N <= fftlens[0]) + return fftlens[0]; + if (N > fftlens[16]) + return next2_to_n(N); + while (N > fftlens[ii]) ii++; + return fftlens[ii]; +} + +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. */ +{ + // The following nummbers were determined using FFTW 3.3.7 on an + // AVX-enabled processor. Metric used was max throughput of good + // correlated data. + if (kernwidth < 6) return 128; + else if (kernwidth < 52) return 256; + else if (kernwidth < 67) return 512; + else if (kernwidth < 378) return 1024; + else if (kernwidth < 664) return 2048; + else if (kernwidth < 1672) return 4096; + else if (kernwidth < 3015) return 10240; + else if (kernwidth < 3554) return 15360; + else if (kernwidth < 6000) return 25600; + else return next2_to_n(kernwidth*5); +} + void spread_with_pad(fcomplex * data, int numdata, fcomplex * result, int numresult, int numbetween, int numpad) /* Prepare the data array for correlation by spreading */ diff --git a/src/exploredat.c b/src/exploredat.c index 5519d6349..f3cdf27bb 100644 --- a/src/exploredat.c +++ b/src/exploredat.c @@ -43,27 +43,27 @@ static int plotstats = 0, usemedian = 0; typedef struct datapart { double tlo; /* Elapsed time (s) of the first point */ - int nlo; /* The sample number of the first point */ - int nn; /* The total number samples in *data */ - float *data; /* Raw data */ + long nlo; /* The sample number of the first point */ + long nn; /* The total number samples in *data */ + float *data; /* Raw data */ } datapart; typedef struct dataview { double vdt; /* Data view time step (2.0**(-zoomlevel))*dt */ - float minval; /* The minimum sample value in this view */ - float maxval; /* The maximum sample value in this view */ - int centern; /* The center sample to plot */ - int lon; /* The lowest sample to plot */ - int zoomlevel; /* Positive = zoomed in, Negative = zoomed out */ - int numsamps; /* The number of samples covered by the display */ - int dispnum; /* The number of points actually plotted */ - int chunklen; /* The length of the chunk of samples used to calculate stats */ - int numchunks; /* The number of chunks that are being displayed */ - float avgmeds[MAXDISPNUM]; /* The average or median samples for each chunk */ - float stds[MAXDISPNUM]; /* The standard deviation of the samples for each chunk */ - float maxs[MAXDISPNUM]; /* The maximum samples for each chunk */ - float mins[MAXDISPNUM]; /* The minimum samples for each chunk */ - float vals[MAXDISPNUM]; /* The raw data values when zoomlevel > 0 */ + float minval; /* The minimum sample value in this view */ + float maxval; /* The maximum sample value in this view */ + long centern; /* The center sample to plot */ + long lon; /* The lowest sample to plot */ + long zoomlevel; /* Positive = zoomed in, Negative = zoomed out */ + long numsamps; /* The number of samples covered by the display */ + long dispnum; /* The number of points actually plotted */ + long chunklen; /* The length of the chunk of samples used to calculate stats */ + long numchunks; /* The number of chunks that are being displayed */ + float avgmeds[MAXDISPNUM]; /* The average or median samples for each chunk */ + float stds[MAXDISPNUM]; /* The standard deviation of the samples for each chunk */ + float maxs[MAXDISPNUM]; /* The maximum samples for each chunk */ + float mins[MAXDISPNUM]; /* The minimum samples for each chunk */ + float vals[MAXDISPNUM]; /* The raw data values when zoomlevel > 0 */ } dataview; @@ -106,7 +106,7 @@ static basicstats *calc_stats(dataview * dv, datapart * dp) static int plot_dataview(dataview * dv, float minval, float maxval, float charhgt) /* The return value is offsetn */ { - int ii, lon, hin, offsetn = 0, tmpn; + long ii, lon, hin, offsetn = 0, tmpn; double lot, hit, offsett = 0.0; float ns[MAXDISPNUM], hiavg[MAXDISPNUM], loavg[MAXDISPNUM]; float scalemin = 0.0, scalemax = 0.0, dscale; @@ -165,7 +165,7 @@ static int plot_dataview(dataview * dv, float minval, float maxval, float charhg char label[50]; offsetn = (lon / 10000) * 10000; - numchar = snprintf(label, 50, "Sample - %d", offsetn); + numchar = snprintf(label, 50, "Sample - %ld", offsetn); cpgmtxt("B", 2.8, 0.5, 0.5, label); } else { cpgmtxt("B", 2.8, 0.5, 0.5, "Sample"); @@ -238,9 +238,9 @@ static int plot_dataview(dataview * dv, float minval, float maxval, float charhg } -static dataview *get_dataview(int centern, int zoomlevel, datapart * dp) +static dataview *get_dataview(long centern, long zoomlevel, datapart * dp) { - int ii, jj, offset; + long ii, jj, offset; double tmpavg, tmpvar; float *tmpchunk; dataview *dv; @@ -254,7 +254,7 @@ static dataview *get_dataview(int centern, int zoomlevel, datapart * dp) (1 << abs(zoomlevel)) : (1 << LOGMINCHUNKLEN); dv->dispnum = (dv->numsamps > MAXDISPNUM) ? MAXDISPNUM : dv->numsamps; if (DEBUGOUT) - printf("zoomlevel = %d numsamps = %d chunklen = %d dispnum %d nn = %d\n", + printf("zoomlevel = %ld numsamps = %ld chunklen = %ld dispnum %ld nn = %ld\n", dv->zoomlevel, dv->numsamps, dv->chunklen, dv->dispnum, dp->nn); dv->numchunks = dv->numsamps / dv->chunklen; dv->centern = centern; @@ -305,7 +305,7 @@ static dataview *get_dataview(int centern, int zoomlevel, datapart * dp) } -static datapart *get_datapart(int nlo, int numn) +static datapart *get_datapart(long nlo, long numn) { datapart *dp; @@ -365,8 +365,8 @@ static void print_help(void) int main(int argc, char *argv[]) { float minval = SMALLNUM, maxval = LARGENUM, inx = 0, iny = 0; - int centern, offsetn; - int zoomlevel, maxzoom = 0, minzoom, xid, psid; + long centern, offsetn; + long zoomlevel, maxzoom = 0, minzoom, xid, psid; char *rootfilenm, inchar; datapart *lodp; dataview *dv; @@ -430,11 +430,11 @@ int main(int argc, char *argv[]) lodp = get_datapart(0, Ndat); #else { - int numsamp; + long numsamp; datfile = chkfopen(argv[1], "rb"); Ndat = chkfilelen(datfile, sizeof(float)); - numsamp = (Ndat > MAXPTS) ? (int) MAXPTS : (int) Ndat; + numsamp = (Ndat > MAXPTS) ? (long) MAXPTS : (long) Ndat; lodp = get_datapart(0, numsamp); } #endif @@ -480,6 +480,7 @@ int main(int argc, char *argv[]) offsetn = plot_dataview(dv, minval, maxval, 1.0); break; case 'M': /* Toggle between median and average */ + /* FALLTHRU */ case 'm': usemedian = (usemedian) ? 0 : 1; free(dv); @@ -488,12 +489,15 @@ int main(int argc, char *argv[]) offsetn = plot_dataview(dv, minval, maxval, 1.0); break; case 'A': /* Zoom in */ + /* FALLTHRU */ case 'a': centern = inx + offsetn; + /* FALLTHRU */ case 'I': + /* FALLTHRU */ case 'i': if (DEBUGOUT) - printf(" Zooming in (zoomlevel = %d)...\n", zoomlevel); + printf(" Zooming in (zoomlevel = %ld)...\n", zoomlevel); if (zoomlevel < maxzoom) { zoomlevel++; free(dv); @@ -501,14 +505,17 @@ int main(int argc, char *argv[]) cpgpage(); offsetn = plot_dataview(dv, minval, maxval, 1.0); } else - printf(" Already at maximum zoom level (%d).\n", zoomlevel); + printf(" Already at maximum zoom level (%ld).\n", zoomlevel); break; case 'X': /* Zoom out */ + /* FALLTHRU */ case 'x': + /* FALLTHRU */ case 'O': + /* FALLTHRU */ case 'o': if (DEBUGOUT) - printf(" Zooming out (zoomlevel = %d)...\n", zoomlevel); + printf(" Zooming out (zoomlevel = %ld)...\n", zoomlevel); if (zoomlevel > minzoom) { zoomlevel--; free(dv); @@ -516,10 +523,11 @@ int main(int argc, char *argv[]) cpgpage(); offsetn = plot_dataview(dv, minval, maxval, 1.0); } else - printf(" Already at minimum zoom level (%d).\n", zoomlevel); + printf(" Already at minimum zoom level (%ld).\n", zoomlevel); break; case '<': /* Shift left 1 full screen */ centern -= dv->numsamps + dv->numsamps / 8; + /* FALLTHRU */ case ',': /* Shift left 1/8 screen */ if (DEBUGOUT) printf(" Shifting left...\n"); @@ -538,6 +546,7 @@ int main(int argc, char *argv[]) break; case '>': /* Shift right 1 full screen */ centern += dv->numsamps - dv->numsamps / 8; + /* FALLTHRU */ case '.': /* Shift right 1/8 screen */ centern += dv->numsamps / 8; if (DEBUGOUT) @@ -643,6 +652,7 @@ int main(int argc, char *argv[]) break; } case 'S': /* Auto-scale */ + /* FALLTHRU */ case 's': printf(" Auto-scaling is on.\n"); minval = SMALLNUM; @@ -651,6 +661,7 @@ int main(int argc, char *argv[]) offsetn = plot_dataview(dv, minval, maxval, 1.0); break; case 'G': /* Goto a time */ + /* FALLTHRU */ case 'g': { char timestr[50]; @@ -664,8 +675,8 @@ int main(int argc, char *argv[]) time = atof(timestr); } offsetn = 0.0; - centern = (int) (time / idata.dt + 0.5); - printf(" Moving to time %.15g (data point %d).\n", time, centern); + centern = (long) (time / idata.dt + 0.5); + printf(" Moving to time %.15g (data point %ld).\n", time, centern); free(dv); dv = get_dataview(centern, zoomlevel, lodp); cpgpage(); @@ -676,6 +687,7 @@ int main(int argc, char *argv[]) print_help(); break; case 'P': /* Print the current plot */ + /* FALLTHRU */ case 'p': { int len; @@ -701,11 +713,12 @@ int main(int argc, char *argv[]) } break; case 'V': /* Show the basic statistics for the current dataview */ + /* FALLTHRU */ case 'v': statvals = calc_stats(dv, lodp); printf("\n Statistics:\n" - " Low sample %d\n" - " Number of samples %d\n" + " Low sample %ld\n" + " Number of samples %ld\n" " Low time (s) %.7g\n" " Duration of samples (s) %.7g\n" " Maximum value %.7g\n" @@ -722,6 +735,7 @@ int main(int argc, char *argv[]) free(statvals); break; case 'Q': /* Quit */ + /* FALLTHRU */ case 'q': printf(" Quitting...\n"); free(dv); diff --git a/src/explorefft.c b/src/explorefft.c index 19a8d2e8d..d31c03b17 100644 --- a/src/explorefft.c +++ b/src/explorefft.c @@ -47,19 +47,19 @@ static int lenzaplist = 0; /* The number of possible lobin/hibin pairs in z static bird *zaplist = NULL; typedef struct fftpart { - int rlo; /* Lowest Fourier freq displayed */ - int numamps; /* Number of raw amplitudes */ + long rlo; /* Lowest Fourier freq displayed */ + long numamps; /* Number of raw amplitudes */ float maxrawpow; /* The highest raw power present */ float *rawpowers; /* The raw powers */ float *medians; /* The local median values (chunks of size LOCALCHUNK bins) */ float *normvals; /* The values to use for normalization (default is median/-log(0.5)) */ - fcomplex *amps; /* Raw FFT amplitudes */ + fcomplex *amps; /* Raw FFT amplitudes */ } fftpart; typedef struct fftview { double dr; /* Fourier frequency stepsize (2.0**(-zoomlevel)) */ double centerr; /* The center Fourier freq to plot */ - int lor; /* The lowest Fourier freq to plot */ + long lor; /* The lowest Fourier freq to plot */ int zoomlevel; /* Positive = zoomed in, Negative = zoomed out */ int numbins; /* The number of full bins from low to high to display */ float maxpow; /* The maximum normalized power in the view */ @@ -232,7 +232,7 @@ static fftview *get_fftview(double centerr, int zoomlevel, fftpart * fp) if (norm_const == 0.0) { for (ii = 0; ii < DISPLAYNUM; ii++) { fv->rs[ii] = fv->lor + ii * fv->dr; - index = (int) ((fv->rs[ii] - fp->rlo) * split + 0.5); + index = (long) ((fv->rs[ii] - fp->rlo) * split + 0.5); fv->powers[ii] = POWER(interp[ii].r, interp[ii].i) * fp->normvals[index]; } @@ -244,13 +244,13 @@ static fftview *get_fftview(double centerr, int zoomlevel, fftpart * fp) } vect_free(interp); } else { /* Down-sampled power spectrum */ - int jj, powindex, normindex, binstocombine; + long jj, powindex, normindex, binstocombine; float *tmprawpwrs, maxpow; binstocombine = (1 << abs(zoomlevel)); fv->numbins = DISPLAYNUM * binstocombine; fv->dr = (double) binstocombine; - fv->lor = (int) floor(centerr - 0.5 * fv->numbins); + fv->lor = (long) floor(centerr - 0.5 * fv->numbins); if (fv->lor + fv->numbins > Nfft) { fv->lor = Nfft - fv->numbins; fv->centerr = fv->lor + fv->numbins / 2; @@ -260,13 +260,13 @@ static fftview *get_fftview(double centerr, int zoomlevel, fftpart * fp) tmprawpwrs = gen_fvect(fv->numbins); if (norm_const == 0.0) { for (ii = 0; ii < fv->numbins; ii++) { - powindex = (int) (fv->lor - fp->rlo + ii + 0.5); - normindex = (int) (powindex * split); + powindex = (long) (fv->lor - fp->rlo + ii + 0.5); + normindex = (long) (powindex * split); tmprawpwrs[ii] = fp->rawpowers[powindex] * fp->normvals[normindex]; } } else { for (ii = 0; ii < fv->numbins; ii++) { - powindex = (int) (fv->lor - fp->rlo + ii + 0.5); + powindex = (long) (fv->lor - fp->rlo + ii + 0.5); tmprawpwrs[ii] = fp->rawpowers[powindex] * norm_const; } } @@ -288,7 +288,7 @@ static fftview *get_fftview(double centerr, int zoomlevel, fftpart * fp) } -static fftpart *get_fftpart(int rlo, int numr) +static fftpart *get_fftpart(long rlo, long numr) { int ii, jj, index; float powargr, powargi, tmppwr, chunk[LOCALCHUNK]; @@ -349,7 +349,7 @@ static void free_fftpart(fftpart * fp) static double find_peak(float inf, fftview * fv, fftpart * fp) { - int ii, lobin, hibin, maxbin = 0; + long ii, lobin, hibin, maxbin = 0; float maxpow = 0.0; double inr, viewfrac = 0.05, newmaxr, newmaxz; rderivs derivs; @@ -414,7 +414,7 @@ static fftview *get_harmonic(double rr, int zoomlevel, fftpart * fp) fftview *harmview; numharmbins = (1 << (LOGDISPLAYNUM - zoomlevel)); - harmpart = get_fftpart((int) (rr - numharmbins), 2 * numharmbins); + harmpart = get_fftpart((long) (rr - numharmbins), 2 * numharmbins); if (harmpart != NULL) { harmview = get_fftview(rr, zoomlevel, harmpart); free_fftpart(harmpart); @@ -504,7 +504,7 @@ static double harmonic_loop(int xid, double rr, int zoomlevel, fftpart * fp) cpgiden(); cpgscr(15, 0.8, 0.8, 0.8); numharmbins = (1 << (LOGDISPLAYNUM - zoomlevel)); - harmpart = get_fftpart((int) (rr - numharmbins), 2 * numharmbins); + harmpart = get_fftpart((long) (rr - numharmbins), 2 * numharmbins); harmview = get_fftview(rr, zoomlevel, harmpart); free_fftpart(harmpart); offsetf = plot_fftview(harmview, 0.0, 1.0, rr, 2); @@ -517,13 +517,13 @@ static double harmonic_loop(int xid, double rr, int zoomlevel, fftpart * fp) printf(" Wrote the plot to the file '%s'.\n", filename); } else if (choice == 'A' || choice == 'a') { if (iny > 1.0) - retval = rr * (int) (inx); + retval = rr * (long) (inx); else if (iny > 0.0) - retval = rr * ((int) (inx) + 4.0); + retval = rr * ((long) (inx) + 4.0); else if (iny > -1.0) - retval = rr / (int) (inx); + retval = rr / (long) (inx); else - retval = rr / ((int) (inx) + 4.0); + retval = rr / ((long) (inx) + 4.0); badchoice = 0; } else { printf(" Option not recognized.\n"); @@ -603,11 +603,11 @@ int main(int argc, char *argv[]) lofp = get_fftpart(0, Nfft); #else { - int numamps; + long numamps; fftfile = chkfopen(argv[1], "rb"); Nfft = chkfilelen(fftfile, sizeof(fcomplex)); - numamps = (Nfft > MAXBINS) ? (int) MAXBINS : (int) Nfft; + numamps = (Nfft > MAXBINS) ? (long) MAXBINS : (long) Nfft; lofp = get_fftpart(0, numamps); } #endif @@ -615,7 +615,7 @@ int main(int argc, char *argv[]) /* Plot the initial data */ { - int initnumbins = INITIALNUMBINS; + long initnumbins = INITIALNUMBINS; if (initnumbins > Nfft) { initnumbins = next2_to_n(Nfft) / 2; @@ -655,9 +655,12 @@ int main(int argc, char *argv[]) switch (inchar) { case 'A': /* Zoom in */ + /* FALLTHRU */ case 'a': centerr = (inx + offsetf) * T; + /* FALLTHRU */ case 'I': + /* FALLTHRU */ case 'i': if (DEBUGOUT) printf(" Zooming in (zoomlevel = %d)...\n", zoomlevel); @@ -671,8 +674,11 @@ int main(int argc, char *argv[]) printf(" Already at maximum zoom level (%d).\n", zoomlevel); break; case 'X': /* Zoom out */ + /* FALLTHRU */ case 'x': + /* FALLTHRU */ case 'O': + /* FALLTHRU */ case 'o': if (DEBUGOUT) printf(" Zooming out (zoomlevel = %d)...\n", zoomlevel); @@ -687,6 +693,7 @@ int main(int argc, char *argv[]) break; case '<': /* Shift left 1 full screen */ centerr -= fv->numbins + fv->numbins / 8; + /* FALLTHRU */ case ',': /* Shift left 1/8 screen */ if (DEBUGOUT) printf(" Shifting left...\n"); @@ -705,6 +712,7 @@ int main(int argc, char *argv[]) break; case '>': /* Shift right 1 full screen */ centerr += fv->numbins - fv->numbins / 8; + /* FALLTHRU */ case '.': /* Shift right 1/8 screen */ if (DEBUGOUT) printf(" Shifting right...\n"); @@ -722,6 +730,7 @@ int main(int argc, char *argv[]) offsetf = plot_fftview(fv, maxpow, 1.0, 0.0, 0); break; case '+': /* Increase height of powers */ + /* FALLTHRU */ case '=': if (maxpow == 0.0) { printf(" Auto-scaling is off.\n"); @@ -732,6 +741,7 @@ int main(int argc, char *argv[]) offsetf = plot_fftview(fv, maxpow, 1.0, 0.0, 0); break; case '-': /* Decrease height of powers */ + /* FALLTHRU */ case '_': if (maxpow == 0.0) { printf(" Auto-scaling is off.\n"); @@ -742,6 +752,7 @@ int main(int argc, char *argv[]) offsetf = plot_fftview(fv, maxpow, 1.0, 0.0, 0); break; case 'S': /* Auto-scale */ + /* FALLTHRU */ case 's': if (maxpow == 0.0) break; @@ -753,6 +764,7 @@ int main(int argc, char *argv[]) break; } case 'G': /* Goto a frequency */ + /* FALLTHRU */ case 'g': { char freqstr[50]; @@ -774,6 +786,7 @@ int main(int argc, char *argv[]) } break; case 'H': /* Show harmonics */ + /* FALLTHRU */ case 'h': { double retval; @@ -793,6 +806,7 @@ int main(int argc, char *argv[]) print_help(); break; case 'D': /* Show details about a selected point */ + /* FALLTHRU */ case 'd': { double newr; @@ -808,6 +822,7 @@ int main(int argc, char *argv[]) } break; case 'L': /* Load a zaplist */ + /* FALLTHRU */ case 'l': { int ii, len; @@ -835,6 +850,7 @@ int main(int argc, char *argv[]) } break; case 'Z': /* Add a birdie to a zaplist */ + /* FALLTHRU */ case 'z': { int badchoice = 2; @@ -884,6 +900,7 @@ int main(int argc, char *argv[]) } break; case 'P': /* Print the current plot */ + /* FALLTHRU */ case 'p': { int len; @@ -910,6 +927,7 @@ int main(int argc, char *argv[]) } break; case 'N': /* Changing power normalization */ + /* FALLTHRU */ case 'n': { float inx2 = 0.0, iny2 = 0.0; @@ -925,6 +943,7 @@ int main(int argc, char *argv[]) cpgcurs(&inx2, &iny2, &choice); switch (choice) { case 'M': + /* FALLTHRU */ case 'm': norm_const = 0.0; maxpow = 0.0; @@ -933,6 +952,7 @@ int main(int argc, char *argv[]) (" Using local median normalization. Autoscaling is on.\n"); break; case 'D': + /* FALLTHRU */ case 'd': norm_const = 1.0 / r0; maxpow = 0.0; @@ -942,6 +962,7 @@ int main(int argc, char *argv[]) r0); break; case 'R': + /* FALLTHRU */ case 'r': norm_const = 1.0; maxpow = 0.0; @@ -950,11 +971,12 @@ int main(int argc, char *argv[]) (" Using raw powers (i.e. no normalization). Autoscaling is on.\n"); break; case 'U': + /* FALLTHRU */ case 'u': { char choice2; float xx = inx, yy = iny; - int lor, hir, numr; + long lor, hir, numr; double avg, var; printf @@ -963,19 +985,19 @@ int main(int argc, char *argv[]) do { cpgcurs(&xx, &yy, &choice2); } while (choice2 != 'A' && choice2 != 'a'); - lor = (int) ((xx + offsetf) * T); + lor = (long) ((xx + offsetf) * T); cpgsci(7); cpgmove(xx, 0.0); cpgdraw(xx, 10.0 * fv->maxpow); do { cpgcurs(&xx, &yy, &choice2); } while (choice2 != 'A' && choice2 != 'a'); - hir = (int) ((xx + offsetf) * T); + hir = (long) ((xx + offsetf) * T); cpgmove(xx, 0.0); cpgdraw(xx, 10.0 * fv->maxpow); cpgsci(1); if (lor > hir) { - int tempr; + long tempr; tempr = hir; hir = lor; lor = tempr; @@ -1006,6 +1028,7 @@ int main(int argc, char *argv[]) } break; case 'Q': /* Quit */ + /* FALLTHRU */ case 'q': printf(" Quitting...\n"); free(fv); diff --git a/src/makewisdom.c b/src/makewisdom.c index 7d1fe008b..5168101c6 100644 --- a/src/makewisdom.c +++ b/src/makewisdom.c @@ -14,8 +14,8 @@ int main(void) fftwf_plan plan; fftwf_complex *inout; int ii, fftlen; - int padlen[13] = { 288, 540, 1080, 2100, 4200, 8232, 16464, 32805, - 65610, 131220, 262440, 525000, 1050000 + int padlen[20] = { 192, 288, 384, 540, 768, 1080, 1280, 2100, 4200, 5120, + 7680, 8232, 10240, 12288, 15360, 16464, 25600, 32805, 65610, 131220 }; fftlen = 2; @@ -34,7 +34,7 @@ int main(void) printf("Generating plans for FFTs of length:\n"); inout = fftwf_malloc(sizeof(fftwf_complex) * BIGFFTWSIZE + 2); - while (fftlen <= 1.1e6) { + while (fftlen <= 1.1e5) { printf(" %d forward\n", fftlen); plan = fftwf_plan_dft_1d(fftlen, inout, inout, FFTW_FORWARD, FFTW_PATIENT); fftwf_destroy_plan(plan); @@ -50,7 +50,7 @@ int main(void) fftlen = 10; - while (fftlen <= 1.1e6) { + while (fftlen <= 1.1e5) { inout = fftwf_malloc(sizeof(fftwf_complex) * fftlen); printf(" %d forward\n", fftlen); plan = fftwf_plan_dft_1d(fftlen, inout, inout, FFTW_FORWARD, FFTW_PATIENT); diff --git a/src/maximize_r.c b/src/maximize_r.c index 75d749400..2a484ec27 100644 --- a/src/maximize_r.c +++ b/src/maximize_r.c @@ -9,7 +9,8 @@ double fminbr(double a, double b, double (*f) (double x), double tol); /* Some Static-Global Variables */ static fcomplex *maxdata; -static int nummaxdata, max_kern_half_width; +static int max_kern_half_width; +static long nummaxdata; static double power_call_r(double r) /* Maximization function used with an array */ @@ -24,7 +25,7 @@ static double power_call_r(double r) } -double max_r_arr(fcomplex * data, int numdata, double rin, +double max_r_arr(fcomplex * data, long numdata, double rin, double *rout, rderivs * derivs) /* Return the Fourier frequency that maximizes the power. */ { diff --git a/src/maximize_rz.c b/src/maximize_rz.c index aa2da89d5..a65a3fa07 100644 --- a/src/maximize_rz.c +++ b/src/maximize_rz.c @@ -2,34 +2,33 @@ #define ZSCALE 4.0 -static fcomplex *maxdata; -static int nummaxdata, max_kern_half_width; - -extern void amoeba(double p[3][2], double y[], double ftol, - double (*funk) (double[]), int *nfunk); - -static double power_call_rz(double rz[]) +#define UNUSED(x) (void)(x) + +void amoeba(double p[3][2], double *y, double ftol, + double (*funk) (double[], fcomplex[], long[], float[], int[], int[]), + int *nfunk, fcomplex data[], long *numdata, float *locpows, int *numharm, int *kernhw); + +static double power_call_rz(double rz[], fcomplex data[], long *numdata, + float *locpows, int *numharm, int *kernhw) /* f-fdot plane power function */ { double powargr, powargi; fcomplex ans; - /* num_funct_calls++; */ - rz_interp(maxdata, nummaxdata, rz[0], rz[1] * ZSCALE, max_kern_half_width, &ans); + UNUSED(locpows); + UNUSED(numharm); + rz_interp(data, *numdata, rz[0], rz[1] * ZSCALE, *kernhw, &ans); return -POWER(ans.r, ans.i); } -double max_rz_arr(fcomplex * data, int numdata, double rin, double zin, +double max_rz_arr(fcomplex * data, long numdata, double rin, double zin, double *rout, double *zout, rderivs * derivs) /* Return the Fourier frequency and Fourier f-dot that */ /* maximizes the power. */ { double y[3], x[3][2], step = 0.4; - float locpow; - int numeval; - - maxdata = data; - nummaxdata = numdata; + float locpow = 0.0; + int numeval = 0, numharm = 1, max_kernhw; /* Now prep the maximization at LOWACC for speed */ @@ -37,7 +36,7 @@ double max_rz_arr(fcomplex * data, int numdata, double rin, double zin, /* the true value of z is a little larger than z. This */ /* keeps a little more accuracy. */ - max_kern_half_width = z_resp_halfwidth(fabs(zin) + 4.0, LOWACC); + max_kernhw = z_resp_halfwidth(fabs(zin) + 4.0, LOWACC); /* Initialize the starting simplex */ @@ -50,18 +49,19 @@ double max_rz_arr(fcomplex * data, int numdata, double rin, double zin, /* Initialize the starting function values */ - y[0] = power_call_rz(x[0]); - y[1] = power_call_rz(x[1]); - y[2] = power_call_rz(x[2]); + y[0] = power_call_rz(x[0], data, &numdata, &locpow, &numharm, &max_kernhw); + y[1] = power_call_rz(x[1], data, &numdata, &locpow, &numharm, &max_kernhw); + y[2] = power_call_rz(x[2], data, &numdata, &locpow, &numharm, &max_kernhw); /* Call the solver: */ numeval = 0; - amoeba(x, y, 1.0e-7, power_call_rz, &numeval); + amoeba(x, y, 1.0e-7, power_call_rz, &numeval, + data, &numdata, &locpow, &numharm, &max_kernhw); /* Restart at minimum using HIGHACC to get a better result */ - max_kern_half_width = z_resp_halfwidth(fabs(x[0][1]) + 4.0, HIGHACC); + max_kernhw = z_resp_halfwidth(fabs(x[0][1]) + 4.0, HIGHACC); /* Re-Initialize some of the starting simplex */ @@ -72,14 +72,15 @@ double max_rz_arr(fcomplex * data, int numdata, double rin, double zin, /* Re-Initialize the starting function values */ - y[0] = power_call_rz(x[0]); - y[1] = power_call_rz(x[1]); - y[2] = power_call_rz(x[2]); + y[0] = power_call_rz(x[0], data, &numdata, &locpow, &numharm, &max_kernhw); + y[1] = power_call_rz(x[1], data, &numdata, &locpow, &numharm, &max_kernhw); + y[2] = power_call_rz(x[2], data, &numdata, &locpow, &numharm, &max_kernhw); /* Call the solver: */ numeval = 0; - amoeba(x, y, 1.0e-10, power_call_rz, &numeval); + amoeba(x, y, 1.0e-10, power_call_rz, &numeval, + data, &numdata, &locpow, &numharm, &max_kernhw); /* The following calculates derivatives at the peak */ @@ -97,14 +98,15 @@ double max_rz_file(FILE * fftfile, double rin, double zin, /* maximizes the power of the candidate in 'fftfile'. */ { double maxz, maxpow, rin_int, rin_frac; - int kern_half_width, filedatalen, startbin, extra = 10; + int kern_half_width, filedatalen, extra = 10; + long startbin; fcomplex *filedata; maxz = fabs(zin) + 4.0; rin_frac = modf(rin, &rin_int); kern_half_width = z_resp_halfwidth(maxz, HIGHACC); filedatalen = 2 * kern_half_width + extra; - startbin = (int) rin_int - filedatalen / 2; + startbin = (long) rin_int - filedatalen / 2; filedata = read_fcomplex_file(fftfile, startbin, filedatalen); maxpow = max_rz_arr(filedata, filedatalen, rin_frac + filedatalen / 2, @@ -113,58 +115,41 @@ double max_rz_file(FILE * fftfile, double rin, double zin, vect_free(filedata); return maxpow; } - - -static int max_num_harmonics; -static fcomplex **maxdata_harmonics; -static float *maxlocpow; -static int *maxr_offset; -static double power_call_rz_harmonics(double rz[]) +static double power_call_rz_harmonics(double rz[], fcomplex data[], long *numdata, + float *locpows, int *numharm, int *kernhw) { - int i; + int ii; double total_power = 0.; double powargr, powargi; fcomplex ans; - for (i = 1; i <= max_num_harmonics; i++) { - rz_interp(maxdata_harmonics[i - 1], nummaxdata, - (maxr_offset[i - 1] + rz[0]) * i - maxr_offset[i - 1], - rz[1] * ZSCALE * i, max_kern_half_width, &ans); - total_power += POWER(ans.r, ans.i) / maxlocpow[i - 1]; + for (ii = 0; ii < *numharm; ii++) { + int n = ii + 1; + rz_interp(data, *numdata, rz[0] * n, rz[1] * ZSCALE * n, *kernhw, &ans); + total_power += POWER(ans.r, ans.i) / locpows[ii]; } return -total_power; } -void max_rz_arr_harmonics(fcomplex * data[], int num_harmonics, - int r_offset[], - int numdata, double rin, double zin, - double *rout, double *zout, rderivs derivs[], - double power[]) +void max_rz_arr_harmonics(fcomplex data[], long numdata, + int num_harmonics, + double rin, double zin, + double *rout, double *zout, + rderivs derivs[], double powers[]) /* Return the Fourier frequency and Fourier f-dot that */ /* maximizes the power. */ { double y[3], x[3][2], step = 0.4; float *locpow; - int numeval; - int i; + int numeval, ii, max_kernhw; locpow = gen_fvect(num_harmonics); - maxlocpow = gen_fvect(num_harmonics); - maxr_offset = r_offset; - maxdata_harmonics = data; - - - //FIXME: z needs to be multiplied by i everywhere - for (i = 1; i <= num_harmonics; i++) { - locpow[i - 1] = - get_localpower3d(data[i - 1], numdata, - (r_offset[i - 1] + rin) * i - r_offset[i - 1], zin * i, - 0.0); - maxlocpow[i - 1] = locpow[i - 1]; + + for (ii = 0; ii < num_harmonics; ii++) { + int n = ii + 1; + locpow[ii] = get_localpower3d(data, numdata, rin * n, zin * n, 0.0); } - nummaxdata = numdata; - max_num_harmonics = num_harmonics; /* Now prep the maximization at LOWACC for speed */ @@ -172,7 +157,7 @@ void max_rz_arr_harmonics(fcomplex * data[], int num_harmonics, /* the true value of z is a little larger than z. This */ /* keeps a little more accuracy. */ - max_kern_half_width = z_resp_halfwidth(fabs(zin * num_harmonics) + 4.0, LOWACC); + max_kernhw = z_resp_halfwidth(fabs(zin * num_harmonics) + 4.0, LOWACC); /* Initialize the starting simplex */ @@ -185,19 +170,20 @@ void max_rz_arr_harmonics(fcomplex * data[], int num_harmonics, /* Initialize the starting function values */ - y[0] = power_call_rz_harmonics(x[0]); - y[1] = power_call_rz_harmonics(x[1]); - y[2] = power_call_rz_harmonics(x[2]); + y[0] = power_call_rz_harmonics(x[0], data, &numdata, locpow, &num_harmonics, &max_kernhw); + y[1] = power_call_rz_harmonics(x[1], data, &numdata, locpow, &num_harmonics, &max_kernhw); + y[2] = power_call_rz_harmonics(x[2], data, &numdata, locpow, &num_harmonics, &max_kernhw); /* Call the solver: */ numeval = 0; - amoeba(x, y, 1.0e-7, power_call_rz_harmonics, &numeval); + amoeba(x, y, 1.0e-7, power_call_rz_harmonics, &numeval, + data, &numdata, locpow, &num_harmonics, &max_kernhw); /* Restart at minimum using HIGHACC to get a better result */ - max_kern_half_width = - z_resp_halfwidth(fabs(x[0][1] * num_harmonics) + 4.0, HIGHACC); + max_kernhw = z_resp_halfwidth(fabs(x[0][1] * num_harmonics) + + 4.0, HIGHACC); /* Re-Initialize some of the starting simplex */ @@ -208,73 +194,27 @@ void max_rz_arr_harmonics(fcomplex * data[], int num_harmonics, /* Re-Initialize the starting function values */ - y[0] = power_call_rz_harmonics(x[0]); - y[1] = power_call_rz_harmonics(x[1]); - y[2] = power_call_rz_harmonics(x[2]); + y[0] = power_call_rz_harmonics(x[0], data, &numdata, locpow, &num_harmonics, &max_kernhw); + y[1] = power_call_rz_harmonics(x[1], data, &numdata, locpow, &num_harmonics, &max_kernhw); + y[2] = power_call_rz_harmonics(x[2], data, &numdata, locpow, &num_harmonics, &max_kernhw); /* Call the solver: */ numeval = 0; - amoeba(x, y, 1.0e-10, power_call_rz_harmonics, &numeval); + amoeba(x, y, 1.0e-10, power_call_rz_harmonics, &numeval, + data, &numdata, locpow, &num_harmonics, &max_kernhw); /* The following calculates derivatives at the peak */ *rout = x[0][0]; *zout = x[0][1] * ZSCALE; - for (i = 1; i <= num_harmonics; i++) { - locpow[i - 1] = - get_localpower3d(data[i - 1], numdata, - (r_offset[i - 1] + *rout) * i - r_offset[i - 1], - (*zout) * i, 0.0); - x[0][0] = (r_offset[i - 1] + *rout) * i - r_offset[i - 1]; - x[0][1] = *zout / ZSCALE * i; - maxdata = data[i - 1]; - power[i - 1] = -power_call_rz(x[0]); - get_derivs3d(data[i - 1], numdata, - (r_offset[i - 1] + *rout) * i - r_offset[i - 1], (*zout) * i, - 0.0, locpow[i - 1], &(derivs[i - 1])); + for (ii = 0; ii < num_harmonics; ii++) { + int n = ii + 1; + locpow[ii] = get_localpower3d(data, numdata, *rout * n, *zout * n, 0.0); + x[0][0] = *rout * n; + x[0][1] = *zout / ZSCALE * n; + powers[ii] = -power_call_rz(x[0], data, &numdata, locpow, &num_harmonics, &max_kernhw); + get_derivs3d(data, numdata, *rout * n, *zout * n, 0.0, locpow[ii], &(derivs[ii])); } - vect_free(locpow); - vect_free(maxlocpow); -} - -void max_rz_file_harmonics(FILE * fftfile, int num_harmonics, - int lobin, - double rin, double zin, - double *rout, double *zout, rderivs derivs[], - double maxpow[]) -/* Return the Fourier frequency and Fourier f-dot that */ -/* maximizes the power of the candidate in 'fftfile'. */ -/* WARNING: not tested */ -{ - int i; - double maxz, rin_int, rin_frac; - int kern_half_width, filedatalen, extra = 10; - int *r_offset; - fcomplex **filedata; - - r_offset = (int *) malloc(sizeof(int) * num_harmonics); - filedata = (fcomplex **) malloc(sizeof(fcomplex *) * num_harmonics); - maxz = fabs(zin * num_harmonics) + 4.0; - kern_half_width = z_resp_halfwidth(maxz, HIGHACC); - filedatalen = 2 * kern_half_width + extra; - - for (i = 1; i <= num_harmonics; i++) { - rin_frac = modf(rin * i, &rin_int); - r_offset[i - 1] = (int) rin_int - filedatalen / 2 + lobin; - filedata[i - 1] = read_fcomplex_file(fftfile, r_offset[i - 1], filedatalen); - } - rin_frac = modf(rin, &rin_int); - max_rz_arr_harmonics(filedata, num_harmonics, - r_offset, - filedatalen, rin_frac + filedatalen / 2, - zin, rout, zout, derivs, maxpow); - - *rout += r_offset[0]; - for (i = 1; i <= num_harmonics; i++) { - vect_free(filedata[i - 1]); - } - free(r_offset); - free(filedata); } diff --git a/src/maximize_rzw.c b/src/maximize_rzw.c index 15c4092f7..97af5dda0 100644 --- a/src/maximize_rzw.c +++ b/src/maximize_rzw.c @@ -1,9 +1,6 @@ #include "accel.h" -// Will need to zap this when we merge with Bridget's code -#define ACCEL_DW 20.0 - -double max_rzw_arr(fcomplex * data, int numdata, double rin, double zin, +double max_rzw_arr(fcomplex * data, long numdata, double rin, double zin, double win, double *rout, double *zout, double *wout, rderivs * derivs) /* Return the Fourier frequency, f-dot, and fdotdot that */ @@ -13,16 +10,18 @@ double max_rzw_arr(fcomplex * data, int numdata, double rin, double zin, { float locpow = get_localpower3d(data, numdata, rin, zin, win); float maxpow = 0, pow, powargr, powargi; - int kern_half_width, extra = 10, startbin = (int)(rin) - 1; - int numz, numw, nextbin, fftlen, numbetween; + int extra = 10, startbin = (int)(rin) - 1; + int numr, numz, numw, nextbin, fftlen, numbetween; // The factor beyond ACCEL_DR, ACCEL_DZ, ACCEL_DW will interpolate int interpfac = 4, wind = 0, zind = 0, rind = 0; - double ifrac = 1.0/interpfac; fcomplex ***vol, amp; + // Use a more conservative length kernel + int kern_half_width = w_resp_halfwidth(fabs(zin)+2, fabs(win)+20, LOWACC); - kern_half_width = w_resp_halfwidth(zin, win, LOWACC); + numr = 3 * interpfac * ACCEL_RDR; numz = numw = 2 * interpfac + 1; - numbetween = interpfac*ACCEL_RDR; + + numbetween = interpfac * ACCEL_RDR; fftlen = next2_to_n(numbetween * (2 * kern_half_width + extra)); // printf("fftlen = %d\n", fftlen); vol = corr_rzw_vol(data, numdata, numbetween, startbin, @@ -33,7 +32,7 @@ double max_rzw_arr(fcomplex * data, int numdata, double rin, double zin, int ii, jj, kk; for (ii = 0; ii < numw; ii++) { for (jj = 0; jj < numz; jj++) { - for (kk = 0; kk < 3*interpfac*ACCEL_RDR; kk++) { + for (kk = 0; kk < numr; kk++) { amp = vol[ii][jj][kk]; pow = POWER(amp.r, amp.i); if (pow > maxpow) { @@ -46,9 +45,9 @@ double max_rzw_arr(fcomplex * data, int numdata, double rin, double zin, } } } - *rout = startbin + rind * ACCEL_DR * ifrac; - *zout = zin-ACCEL_DZ + zind * ACCEL_DZ * ifrac; - *wout = win-ACCEL_DW + wind * ACCEL_DW * ifrac; + *rout = startbin + rind * ACCEL_DR / (double) interpfac; + *zout = zin-ACCEL_DZ + zind * ACCEL_DZ / (double) interpfac; + *wout = win-ACCEL_DW + wind * ACCEL_DW / (double) interpfac; vect_free(vol[0][0]); vect_free(vol[0]); vect_free(vol); @@ -56,19 +55,22 @@ double max_rzw_arr(fcomplex * data, int numdata, double rin, double zin, return maxpow; } + double max_rzw_file(FILE * fftfile, double rin, double zin, double win, double *rout, double *zout, double *wout, rderivs * derivs) /* Return the Fourier frequency, f-dot, and fdotdot that */ /* maximizes the power of the candidate in 'fftfile'. */ { double maxpow, rin_int, rin_frac; - int kern_half_width, filedatalen, startbin, extra = 10; + int filedatalen, extra = 10; + long startbin; fcomplex *filedata; + // Use a more conservative length kernel + int kern_half_width = w_resp_halfwidth(fabs(zin)+2, fabs(win)+20, LOWACC); rin_frac = modf(rin, &rin_int); - kern_half_width = w_resp_halfwidth(zin, win, HIGHACC); filedatalen = 2 * kern_half_width + extra; - startbin = (int) rin_int - filedatalen / 2; + startbin = (long) rin_int - filedatalen / 2; filedata = read_fcomplex_file(fftfile, startbin, filedatalen); maxpow = max_rzw_arr(filedata, filedatalen, rin_frac + filedatalen / 2, @@ -77,3 +79,113 @@ double max_rzw_file(FILE * fftfile, double rin, double zin, double win, vect_free(filedata); return maxpow; } + + +void max_rzw_arr_harmonics(fcomplex data[], long numdata, + int num_harmonics, + double rin, double zin, double win, + double *rout, double *zout, double *wout, + rderivs derivs[], double powers[]) +/* Return the Fourier frequency, f-dot, and f-dotdot that */ +/* maximizes the *summed* power of the multi-harmonic candidate */ +{ + float maxpow = 0, pow, powargr, powargi, locpow; + long lobin, hhlobin = 0; + int hind, ii, jj, kk, extra = 10, nextbin, fftlen; + int interpfac = 4; // Factor beyond ACCEL_D[RZW] we will interpolate + int numwid = 2; // The number of full peak widths we will search + int numbetween = interpfac * ACCEL_RDR; + int numr = 3 * numbetween; // Search 3 full Fourier bins around high harm + int numz = 2 * numwid * interpfac + 1; + int numw = numz; + int rind = 0, zind = 0, wind = 0; + double dr = 1.0 / (double) numbetween; + + // The summed power spectrum, initialized to zeros + float ***powsum = gen_f3Darr(numw, numz, numr); + for (ii = 0; ii < numw * numz * numr; ii++) + powsum[0][0][ii] = 0.0; + + for (hind = 0; hind < num_harmonics; hind++) { + int n = num_harmonics - hind; // harmonic number, starting from highest + double rh = rin * n, zh = zin * n, wh = win * n; + // Use a more conservative length kernel + int kern_half_width = w_resp_halfwidth(fabs(zh)+2, fabs(wh)+20, LOWACC); + fcomplex ***vol, amp; + double rh_int, rh_frac, hfrac = n / (double) num_harmonics; + + locpow = get_localpower3d(data, numdata, rh, zh, wh); + rh_frac = modf(rh, &rh_int); + // Will do 1+ bins below and 1+ bins above rin + lobin = (long) rh_int - 1; + if (hind==0) hhlobin = lobin; + fftlen = next2_to_n(numbetween * (2 * kern_half_width + extra)); + // Create the RZW volume for the harmonic. + // Note that we are computing the z and w values in exact harmonic + // ratios. But the r values are on a power-of-two grid. + vol = corr_rzw_vol(data, numdata, numbetween, lobin, + zh-numwid*hfrac*ACCEL_DZ, + zh+numwid*hfrac*ACCEL_DZ, numz, + wh-numwid*hfrac*ACCEL_DW, + wh+numwid*hfrac*ACCEL_DW, numw, + fftlen, LOWACC, &nextbin); + // Calculate and sum the powers + for (ii = 0; ii < numw; ii++) { + for (jj = 0; jj < numz; jj++) { + for (kk = 0; kk < numr; kk++) { + rind = (long) round(((hhlobin + kk * dr) * hfrac + - lobin) * numbetween); + amp = vol[ii][jj][rind]; + powsum[ii][jj][kk] += POWER(amp.r, amp.i) / locpow; + } + } + } + vect_free(vol[0][0]); + vect_free(vol[0]); + vect_free(vol); + } + + // Now search the power sums for the highest value + rind = zind = wind = 0; + for (ii = 0; ii < numw; ii++) { + for (jj = 0; jj < numz; jj++) { + for (kk = 0; kk < numr; kk++) { + pow = powsum[ii][jj][kk]; + if (pow > maxpow) { + maxpow = pow; + wind = ii; + zind = jj; + rind = kk; + } + } + } + } + + // Calculate the proper r, z, and w peak values + { + double nh = num_harmonics; + *rout = (hhlobin + rind * dr) / nh; + *zout = ((zin * nh) - numwid * ACCEL_DZ + + zind * ACCEL_DZ / (double) interpfac) / nh; + *wout = ((win * nh) - numwid * ACCEL_DW + + wind * ACCEL_DW / (double) interpfac) / nh; + } + + // Now calculate the derivatives at the peak + for (ii = 0; ii < num_harmonics; ii++) { + int hh = ii + 1; + double rr = *rout * hh, zz = *zout * hh, ww = *wout * hh; + locpow = get_localpower3d(data, numdata, rr, zz, ww); + get_derivs3d(data, numdata, rr, zz, ww, locpow, &(derivs[ii])); + powers[ii] = derivs[ii].pow; + } + /* + printf("numr = %d numz = %d numw = %d\n", numr, numz, numw); + printf("rind = %d zind = %d wind = %d\n", rind, zind, wind); + printf("rin = %f zin = %f win = %f\n", rin , zin , win); + printf("rout = %f zout = %f wout = %f\n", *rout, *zout, *wout); + */ + vect_free(powsum[0][0]); + vect_free(powsum[0]); + vect_free(powsum); +} diff --git a/src/prepfold.c b/src/prepfold.c index 5aef2990a..0c3fcce6e 100644 --- a/src/prepfold.c +++ b/src/prepfold.c @@ -162,7 +162,6 @@ int main(int argc, char *argv[]) if (!cmd->outfileP) { char *tmprootnm, *suffix; - printf("XXX %s\n", cmd->argv[0]); split_root_suffix(cmd->argv[0], &tmprootnm, &suffix); if ((cmd->startT != 0.0) || (cmd->endT != 1.0)) { rootnm = (char *) calloc(strlen(tmprootnm) + 11, sizeof(char)); @@ -346,38 +345,20 @@ int main(int argc, char *argv[]) /* The following allows using inf files from searches of a subset */ /* of events from an event file. */ - if (cmd->rzwcandP || cmd->accelcandP) { + if (cmd->accelcandP) { infodata rzwidata; char *cptr = NULL; - if (cmd->rzwcandP) { - if (!cmd->rzwfileP) { - printf("\nYou must enter a name for the rzw candidate "); - printf("file (-rzwfile filename)\n"); - printf("Exiting.\n\n"); - exit(1); - } else if (NULL != (cptr = strstr(cmd->rzwfile, "_rzw"))) { - ii = (long) (cptr - cmd->rzwfile); - } else if (NULL != (cptr = strstr(cmd->rzwfile, "_ACCEL"))) { - ii = (long) (cptr - cmd->rzwfile); - } - cptr = (char *) calloc(ii + 1, sizeof(char)); - strncpy(cptr, cmd->rzwfile, ii); - } - if (cmd->accelcandP) { - if (!cmd->accelfileP) { - printf("\nYou must enter a name for the ACCEL candidate "); - printf("file (-accelfile filename)\n"); - printf("Exiting.\n\n"); - exit(1); - } else if (NULL != (cptr = strstr(cmd->accelfile, "_rzw"))) { - ii = (long) (cptr - cmd->accelfile); - } else if (NULL != (cptr = strstr(cmd->accelfile, "_ACCEL"))) { - ii = (long) (cptr - cmd->accelfile); - } - cptr = (char *) calloc(ii + 1, sizeof(char)); - strncpy(cptr, cmd->accelfile, ii); + if (!cmd->accelfileP) { + printf("\nYou must enter a name for the ACCEL candidate "); + printf("file (-accelfile filename)\n"); + printf("Exiting.\n\n"); + exit(1); + } else if (NULL != (cptr = strstr(cmd->accelfile, "_ACCEL"))) { + ii = (long) (cptr - cmd->accelfile); } + cptr = (char *) calloc(ii + 1, sizeof(char)); + strncpy(cptr, cmd->accelfile, ii); readinf(&rzwidata, cptr); free(cptr); idata.mjd_i = rzwidata.mjd_i; @@ -403,7 +384,7 @@ int main(int argc, char *argv[]) events = read_events(s.files[0], cmd->doubleP, eventtype, &numevents, idata.mjd_i + idata.mjd_f, idata.N * idata.dt, cmd->startT, cmd->endT, cmd->offset); - if (cmd->rzwcandP || cmd->accelcandP) + if (cmd->accelcandP) T = idata.N * idata.dt; else { /* The 1e-8 prevents floating point rounding issues from @@ -437,14 +418,14 @@ int main(int argc, char *argv[]) get_psr_from_parfile(cmd->parname, 51000.0, &psr); search.candnm = (char *) calloc(strlen(psr.jname) + 5, sizeof(char)); sprintf(search.candnm, "PSR_%s", psr.jname); - } else if (cmd->rzwcandP) { - slen = 20; - search.candnm = (char *) calloc(slen, sizeof(char)); - sprintf(search.candnm, "RZW_Cand_%d", cmd->rzwcand); } else if (cmd->accelcandP) { + char *cptr = NULL; slen = 22; search.candnm = (char *) calloc(slen, sizeof(char)); - sprintf(search.candnm, "ACCEL_Cand_%d", cmd->accelcand); + if (NULL != (cptr = strstr(cmd->accelfile, "_JERK"))) + sprintf(search.candnm, "JERK_Cand_%d", cmd->accelcand); + else + sprintf(search.candnm, "ACCEL_Cand_%d", cmd->accelcand); } else { slen = 20; search.candnm = (char *) calloc(slen, sizeof(char)); @@ -740,49 +721,35 @@ int main(int argc, char *argv[]) search.orb.w = (cmd->w + dtmp * cmd->wdot / SECPERJULYR); binary = 1; - } else if (cmd->rzwcandP || cmd->accelcandP) { + } else if (cmd->accelcandP) { fourierprops rzwcand; infodata rzwidata; char *cptr = NULL; + double T, r0, z0; - if (cmd->rzwcandP) { - if (!cmd->rzwfileP) { - printf("\nYou must enter a name for the rzw candidate "); - printf("file (-rzwfile filename)\n"); - printf("Exiting.\n\n"); - exit(1); - } else if (NULL != (cptr = strstr(cmd->rzwfile, "_rzw"))) { - ii = (long) (cptr - cmd->rzwfile); - } else if (NULL != (cptr = strstr(cmd->rzwfile, "_ACCEL"))) { - ii = (long) (cptr - cmd->rzwfile); - } - cptr = (char *) calloc(ii + 1, sizeof(char)); - strncpy(cptr, cmd->rzwfile, ii); - } - if (cmd->accelcandP) { - if (!cmd->accelfileP) { - printf("\nYou must enter a name for the ACCEL candidate "); - printf("file (-accelfile filename)\n"); - printf("Exiting.\n\n"); - exit(1); - } else if (NULL != (cptr = strstr(cmd->accelfile, "_rzw"))) { - ii = (long) (cptr - cmd->accelfile); - } else if (NULL != (cptr = strstr(cmd->accelfile, "_ACCEL"))) { - ii = (long) (cptr - cmd->accelfile); - } - cptr = (char *) calloc(ii + 1, sizeof(char)); - strncpy(cptr, cmd->accelfile, ii); + if (!cmd->accelfileP) { + printf("\nYou must enter a name for the ACCEL candidate "); + printf("file (-accelfile filename)\n"); + printf("Exiting.\n\n"); + exit(1); + } else if (NULL != (cptr = strstr(cmd->accelfile, "_ACCEL"))) { + ii = (long) (cptr - cmd->accelfile); } + cptr = (char *) calloc(ii + 1, sizeof(char)); + strncpy(cptr, cmd->accelfile, ii); printf("\nAttempting to read '%s.inf'. ", cptr); readinf(&rzwidata, cptr); free(cptr); printf("Successful.\n"); - if (cmd->rzwfileP) - get_rzw_cand(cmd->rzwfile, cmd->rzwcand, &rzwcand); - if (cmd->accelfileP) - get_rzw_cand(cmd->accelfile, cmd->accelcand, &rzwcand); - f = (rzwcand.r - 0.5 * rzwcand.z) / (rzwidata.dt * rzwidata.N); - fd = rzwcand.z / ((rzwidata.dt * rzwidata.N) * (rzwidata.dt * rzwidata.N)); + get_rzw_cand(cmd->accelfile, cmd->accelcand, &rzwcand); + T = rzwidata.dt * rzwidata.N; + // fourier props file reports average r and average z. + // We need the starting values. + z0 = rzwcand.z - 0.5 * rzwcand.w; + r0 = rzwcand.r - 0.5 * z0 - rzwcand.w / 6.0; + f = r0 / T; + fd = z0 / (T * T); + fdd = rzwcand.w / (T * T * T); /* Now correct for the fact that we may not be starting */ /* to fold at the same start time as the rzw search. */ @@ -802,7 +769,7 @@ int main(int argc, char *argv[]) /* Determine the pulsar parameters to fold if we are not getting */ /* the data from a .cand file, the pulsar database, or a makefile. */ - if (!cmd->rzwcandP && !cmd->accelcandP && !cmd->psrnameP && !cmd->parnameP) { + if (!cmd->accelcandP && !cmd->psrnameP && !cmd->parnameP) { double p = 0.0, pd = 0.0, pdd = 0.0; if (cmd->pP) { diff --git a/src/read_fft.c b/src/read_fft.c index 0c837f146..a2cf16675 100644 --- a/src/read_fft.c +++ b/src/read_fft.c @@ -1,7 +1,7 @@ #include "chkio.h" #include "ransomfft.h" -fcomplex *read_fcomplex_file(FILE * file, int firstpt, int numpts) +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. */ @@ -12,11 +12,11 @@ fcomplex *read_fcomplex_file(FILE * file, int firstpt, int numpts) /* If the number of bins to read takes us past the end of */ /* file, the returned vector will be zero padded. */ { - int ii, startpt = 0, lesspad = 0; + long ii, startpt = 0, lesspad = 0; fcomplex *result, *fcptr = NULL, zeros = { 0.0, 0.0 }; if (numpts < 0) { - printf("\n\n numpts = %d (out-of-bounds) in read_fcomplex_file().", numpts); + printf("\n\n numpts = %ld (out-of-bounds) in read_fcomplex_file().", numpts); printf(" Exiting.\n\n"); exit(1); } @@ -47,7 +47,7 @@ fcomplex *read_fcomplex_file(FILE * file, int firstpt, int numpts) } -float *read_float_file(FILE * file, int firstpt, int numpts) +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. */ @@ -58,11 +58,11 @@ float *read_float_file(FILE * file, int firstpt, int numpts) /* If the number of bins to read takes us past the end of */ /* file, the returned vector will be zero padded. */ { - int ii, startpt = 0, lesspad = 0; + long ii, startpt = 0, lesspad = 0; float *result, *fptr = NULL; if (numpts < 0) { - printf("\n\n numpts = %d (out-of-bounds) in read_float_file().", numpts); + printf("\n\n numpts = %ld (out-of-bounds) in read_float_file().", numpts); printf(" Exiting.\n\n"); exit(1); } diff --git a/src/zapping.c b/src/zapping.c index af9a6f4ae..3aeecf48c 100644 --- a/src/zapping.c +++ b/src/zapping.c @@ -17,7 +17,7 @@ float calc_median_powers(fcomplex * amplitudes, int numamps) } fcomplex *get_rawbins(FILE * fftfile, double bin, - int numtoget, float *med, int *lobin) + int numtoget, float *med, long *lobin) { fcomplex *result; *lobin = (int) bin - numtoget / 2; @@ -28,13 +28,13 @@ fcomplex *get_rawbins(FILE * fftfile, double bin, void zapbirds(double lobin, double hibin, FILE * fftfile, fcomplex * fft) { - int ii, ilobin, ihibin, binstozap, lodatabin; + long ii, ilobin, ihibin, binstozap, lodatabin; float median_lo, median_hi, avgamp; /* double phase, radargr, radargi, radtmp; */ fcomplex *data; - ilobin = (int) floor(lobin); - ihibin = (int) ceil(hibin); + ilobin = (long) floor(lobin); + ihibin = (long) ceil(hibin); binstozap = ihibin - ilobin; if (lobin - 1.5 * MEDIANBINS > 1) { if (fftfile) { /* If we are reading a file */ diff --git a/tests/test_max_rzw_harmonics.py b/tests/test_max_rzw_harmonics.py index 1be68b274..2b30270d9 100644 --- a/tests/test_max_rzw_harmonics.py +++ b/tests/test_max_rzw_harmonics.py @@ -2,12 +2,11 @@ import presto from numpy.random import standard_normal as norm from numpy.random import uniform -import matplotlib.pyplot as plt import time N = 2**17 -noiseamp = 40.0 -numharm = 4 +noiseamp = 1.0 +numharm = 1 numtrials = 100 us = np.arange(N, dtype=np.float64) / N # normalized time coordinate @@ -22,7 +21,7 @@ r = N/(4*numharm) + uniform(0.0, 1.0, 1)[0] # average freq over "observation" z = uniform(-100, 100, 1)[0] # average fourier f-dot w = uniform(-600, 600, 1)[0] # fourier freq double deriv - # w = 0.0 # fourier freq double deriv + w = 0.0 # fourier freq double deriv data = np.zeros_like(us) for ii in range(numharm): rh = r * (ii + 1) @@ -35,8 +34,8 @@ data += noiseamp * norm(N) ft = presto.rfft(data) - offset = uniform(-1.0, 1.0, 3) * np.array([0.5, 2.0, 20.0]) / numharm - + offset = uniform(-1.0, 1.0, 3) * np.array([0.5, 2.0, 20.0]) / (0.5 * numharm) + a = time.clock() if (numharm > 1): [maxpow, rmax, zmax, rds] = presto.maximize_rz_harmonics(ft, r+offset[0], @@ -73,7 +72,3 @@ print "rzwerrs:" print " avg: %6.3f %6.3f %6.3f %6.3f" % tuple(rzwerrs.mean(axis=0)) print " std: %6.3f %6.3f %6.3f %6.3f" % tuple(rzwerrs.std(axis=0)) - -plt.hist(rzerrs[:, 0], bins=50, color='blue', alpha=0.5) -plt.hist(rzwerrs[:,0], bins=50, color='red', alpha=0.5) -plt.show()