diff --git a/src/data/timing/colacal.sl b/src/data/timing/colacal.sl index 8ab0dbdea40d63aeb4daef6632c6c2eefdb151cd..9a4f6bd9cfae17d6c13084a2f45cb8856bec3700 100644 --- a/src/data/timing/colacal.sl +++ b/src/data/timing/colacal.sl @@ -4,7 +4,7 @@ define colacal(freq, cpd, noicpd, lopsd, hipsd, noilopsd, noihipsd, siglopsd, si %!%+ %\function{colacal} %\synopsis{Timing Tools: Coherence and Lag Calculation} -%\usage{ ([6 outputs]) = colacal ([10 required inputs]);} +%\usage{ Struct_Type cola = colacal ([10 required inputs]);} %\description % Input:\n % freq - Fourier Frequency Array\n @@ -17,12 +17,13 @@ define colacal(freq, cpd, noicpd, lopsd, hipsd, noilopsd, noihipsd, siglopsd, si % siglopsd - Noise-corrected PSD (low channel)\n % sighipsd - Noise-corrected PSD (high channel)\n % alln - number of averaged segments in each frequency bin\n -% Output:\n +% Output: Struct containing\n % rawcof - non-noise-corrected coherence function [Vaughan & Nowak, 1997, ApJ, 474, L43 (Eqn. 2)]\n % cof - noise-corrected coherence function [Vaughan & Nowak, 1997, ApJ, 474, L43 (Eqn. 8, Part 1)]\n % errcof - one-sigma uncertainty of cof [Vaughan & Nowak, 1997, ApJ, 474, L43 (Eqn. 8, Part 2)]\n % lag - time lag [Nowak et al., 1999, ApJ, 510, 874 (Sect. 4)]\n % errlag - one-sigma uncertainty of lag [Nowak et al., 1999, ApJ, 510, 874 (Eqn. 16)]\n +% errlag_Uttley14 - one-sigma uncertainty of lag [Uttley et al., 2014, A&A Rv. 22, 72 (Eq. 11-12)] % sigcpd - noise-corrected cross-power-density [sigcpd = cpd - noicpd]\n %!%- { @@ -35,5 +36,10 @@ define colacal(freq, cpd, noicpd, lopsd, hipsd, noilopsd, noihipsd, siglopsd, si ); variable lag = atan2(Imag(cpd), Real(cpd)) / (2.*PI*freq); variable errlag = sqrt( (1.-rawcof)/(2.*rawcof)/alln ) / (2.*PI*freq); - return rawcof, cof, errcof, lag, errlag, sigcpd; + variable sqr_bias = ((lopsd-noilopsd)*noihipsd+(hipsd-noihipsd)*noilopsd+noilopsd*noihipsd)/alln; % Uttley+14, Footnote 4 + variable rawcof_w_bias = (sqr(abs(cpd))-sqr_bias) / (lopsd * hipsd); % Uttley+14 Eq. 11 + variable errlag_uttley14 = sqrt( (1.-rawcof_w_bias)/(2.*rawcof)/alln ) / (2.*PI*freq); % Uttley+14 Eq. 12 + return struct{rawcof=rawcof, cof=cof, errcof=errcof, + lag=lag, errlag=errlag, errlag_uttley14=errlag_uttley14, + sigcpd=sigcpd}; } diff --git a/src/data/timing/cross_power_density.sl b/src/data/timing/cross_power_density.sl new file mode 100644 index 0000000000000000000000000000000000000000..a4f9734aa70fd198dde64de9ad544a06315e0e6c --- /dev/null +++ b/src/data/timing/cross_power_density.sl @@ -0,0 +1,44 @@ +private define cross_power_density(rate1, rate2, numseg, dimseg){ +%!%+ +%\function{cross_power_density} +%\synopsis{Timing Tools: Cross Power Density (aka. Cross-Spectrum)} +%\usage{(cpd, errcpd) = cross_power_density(rate1, rate2, numseg, dimseg)} +%\description +% Calculates the CPD, averaged over all segments, as +% Conj(DFT(rate1))*DFT(rate2) +% using the dorDFT function. The return values are of Complex_Type. +% The error is the standard error on the mean of the averaged CPD. +%\seealso{dorDFT} +%!%- + variable startbin = 0; + variable endbin = dimseg-1; + + variable cpd = Complex_Type[dimseg/2]; % = 0 + 0i + variable cpds = Complex_Type[numseg, dimseg/2]; % Record of CPDs of all segments for error calculation + + variable seg=0; + loop(numseg) + { + variable cpd_of_segment = Conj(dorDFT(rate1[[startbin:endbin]])) * dorDFT(rate2[[startbin:endbin]]); + + cpd += cpd_of_segment; + startbin += dimseg; + endbin += dimseg; + + % Divide by dimseg because dorDFT has normalization sqrt(N) and is + % applied twice: once for DFT(rate1), once for DFT(rate2). + cpds[seg,*] = cpd_of_segment / dimseg; + seg++; + } + + cpd /= (numseg * dimseg); % Average over all segments + + % ====== Error calculation ======================================== + variable errcpd = Complex_Type[dimseg/2]; + variable freq; + _for freq (0, dimseg/2-1, 1){ + errcpd[freq] = moment(Real(cpds[*,freq])).sdom + moment(Imag(cpds[*,freq])).sdom*1i; + } + + return cpd, errcpd; +} diff --git a/src/data/timing/dorDFT.sl b/src/data/timing/dorDFT.sl index 1ad36291db85e086016a842633cfde12c59bc670..793cabf193ab851634714fbf5b795dd2e2420eb9 100644 --- a/src/data/timing/dorDFT.sl +++ b/src/data/timing/dorDFT.sl @@ -2,11 +2,15 @@ private define dorDFT(rate) %!%+ %\function{dorDFT} %\synopsis{Timing Tools: Discrete Fourier Transform} -%\usage{(dft) = dorDFT(rate);} +%\usage{dft = dorDFT(rate);} %\description -% Performs a Fast Fourier Transform on a given rate array with -% the same properties as IDL FFT (Normalization = 1/sqrt(2*PI)) AND -% returns only the meaningful bins for PSD calculations. +% Performs a Fast Fourier Transform on a given rate array with the +% same properties as IDL FFT: +% * Normalization = sqrt(length(rate)) +% * renormalize rate around 0 by subtracting mean such that +% variability is calculated wrt. mean rate +% * return only meaningful bins for PSD calculations, i.e. first +% bin up to Nyquist frequency %!%- { variable l = length(rate); diff --git a/src/data/timing/foucalc.sl b/src/data/timing/foucalc.sl index ff465229990e953e99e1a0c65bdbbc669810e40c..ec2d8ceb5c9ef6b5456b047fb233492e4a514f3e 100644 --- a/src/data/timing/foucalc.sl +++ b/src/data/timing/foucalc.sl @@ -26,8 +26,9 @@ define foucalc() % \code{lc = struct { time=[t_1, t_2, ...], rate1=[r1_1, r1_2, ...], rate2=... };}. % [However, one can also use an array of structures (with an enormous overhead), % \code{lc = [ struct { time=t_1, rate1=r1_1, rate2=...}, struct { time=t_2, rate1=r1_2, rate2=... }, ... ];}.] +% Also specifying "rate" instead of "rate1" is possible. % -% \code{dimseg} is the segmentsize used for the FFTs, whihch should therefore be a power of 2. +% \code{dimseg} is the segment size used for the FFTs, which should therefore be a power of 2. % % The returned structure contains the following fields:\n % - \code{freq}: @@ -51,11 +52,14 @@ define foucalc() % - Cross power spectra, coherence and time lag functions for each pair of energy bands:\n % * \code{realcpd}, \code{imagcpd}: % real and imaginary part of the cross power density\n +% * \code{errrealcpd}, \code{errimagcpd}: +% standard error on the mean of the averaged cross power density\n % * \code{noicpd = ( sigpsd_lo * noipsd_hi + sigpsd_hi * noipsd_lo + noipsd_lo * noipsd_hi ) / numseg}\n % * \code{rawcof}, \code{cof}, \code{errcof}: % non-noise-corrected (raw) and noise-corrected coherence function and its one-sigma uncertainty\n % * \code{lag}, \code{errlag}: % time lag and its one-sigma uncertainty\n +%\seealso{makepsd, cross_power_density, colacal} %!%- { variable lc, dimseg; @@ -84,7 +88,7 @@ define foucalc() fieldnames = [fieldnames, ["rawpsd", "errpsd", "normpsd", "errnormpsd", "noipsd", "sigpsd", "noinormpsd", "signormpsd", "effnoinormpsd"]+string(i)]; _for i (1, numchan, 1) _for j (i+1, numchan, 1) - fieldnames = [fieldnames, ["realcpd", "imagcpd", "noicpd", "rawcof", "cof", "errcof", "lag", "errlag"]+string(i)+string(j)]; % assuming i < j < 100 ... + fieldnames = [fieldnames, ["realcpd", "imagcpd", "errrealcpd", "errimagcpd", "noicpd", "rawcof", "cof", "errcof", "lag", "errlag"]+string(i)+string(j)]; % assuming i < j < 100 ... variable fouquant = @Struct_Type(fieldnames); % create structure % Calculating Fourier Frequency Array @@ -113,9 +117,10 @@ define foucalc() _for i (1, numchan, 1) { variable istr = string(i); - if(verbose) ()=printf("Calculating power spectrum for rate%d. ", i); + if (verbose) ()=printf("Calculating power spectrum for rate%d. ", i); - variable rate = array_struct_field(lc, "rate"+istr); + % If "rate1" is not available in struct, use "rate" + variable rate = struct_field_exists(lc, "rate"+istr) ? array_struct_field(lc, "rate"+istr) : array_struct_field(lc, "rate"); avgrate[i-1] = mean(rate); % Calculating PSDs @@ -189,23 +194,17 @@ define foucalc() variable lostr = string(lo); variable histr = string(hi); - variable startbin = 0; - variable endbin = dimseg-1; - if(verbose) vmessage("Calculating cross power spectrum (=> coherence and timelags) for rate%d and rate%d.", lo, hi); - variable cpd = Complex_Type[dimseg/2]; % = 0 + 0i - loop(numseg) - { - cpd += Conj( dorDFT(array_struct_field(lc, "rate"+lostr)[[startbin:endbin]]) ) - * dorDFT(array_struct_field(lc, "rate"+histr)[[startbin:endbin]]); - startbin += dimseg; - endbin += dimseg; - } - cpd /= (numseg * dimseg); - % write Real(cpd) and Imag(cpd) to output structures - variable Re_cpd = Real(cpd); - variable Im_cpd = Imag(cpd); - set_struct_field(fouquant, "realcpd"+lostr+histr, Re_cpd); - set_struct_field(fouquant, "imagcpd"+lostr+histr, Im_cpd); + if (verbose) vmessage("Calculating cross power spectrum (=> coherence and timelags) for rate%d and rate%d.", lo, hi); + + % If "rate1" is not available in struct, use "rate" + variable rate1 = struct_field_exists(lc, "rate"+lostr) ? array_struct_field(lc, "rate"+lostr) : array_struct_field(lc, "rate"); + variable rate2 = array_struct_field(lc, "rate"+histr); + variable cpd, errcpd; + (cpd, errcpd) = cross_power_density(rate1, rate2, numseg, dimseg); + set_struct_field(fouquant, "realcpd"+lostr+histr, Real(cpd)); + set_struct_field(fouquant, "imagcpd"+lostr+histr, Imag(cpd)); + set_struct_field(fouquant, "errrealcpd"+lostr+histr, Real(errcpd)); + set_struct_field(fouquant, "errimagcpd"+lostr+histr, Imag(errcpd)); % use SIGPSD to calculate the cross power density noise variable sigpsd_lo = get_struct_field(fouquant, "sigpsd"+lostr); @@ -218,17 +217,15 @@ define foucalc() variable psd_lo = get_struct_field(fouquant,"rawpsd"+lostr); variable psd_hi = get_struct_field(fouquant,"rawpsd"+histr); - variable rawcof, cof, errcof, lag, errlag; - (rawcof, cof, errcof, lag, errlag,) - = colacal(fouquant.freq, cpd, noicpd, psd_lo, psd_hi, noipsd_lo, noipsd_hi, sigpsd_lo, sigpsd_hi, numseg); % fouquant.numavgall - set_struct_field(fouquant, "rawcof"+lostr+histr, rawcof); - set_struct_field(fouquant, "cof" +lostr+histr, cof); - set_struct_field(fouquant, "errcof"+lostr+histr, errcof); - set_struct_field(fouquant, "lag" +lostr+histr, lag); - set_struct_field(fouquant, "errlag"+lostr+histr, errlag); + variable cola = colacal(fouquant.freq, cpd, noicpd, psd_lo, psd_hi, noipsd_lo, noipsd_hi, sigpsd_lo, sigpsd_hi, numseg); % fouquant.numavgall + set_struct_field(fouquant, "rawcof"+lostr+histr, cola.rawcof); + set_struct_field(fouquant, "cof" +lostr+histr, cola.cof); + set_struct_field(fouquant, "errcof"+lostr+histr, cola.errcof); + set_struct_field(fouquant, "lag" +lostr+histr, cola.lag); + set_struct_field(fouquant, "errlag"+lostr+histr, cola.errlag); %%% set_struct_field(fouquant, "sigcpd"+lostr+histr, sigcpd); - rawcof=NULL; cof=NULL; errcof=NULL; lag=NULL; errlag=NULL; - cpd = NULL; + cola= NULL; + cpd = NULL; errcpd=NULL; } } diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl new file mode 100644 index 0000000000000000000000000000000000000000..3f2491c4c87f44a98ee5be0f52fa8ef9203f7bf6 --- /dev/null +++ b/src/data/timing/lag_energy_spectrum.sl @@ -0,0 +1,307 @@ +private define filterFouQuantForFrequencyRange(fouquan, f_min, f_max) { + variable filt = where(f_min < fouquan.freq <= f_max); + if (qualifier("verbose",0)>1){ + vmessage(" Minimal frequency bin %d at %g Hz", filt[0], fouquan.freq[filt[0]]); + vmessage(" Maximal frequency bin %d at %g Hz", filt[-1], fouquan.freq[filt[-1]]); + vmessage(" Number of frequency bins: %d", length(filt)); + } + return struct_filter(fouquan, filt ; copy); +} + +private define verifyLightcurveTime(lc1, lc2){ + %% Test whether the time interval of reference and energy-resolved + %% lightcurve is the same + if (length(lc1.time)!=length(lc2.time) or (lc1.time[0]!=lc2.time[0]) or (lc1.time[-1]!=lc2.time[-1])){ + exit("ERROR: The time arrays of the reference and small band lightcurve seem to be different!"); + } +} + +private define getLightcurve(lc, dt, dimseg){ + return segment_lc_for_psd(fits_read_table(lc), dt, dimseg ;; __qualifiers); +} + +private define calculateSummedReferenceLightcurve(lc_list, index){ + %% As reference lightcurve, take the summed lightcurves of all + %% energy bands but the one we do the cross-spectrum with (the + %% lightcurve at index, often called the "channel of interest"). + variable channel_of_interest = fits_read_table(lc_list[index]); + + variable jj; + variable summed_rate = Double_Type[length(channel_of_interest.rate)]; + _for jj (0, length(lc_list)-1, 1){ + if (jj==index){ + continue; + } else { + variable lc = fits_read_table(lc_list[jj]); + verifyLightcurveTime(channel_of_interest, lc); + summed_rate += lc.rate; + % Error is not needed for cross-spectrum calculation + } + } + return struct{time=channel_of_interest.time, rate=summed_rate}; +} + +private define getReferenceLightcurve(lc_ref, cof_lc, index, dt, dimseg){ + % Default: calculate the reference lightcurve as sum of all + % energy-resolved lightcurves but the current channel-of-interest. + % The Poisson error will not be correlated and therefore we do not + % need to subtract off the channel-of-interest. In this case lc_ref + % is a list of lightcurves. + % + % If we use a broad reference lightcurve, which contains the + % channel-of-interest, we need to subtract off the small band LC + % from the reference lightcurve. Uttley+14 call this a + % channel-of-interest (CI) correction. This is needed to account + % for the Poisson noise which is correlated for both LCs. See + % also Uttley+14 Sect 3.2. In this case the lc_ref parameter is an + % already loaded, segmented lightcurve. + variable rate; + if (qualifier_exists("lc_ref")==1){ + rate = lc_ref.rate - cof_lc.rate; + } else { + variable summed_lc_ref = calculateSummedReferenceLightcurve(lc_ref, index); + lc_ref = segment_lc_for_psd(summed_lc_ref, dt, dimseg ;; __qualifiers); + rate = lc_ref.rate; + } + verifyLightcurveTime(lc_ref, cof_lc); + return struct{time=lc_ref.time, rate=rate}; +} + +define lag_energy_spectrum(lc_list, dimseg) { +%!%+ +%\function{lag_energy_spectrum} +%\synopsis{Calculates the time lag of energy-resolved lightcurves} +%\usage{Struct_Type les = lag_energy_spectrum(String_Type[lcs], dimseg)} +%\description +% * lc_list: Array of file names of the energy-resolved +% lightcurves (must have same length and time binning) +% * dimseg: Segmentation length in bins, needed to segment the +% lightcurve and do the Fourier calculation +% +% returns: lag-energy-spectrum as a struct with fields lag, errlag +%\qualifiers{ +% lc_ref - Filename of a (broad) reference lightcurve. If not +% set, the reference lightcurve will consist of the +% summed up energy lightcurves except for the current +% channel-of-interest. +% dt - Time resolution of the lightcurve in seconds +% [default: TIMEDEL keyword] +% f_min - For frequency-resolved. Default: lowest sampled frequency +% [1/(dt*dimseg)]. Unit: Hz +% f_max - See f_min. Default: Nyquist frequency [1/(2*dt)]. Unit: Hz +% verbose - Increase verbosity [default=0] +% outfile - String_Type or Fits_File_Type. If set, write output +% either to open FITS file pointer or create new file +% (which then contains only the current frequency +% interval) +% elo, ehi - Double_Type: Energy grid, only used for writing outfile +% +% (Furthermore, the qualifiers of segment_lc_for_psd can be specified) +%} +%\notes +% * The default error computation is Nowak et al., 1999, ApJ, 510, +% 874 (Sect. 4), i.e., from coherence without bias factor. +% * All input lightcurves are accounted for their gaps and segmented +% by segment_lc_for_psd. +% * The cross-spectrum (CPD) and lags are calculated with foucalc. +% * The PSD normalization is fixed to Miyamoto. +% * The average time lag in the frequency interval is calculated by +% the normal mean on the imaginary and real part of the +% cross-spectrum: +% atan2(mean(Imag(CPD)), mean(Real(CPD)))/(PI*(f_min+f_max)) +% * Reference lightcurve: For each energy bin, the function takes the +% summed lightcurves of all but the current energy bands (always +% excluding the lightcurve for the cross-spectrum). +% * The function does *not* use any information about the energy +% grid. This has to be created by yourself (and must be the same as +% the lightcurves in lc_list)! The elo, ehi qualifiers are only +% for the records in the output file +% +% If outfile is a string, a new FITS file will be created, but you +% can also give a file pointer and the function will only create a +% new extension called BAND__Hz. This way you can +% store frequency-resolved lag-energy spectra in one output file by +% always passing the same file pointer (see example). If you want +% to add energy information to the output FITS file, you can pass +% two Double_Type arrays (same length as lc_list) via the elo, ehi +% qualifiers. The columns of the created FITS file are: +% +% ENERGY_LO - Energy channel (keV), taken from elo qualifier +% ENERGY_HI - Energy channel (keV), taken from ehi qualifier +% MEAN_PHASE - Phase (radians), computed as atan2(mean(Im(CPD)), mean(Re(CPD))) +% LAG - Time lag (s), computed as phase/(PI*(f_min+f_max)) +% LAG_ERR_UTTLEY+14 - Error computed from coherence with bias +% Uttley et al., 2014, A&A Rv. 22, 72 (Eq. 11-12) +% LAG_ERR_NOWAK+99 - Error computed from coherence without bias (default) +% Nowak et al., 1999, ApJ, 510, 874 (Sect. 4) +% LAG_ERR_PROPERR - Error computed from Gaussian error propagation +% LAG_ERR_INGRAM19 - (to be implemented) +% MEAN_REAL/IMAG_CPD - real/imaginary part of the mean of the +% segment-averaged cross-spectrum +% MEAN_REAL/IMAG_CPD_ERR - standard error on the mean of the +% real/imaginary segment-averaged cross-spectrum +% N_SEGMENTS - Number of segments averaged over (foucalc.numavgall) +% N_FREQUENCIES - Number of frequencies averaged over (as in +% "mean of the CPD") +% MIN/MAX_FREQUENCY - Minimal/Maximal frequency of the frequency-filtered +% Fourier products - this can be different from +% f_min due to the discretization of the +% *Discrete* Fourier Transform +% RMS - Root Mean Squared variability of the +% channel-of-interest lightcurve +% RMS_REF - RMS variability of the reference lightcurve +% +% The last three fields should be the same for all lightcurves. +% +% If you find any bugs or have questions, please contact ole.koenig@fau.de +%\example +% %% Example for frequency-resolved lag-energy spectrum +% (flo, fhi) = log_grid(1, 50, n_freqs); % (Hz) +% % Energy grid must be the same as extracted lightcurves +% (elo, ehi) = log_grid(0.5, 10, n_energies); % (keV) +% variable fp = fits_open_file("test.fits", "c"); +% _for ii (0, n_freqs-1, 1) { % going through frequency range +% variable les = lag_energy_spectrum(lc_names, dimseg +% ; f_min = flo[ii], f_max=fhi[ii], +% verbose=2, outfile=fp, elo=elo, ehi=ehi); +% hplot_with_err(elo, ehi, les.lag, les.errlag); +% } +% fits_close_file(fp); +%\seealso{segment_lc_for_psd, foucalc, colacal} +%!%- + variable dt = qualifier("dt", fits_read_key(lc_list[0], "TIMEDEL")); + variable f_min = qualifier("f_min", 1/(dt*dimseg)); + variable f_max = qualifier("f_max", 1/(2*dt)); % Nyquist frequency + variable verbose = qualifier("verbose", 0); + variable normtype = "Miyamoto"; + + if (dt==NULL) { + return "ERROR: Couldn't receive the time resolution from TIMEDEL, need dt qualifier!"; + } + + if (verbose>0) vmessage("Frequency range %g-%g Hz", f_min, f_max); + + variable N_energies = length(lc_list); + variable rmss = Double_Type[N_energies]; + variable rms_refs = Double_Type[N_energies]; + variable lags = Double_Type[N_energies]; + variable errlag_nowak99 = Double_Type[N_energies]; + variable errlag_uttley14 = Double_Type[N_energies]; + + variable lc_ref; + if (qualifier_exists("lc_ref")==1){ + lc_ref = getLightcurve(qualifier("lc_ref"), dt, dimseg ;; __qualifiers); + } else { + lc_ref = lc_list; + } + + if (qualifier_exists("outfile")==1){ + variable fp; + if (typeof(qualifier("outfile")) == String_Type){ + fp = fits_open_file(qualifier("outfile"), "c"); + } else { + fp = qualifier("outfile"); + } + variable min_freqs = Double_Type[N_energies]; + variable max_freqs = Double_Type[N_energies]; + variable mean_realcpds = Double_Type[N_energies]; + variable mean_realcpd_errs = Double_Type[N_energies]; + variable mean_imagcpds = Double_Type[N_energies]; + variable mean_imagcpd_errs = Double_Type[N_energies]; + variable mean_phases = Double_Type[N_energies]; + variable errlag_properr = Double_Type[N_energies]; + variable errlag_ingram19 = Double_Type[N_energies]; + variable n_segments = Integer_Type[N_energies]; + variable n_frequencies = Integer_Type[N_energies]; + } + + %% Go through the energy-resolved lightcurves (channel-of-interests) + variable ii; + _for ii (0,N_energies-1,1){ + + variable channel_of_interest = getLightcurve (lc_list[ii], dt, dimseg ;; __qualifiers); + variable reference = getReferenceLightcurve (lc_ref, channel_of_interest, ii, dt, dimseg ;; __qualifiers); + + variable rms; + variable res = foucalc( struct{ time = channel_of_interest.time, + rate1 = reference.rate, + rate2 = channel_of_interest.rate}, + dimseg ; normtype=normtype, RMS=&rms); + (rmss[ii], rms_refs[ii]) = rms[0], rms[1]; + + variable res_f = filterFouQuantForFrequencyRange(res, f_min, f_max ; verbose=verbose); + + %% Coherence and lag calculation on averaged cross-spectrum + variable cola = colacal(0.5*(f_min+f_max), % mid frequency + mean(res_f.realcpd12)+mean(res_f.imagcpd12)*1i, + mean(res_f.noicpd12), % not needed for lag+error + mean(res_f.rawpsd1), mean(res_f.rawpsd2), + mean(res_f.noipsd1), mean(res_f.noipsd2), + mean(res_f.sigpsd1), mean(res_f.sigpsd2), % not needed for lag+err + res_f.numavgall[0]*length(res_f.freq)); % K*M with M=number of frequencies averaged over and K=number of segments + (lags[ii], errlag_nowak99[ii], errlag_uttley14[ii]) = cola.lag, cola.errlag, cola.errlag_uttley14; + + if (qualifier_exists("outfile")==1){ + variable m_ReCPD = mean(res_f.realcpd12); + variable m_ReCPD_err = mean(res_f.errrealcpd12); + variable m_ImCPD = mean(res_f.imagcpd12); + variable m_ImCPD_err = mean(res_f.errimagcpd12); + + %% Calculation from Gaussian error propagation + variable m_phase = atan2(m_ImCPD, m_ReCPD); % radians + errlag_properr[ii]=1/(PI*(f_min+f_max))* + sqrt((m_ReCPD/(m_ImCPD^2+m_ReCPD^2))^2*m_ImCPD_err^2 + + (m_ImCPD/(m_ImCPD^2+m_ReCPD^2))^2*m_ReCPD_err^2 + ); + + %% Implement error calculation from Ingram 2019, Eq. 19 + + min_freqs[ii] = min(res_f.freq); + max_freqs[ii] = max(res_f.freq); + mean_realcpds[ii] = m_ReCPD; + mean_realcpd_errs[ii] = m_ReCPD_err; + mean_imagcpds[ii] = m_ImCPD; + mean_imagcpd_errs[ii] = m_ImCPD_err; + mean_phases[ii] = m_phase; + n_segments[ii] = res_f.numavgall[0]; + n_frequencies[ii] = length(res_f.freq); + } + } + + if (qualifier_exists("outfile")==1){ + variable extname = sprintf("BAND_%g_%gHz", f_min, f_max); + variable sdata = struct{ + "ENERGY_LO"=qualifier("elo", Double_Type[N_energies]), + "ENERGY_HI"=qualifier("ehi", Double_Type[N_energies]), + "LAG"=lags, + "LAG_ERR_UTTLEY+14"=errlag_uttley14, + "LAG_ERR_NOWAK+99"=errlag_nowak99, + "LAG_ERR_PROPERR"=errlag_properr, + "LAG_ERR_INGRAM19"=errlag_ingram19, + "MEAN_IMAG_CPD"=mean_imagcpds, + "MEAN_IMAG_CPD_ERR"=mean_imagcpd_errs, + "MEAN_REAL_CPD"=mean_realcpds, + "MEAN_REAL_CPD_ERR"=mean_realcpd_errs, + "MEAN_PHASE"=mean_phases, + "N_SEGMENTS"=n_segments, + "N_FREQUENCIES"=n_frequencies, + "MIN_FREQUENCY"=min_freqs, + "MAX_FREQUENCY"=max_freqs, + "RMS"=rmss, + "RMS_REF"=rms_refs + }; + variable keys = struct{ + F_MIN=f_min, F_MAX=f_max, DIMSEG=dimseg, TIMEDEL=dt, NORMTYPE=normtype + }; + variable hist = struct{ + history = [sprintf("Wrote lag-energy spectrum in frequency range %f-%fHz", f_min, f_max)], + comment = ["Time unit is seconds, energy unit is keV, frequency unit is Hertz", + "Analyzed lightcurves:", lc_list] + }; + fits_write_binary_table(fp, extname, sdata, keys, hist); + + if (typeof(qualifier("outfile")) == String_Type) fits_close_file(fp); + } + + return struct{ lag = lags, errlag = errlag_nowak99 }; +} diff --git a/src/data/timing/segment_lc_for_psd.sl b/src/data/timing/segment_lc_for_psd.sl index c1916c494d7c6ae35acc4bab2bbc4e77eb97e284..27c5b735c9607e77c92f9658b4f9da3f299e5b62 100644 --- a/src/data/timing/segment_lc_for_psd.sl +++ b/src/data/timing/segment_lc_for_psd.sl @@ -1,4 +1,4 @@ -define segment_lc_for_psd(lc,dt,dimseg) { %{{{ +define segment_lc_for_psd(lc, dt, dimseg) { %{{{ %!%+ %\function{segment_lc_for_psd} %\synopsis{Function to segment a lightcurve in preparation for a PSD, @@ -10,7 +10,7 @@ define segment_lc_for_psd(lc,dt,dimseg) { %{{{ % to dt [default: 1.1] % ratefield - Name of rate field in lightcurve [default: rate] % verbose - Increase verbosity to display which parts of the -% lightcurve and what proportion was rejected. +% lightcurve and what proportion was rejected [default=0]. %} %\description % Given a lightcurve with gaps, this function segments it such that @@ -42,17 +42,14 @@ define segment_lc_for_psd(lc,dt,dimseg) { %{{{ % dimseg=16; % (bins) % lc_split = segment_lc_for_psd(lc_gaps, dt, dimseg); % res=foucalc(struct{time=lc_split.time,rate1=lc_split.rate},dimseg); -%\seealso{split_lc_at_gaps, get_smallest_segment_in_lc, foucalc} +%\seealso{split_lc_at_gaps, foucalc} %!%- - % define relative factor, telling how much the gap is relative to dt; default = 1.1; - variable gapfactor = qualifier("gapfactor",1.1); + variable gapfactor = qualifier("gapfactor", 1.1); + variable ratefield = qualifier("ratefield", "rate"); + variable verbose = qualifier("verbose", 0); - %allow qualifier that defines the name of the "rate" field; if not set, default is "rate" - variable ratefield = qualifier("ratefield","rate"); - - % split lightcurve at gaps - variable split_lc = split_lc_at_gaps(lc,dt*gapfactor); + variable split_lc = split_lc_at_gaps(lc, dt*gapfactor); variable segmented_lc = struct{ time=Double_Type[0], @@ -62,7 +59,8 @@ define segment_lc_for_psd(lc,dt,dimseg) { %{{{ %% keep a record of the total time and the time thrown away variable total_time=0.; % (s) - variable rejected_time=0.; + variable rejected_time=0.; + variable n_segments=0; variable ii; _for ii (0,length(split_lc)-1,1) { @@ -76,7 +74,8 @@ define segment_lc_for_psd(lc,dt,dimseg) { %{{{ % only do this if there is at least one segment if (nn > 0) { % note that each segment goes from 0 to nn*dimseg-1 - % (the -1 is the important part!) + % (the -1 is the important part!) + n_segments++; segmented_lc.time=[segmented_lc.time,split_lc[ii].time[[0:nn*dimseg-1]]]; segmented_lc.rate=[segmented_lc.rate, get_struct_field(split_lc[ii],ratefield)[[0:nn*dimseg-1]]]; @@ -85,15 +84,17 @@ define segment_lc_for_psd(lc,dt,dimseg) { %{{{ } } else { rejected_time += exposure; - if (qualifier_exists("verbose")) { - vmessage(" ***Skipping lightcurve segment with index %i, length %.2fs",ii,exposure); + if (verbose>1){ + vmessage(" ***Skipping lightcurve segment with index %i, length %gs",ii,exposure); } } } - if (qualifier_exists("verbose")) { - vmessage("\nRejected time: %.2f s (%.3f percent of total data time=%.2fs)\n", + if (verbose>0) { + vmessage("\nRejected time: %g s (%g percent of total data time=%gs)", rejected_time, rejected_time/total_time*100, total_time); + vmessage("Number of segments with length>%gs (dimseg=%d, dt=%gs): %d", + dt*dimseg, dimseg, dt, n_segments); } return segmented_lc;