From 54e19ea0015f68f5e71712634f97030a4125472e Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Tue, 15 Dec 2020 09:52:08 +0100 Subject: [PATCH 1/9] Add new function lag_energy_spectrum * Still needs a bit of debugging and testing --- src/data/timing/lag_energy_spectrum.sl | 119 +++++++++++++++++++++++++ 1 file changed, 119 insertions(+) create mode 100644 src/data/timing/lag_energy_spectrum.sl diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl new file mode 100644 index 00000000..9379ca11 --- /dev/null +++ b/src/data/timing/lag_energy_spectrum.sl @@ -0,0 +1,119 @@ +require("isisscripts"); + +define filter_lc_for_timeseg(lc, t_lo, t_hi){ + variable filt = where(t_lo <= lc.time < t_hi); + return struct_filter(lc, filt ; copy); +} + +private define filter_fouquan_for_frequency_range(fouquan, f_min, f_max) { + variable filt = where(f_min <= fouquan.freq <= f_max); + return struct_filter(fouquan, filt ; copy); +} + +define lag_energy_spectrum(lc_list, lc_ref, dimseg) { +%!%+ +%\function{lag_energy_spectrum} +%\synopsis{Calculates the time lag-energy spectrum between +%energy-resolved lightcurves and a broad reference lightcurve.} +%\usage{} +%\qualifiers{ +% tstart - The lightcurves can be constrained to specific time +% intervals (e.g., when the lightcurve changes +% dramatically in flux or shows really large gaps). Unit: s +% tstop - See tstart. Unit: s +% dt - The lightcurves are parsed to segment_lc_for_psd and +% split at the gaps. In order to do so, one needs the time +% resolution at is segmented and accounted for the +% gaps. dt is the time resolution of the lightcurve in +% seconds. +% f_min - The lag-energy spectrum can also be computed only in +% a small frequency band (for a frequency-resolved +% lag-energy spectrum). Default: lowest sampled frequency +% [1/(dt*dimseg)]. Unit: Hz +% f_max - See f_min. Default: Nyquist frequency [1/(2*dt)]. Unit: Hz +% normtype - Normalization type of PSD [default: Miyamoto]. +% raw - return an additional list of the raw lag-frequency spectra. +%} +%\description +% * lc_list: List of names of the energy-resolved lightcurves +% * lc_ref: Name of the (broad) reference lightcurve +% * dimseg: Segmentation length in bins, needed to segment the +% lightcurve and do the Fourier calculation +% +% returns: the lag-energy-spectrum as a struct with the fields +% lags, errs +%\notes +% The function does *not* use and will not return any information +% about the energy grid. This has to be created by yourself (and +% must be the same as the lightcurves in lc_list)! +%\example +%\seealso{segment_lc_for_psd, foucalc} +%!%- + + variable dt = qualifier("dt", fits_read_key(lc_ref, "TIMEDEL")); + if (dt==NULL) { + return "ERROR: Couldn't receive the time resolution from TIMEDEL lightcurve, need dt qualifier!"; + } + + variable f_min = qualifier("f_min", 1/(dt*dimseg)); + variable f_max = qualifier("f_max", 1/(2*dt)); % Nyquist frequency + variable normtype = qualifier("normtype","Miyamoto"); + + variable lc_ref_segm = segment_lc_for_psd(fits_read_table(lc_ref), dt, dimseg); + + variable tstart = qualifier("tstart", min(lc_ref_segm.time)); + variable tstop = qualifier("tstop", max(lc_ref_segm.time)); + + variable ref = filter_lc_for_timeseg(lc_ref_segm, tstart, tstop); + + variable tlags = Float_Type[length(lc_list)]; + variable tlagerrs = Float_Type[length(lc_list)]; + + if (qualifier_exists("raw")) variable raw = {}; + + %% Go through the energy-resolved lightcurves + variable ii; + _for ii (0,length(lc_list)-1,1){ + + variable lc_segm = segment_lc_for_psd(fits_read_table(lc_list[ii]), dt, dimseg); + + variable lc = filter_lc_for_timeseg(lc_segm, tstart, tstop); + + %% Test whether the time interval of reference and energy-resolved + %% lightcurve is the same + if (length(ref.time)!=length(lc.time) + or (ref.time[0]!=lc.time[0]) + or (ref.time[-1]!=lc.time[-1])){ + return sprintf("ERROR: The time arrays of %s and %s seem to be different!", + lc_ref, lc_list[ii]); + } + + %% Calculate the cross-spectrum (CPD) and lags with foucalc Note + %% that we subtract off the small band LC from the reference LC. + %% Uttley+11 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+11 Sect 3.2. + %% The time lag is computed from the CPD as: + %% phase = atan2(r.imagcpd12,r.realcpd12); + %% tlag = phase/(2*PI*r.freq); + variable res = foucalc( struct{ time = lc.time, + rate1 = ref.rate - lc.rate, + rate2 = lc.rate}, + dimseg ; normtype=normtype); + + variable r_filtered = filter_fouquan_for_frequency_range(res, f_min, f_max); + + if (qualifier_exists("raw")) list_append(raw, r_filtered); + + variable mom = moment(r_filtered.lag12); + tlags[ii] = mom.ave; + tlagerrs[ii] = mom.sdev/sqrt(mom.num); % standard error on the mean + + } + + if (qualifier_exists("raw")){ + return struct{ lags = tlags, errs = tlagerrs, raw=raw}; + } else { + return struct{ lags = tlags, errs = tlagerrs}; + } +} -- GitLab From de99b20667977eea1750b3b67e11adfcc868db61 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Thu, 21 Jan 2021 10:30:00 +0100 Subject: [PATCH 2/9] Foucalc: Export CPD calculation into new function * Clean-up of foucalc code: export the cross-power-density computation into new function cross_power_density * Also update help of dorDFT: It previously said that the normalization is 1/sqrt(2*PI) but I think this is incorrect. The normalization of the ISIS FFT is 1/sqrt(N) and dorDFT multiplies the output by N, thus the normalization is sqrt(N) --- src/data/timing/cross_power_density.sl | 25 +++++++++++++++++++++++++ src/data/timing/dorDFT.sl | 12 ++++++++---- src/data/timing/foucalc.sl | 18 +++++------------- 3 files changed, 38 insertions(+), 17 deletions(-) create mode 100644 src/data/timing/cross_power_density.sl diff --git a/src/data/timing/cross_power_density.sl b/src/data/timing/cross_power_density.sl new file mode 100644 index 00000000..b1b78dc0 --- /dev/null +++ b/src/data/timing/cross_power_density.sl @@ -0,0 +1,25 @@ +private define cross_power_density(lc1, lc2, numseg, dimseg){ +%!%+ +%\function{cross_power_density} +%\synopsis{Timing Tools: Cross Power Density (aka. Cross-Spectrum)} +%\usage{cpd = cross_power_density(lc1, lc2, numseg, dimseg)} +%\description +% Calculates the CPD, averaged over all segments, as +% Conj(DFT(rate1))*DFT(rate2) +%\seealso{dorDFT} +%!%- + variable startbin = 0; + variable endbin = dimseg-1; + + variable cpd = Complex_Type[dimseg/2]; % = 0 + 0i + loop(numseg) + { + variable cpd_of_segment = Conj(dorDFT(lc1[[startbin:endbin]])) * dorDFT(lc2[[startbin:endbin]]); + cpd += cpd_of_segment; + startbin += dimseg; + endbin += dimseg; + } + + cpd /= (numseg * dimseg); % Average over all segments + return cpd; +} diff --git a/src/data/timing/dorDFT.sl b/src/data/timing/dorDFT.sl index 1ad36291..b7854c15 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 %!%- { variable l = length(rate); diff --git a/src/data/timing/foucalc.sl b/src/data/timing/foucalc.sl index ff465229..ed38cfa9 100644 --- a/src/data/timing/foucalc.sl +++ b/src/data/timing/foucalc.sl @@ -56,6 +56,7 @@ define foucalc() % 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; @@ -189,19 +190,10 @@ 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 + if (verbose) vmessage("Calculating cross power spectrum (=> coherence and timelags) for rate%d and rate%d.", lo, hi); + variable cpd = cross_power_density(array_struct_field(lc, "rate"+lostr), + array_struct_field(lc, "rate"+histr), + numseg, dimseg); variable Re_cpd = Real(cpd); variable Im_cpd = Imag(cpd); set_struct_field(fouquant, "realcpd"+lostr+histr, Re_cpd); -- GitLab From 56dfa43d430c602c27e64316d0e5a1192bf22ed7 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Fri, 22 Jan 2021 13:57:54 +0100 Subject: [PATCH 3/9] Foucalc: Enable use of "rate" (stands for "rate1") The old behavior, i.e. using a struct of the form struct{time, rate1} is still possible, but now also struct{time, rate} can be used. --- src/data/timing/cross_power_density.sl | 6 +++--- src/data/timing/foucalc.sl | 16 ++++++++++------ 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/data/timing/cross_power_density.sl b/src/data/timing/cross_power_density.sl index b1b78dc0..3a904f10 100644 --- a/src/data/timing/cross_power_density.sl +++ b/src/data/timing/cross_power_density.sl @@ -1,8 +1,8 @@ -private define cross_power_density(lc1, lc2, numseg, dimseg){ +private define cross_power_density(rate1, rate2, numseg, dimseg){ %!%+ %\function{cross_power_density} %\synopsis{Timing Tools: Cross Power Density (aka. Cross-Spectrum)} -%\usage{cpd = cross_power_density(lc1, lc2, numseg, dimseg)} +%\usage{cpd = cross_power_density(rate1, rate2, numseg, dimseg)} %\description % Calculates the CPD, averaged over all segments, as % Conj(DFT(rate1))*DFT(rate2) @@ -14,7 +14,7 @@ private define cross_power_density(lc1, lc2, numseg, dimseg){ variable cpd = Complex_Type[dimseg/2]; % = 0 + 0i loop(numseg) { - variable cpd_of_segment = Conj(dorDFT(lc1[[startbin:endbin]])) * dorDFT(lc2[[startbin:endbin]]); + variable cpd_of_segment = Conj(dorDFT(rate1[[startbin:endbin]])) * dorDFT(rate2[[startbin:endbin]]); cpd += cpd_of_segment; startbin += dimseg; endbin += dimseg; diff --git a/src/data/timing/foucalc.sl b/src/data/timing/foucalc.sl index ed38cfa9..8ce3e5c8 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}: @@ -114,9 +115,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 @@ -191,9 +193,11 @@ define foucalc() variable histr = string(hi); if (verbose) vmessage("Calculating cross power spectrum (=> coherence and timelags) for rate%d and rate%d.", lo, hi); - variable cpd = cross_power_density(array_struct_field(lc, "rate"+lostr), - array_struct_field(lc, "rate"+histr), - numseg, dimseg); + + % 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 = cross_power_density(rate1, rate2, numseg, dimseg); variable Re_cpd = Real(cpd); variable Im_cpd = Imag(cpd); set_struct_field(fouquant, "realcpd"+lostr+histr, Re_cpd); -- GitLab From 9da67e678475dad98e657d4c89a3628d5ef9e195 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Fri, 22 Jan 2021 16:12:39 +0100 Subject: [PATCH 4/9] Foucalc: Add error on cross spectrum * The (complex) error on the cross spectrum is the standard error on the mean, i.e. stddev/sqrt(N), computed from the cross-spectra of the individual segments. * The foucalc return struct contains the error as fields errrealcpd and errimagcpd --- src/data/timing/cross_power_density.sl | 23 +++++++++++++++++++++-- src/data/timing/dorDFT.sl | 2 +- src/data/timing/foucalc.sl | 15 +++++++++------ 3 files changed, 31 insertions(+), 9 deletions(-) diff --git a/src/data/timing/cross_power_density.sl b/src/data/timing/cross_power_density.sl index 3a904f10..a4f9734a 100644 --- a/src/data/timing/cross_power_density.sl +++ b/src/data/timing/cross_power_density.sl @@ -2,24 +2,43 @@ private define cross_power_density(rate1, rate2, numseg, dimseg){ %!%+ %\function{cross_power_density} %\synopsis{Timing Tools: Cross Power Density (aka. Cross-Spectrum)} -%\usage{cpd = cross_power_density(rate1, rate2, numseg, dimseg)} +%\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 - return cpd; + + % ====== 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 b7854c15..793cabf1 100644 --- a/src/data/timing/dorDFT.sl +++ b/src/data/timing/dorDFT.sl @@ -10,7 +10,7 @@ private define dorDFT(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 +% bin up to Nyquist frequency %!%- { variable l = length(rate); diff --git a/src/data/timing/foucalc.sl b/src/data/timing/foucalc.sl index 8ce3e5c8..bebfbaff 100644 --- a/src/data/timing/foucalc.sl +++ b/src/data/timing/foucalc.sl @@ -52,6 +52,8 @@ 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 @@ -86,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 @@ -197,11 +199,12 @@ define foucalc() % 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 = cross_power_density(rate1, rate2, numseg, dimseg); - 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); + 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); -- GitLab From b4290fef3406b80ef0753aca4d3e86f033fb748d Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Fri, 22 Jan 2021 16:18:27 +0100 Subject: [PATCH 5/9] segment_lc_for_psd: Print number of segments if verbose --- src/data/timing/segment_lc_for_psd.sl | 31 ++++++++++++++------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/src/data/timing/segment_lc_for_psd.sl b/src/data/timing/segment_lc_for_psd.sl index c1916c49..18ed05bb 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")) { + if (verbose>1){ vmessage(" ***Skipping lightcurve segment with index %i, length %.2fs",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: %.2f s (%.3f percent of total data time=%.2fs)", 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; -- GitLab From 4a64f55c3cb6eb6f13385b0d658abe7f7ea9ef20 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Mon, 25 Jan 2021 11:59:02 +0100 Subject: [PATCH 6/9] Foucalc: return struct in colacal colacal used to return 6 parameters. It is more practical and nicer to read if a struct is returned instead. Also, because there are different ways to calculate the error on the time lag than Nowak+99. If they're added, I struct will prevent the code line from exploding when more and more parameters are added. --- src/data/timing/colacal.sl | 3 ++- src/data/timing/foucalc.sl | 18 ++++++++---------- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/src/data/timing/colacal.sl b/src/data/timing/colacal.sl index 8ab0dbde..14524268 100644 --- a/src/data/timing/colacal.sl +++ b/src/data/timing/colacal.sl @@ -35,5 +35,6 @@ 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; + return struct{rawcof=rawcof, cof=cof, errcof=errcof, + lag=lag, errlag=errlag, sigcpd=sigcpd}; } diff --git a/src/data/timing/foucalc.sl b/src/data/timing/foucalc.sl index bebfbaff..ec2d8ceb 100644 --- a/src/data/timing/foucalc.sl +++ b/src/data/timing/foucalc.sl @@ -217,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; } } -- GitLab From 1c1ac030d8aaa09356d26d4edfaa618fe886ad31 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Tue, 2 Feb 2021 09:49:33 +0100 Subject: [PATCH 7/9] colacal: Add error calc. from Uttley+14 --- src/data/timing/colacal.sl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/data/timing/colacal.sl b/src/data/timing/colacal.sl index 14524268..9a4f6bd9 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,6 +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); + 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, sigcpd=sigcpd}; + lag=lag, errlag=errlag, errlag_uttley14=errlag_uttley14, + sigcpd=sigcpd}; } -- GitLab From 108e0ada2f479da3636193d294294353938c5bba Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Tue, 2 Feb 2021 09:49:56 +0100 Subject: [PATCH 8/9] lag_energy_spectrum: Add outfile, lc_ref qualifier * Via the lc_ref qualifier one can use a broad lightcurve as reference band. If not set, the default is to calculate the reference lightcurve as the sum of all energy-resolved lightcurves but the current channel-of-interest. * Add outfile qualifier: If set, a FITS file is written with information about CPD, phase, averaging, etc. --- src/data/timing/lag_energy_spectrum.sl | 342 +++++++++++++++++++------ src/data/timing/segment_lc_for_psd.sl | 4 +- 2 files changed, 263 insertions(+), 83 deletions(-) diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl index 9379ca11..3afd5fce 100644 --- a/src/data/timing/lag_energy_spectrum.sl +++ b/src/data/timing/lag_energy_spectrum.sl @@ -1,119 +1,299 @@ -require("isisscripts"); +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); +} -define filter_lc_for_timeseg(lc, t_lo, t_hi){ - variable filt = where(t_lo <= lc.time < t_hi); - return struct_filter(lc, 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 filter_fouquan_for_frequency_range(fouquan, f_min, f_max) { - variable filt = where(f_min <= fouquan.freq <= f_max); - return struct_filter(fouquan, filt ; copy); +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}; } -define lag_energy_spectrum(lc_list, lc_ref, dimseg) { +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-energy spectrum between -%energy-resolved lightcurves and a broad reference lightcurve.} +% energy-resolved lightcurves.} %\usage{} %\qualifiers{ -% tstart - The lightcurves can be constrained to specific time -% intervals (e.g., when the lightcurve changes -% dramatically in flux or shows really large gaps). Unit: s -% tstop - See tstart. Unit: s -% dt - The lightcurves are parsed to segment_lc_for_psd and -% split at the gaps. In order to do so, one needs the time -% resolution at is segmented and accounted for the -% gaps. dt is the time resolution of the lightcurve in -% seconds. -% f_min - The lag-energy spectrum can also be computed only in -% a small frequency band (for a frequency-resolved -% lag-energy spectrum). Default: lowest sampled frequency +% 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 -% normtype - Normalization type of PSD [default: Miyamoto]. -% raw - return an additional list of the raw lag-frequency spectra. +% 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) %} %\description % * lc_list: List of names of the energy-resolved lightcurves -% * lc_ref: Name of the (broad) reference lightcurve % * dimseg: Segmentation length in bins, needed to segment the % lightcurve and do the Fourier calculation -% +% % returns: the lag-energy-spectrum as a struct with the fields -% lags, errs +% lag, lagerrs %\notes -% The function does *not* use and will not return 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 default energy lag 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(Im(CPD)), mean(Re(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. 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. +% +% 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 +% +% 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 -%\seealso{segment_lc_for_psd, foucalc} +% %% 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_ref, "TIMEDEL")); - if (dt==NULL) { - return "ERROR: Couldn't receive the time resolution from TIMEDEL lightcurve, need dt qualifier!"; - } - + 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 normtype = qualifier("normtype","Miyamoto"); + variable verbose = qualifier("verbose", 0); + variable normtype = "Miyamoto"; - variable lc_ref_segm = segment_lc_for_psd(fits_read_table(lc_ref), dt, dimseg); - - variable tstart = qualifier("tstart", min(lc_ref_segm.time)); - variable tstop = qualifier("tstop", max(lc_ref_segm.time)); + 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 ref = filter_lc_for_timeseg(lc_ref_segm, tstart, tstop); + variable lags = Double_Type[length(lc_list)]; + variable errlag_nowak99 = Double_Type[length(lc_list)]; + variable errlag_uttley14 = Double_Type[length(lc_list)]; - variable tlags = Float_Type[length(lc_list)]; - variable tlagerrs = Float_Type[length(lc_list)]; + 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("raw")) variable raw = {}; + 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[length(lc_list)]; + variable max_freqs = Double_Type[length(lc_list)]; + variable mean_realcpds = Double_Type[length(lc_list)]; + variable mean_realcpd_errs = Double_Type[length(lc_list)]; + variable mean_imagcpds = Double_Type[length(lc_list)]; + variable mean_imagcpd_errs = Double_Type[length(lc_list)]; + variable mean_phases = Double_Type[length(lc_list)]; + variable errlag_properr = Double_Type[length(lc_list)]; + variable errlag_ingram19 = Double_Type[length(lc_list)]; + variable n_segments = Integer_Type[length(lc_list)]; + variable n_frequencies = Integer_Type[length(lc_list)]; + + } - %% Go through the energy-resolved lightcurves + %% Go through the energy-resolved lightcurves (channel-of-interests) variable ii; _for ii (0,length(lc_list)-1,1){ - variable lc_segm = segment_lc_for_psd(fits_read_table(lc_list[ii]), dt, dimseg); - - variable lc = filter_lc_for_timeseg(lc_segm, tstart, tstop); - - %% Test whether the time interval of reference and energy-resolved - %% lightcurve is the same - if (length(ref.time)!=length(lc.time) - or (ref.time[0]!=lc.time[0]) - or (ref.time[-1]!=lc.time[-1])){ - return sprintf("ERROR: The time arrays of %s and %s seem to be different!", - lc_ref, lc_list[ii]); - } - - %% Calculate the cross-spectrum (CPD) and lags with foucalc Note - %% that we subtract off the small band LC from the reference LC. - %% Uttley+11 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+11 Sect 3.2. - %% The time lag is computed from the CPD as: - %% phase = atan2(r.imagcpd12,r.realcpd12); - %% tlag = phase/(2*PI*r.freq); - variable res = foucalc( struct{ time = lc.time, - rate1 = ref.rate - lc.rate, - rate2 = lc.rate}, + variable channel_of_interest = getLightcurve (lc_list[ii], dt, dimseg ;; __qualifiers); + variable reference = getReferenceLightcurve (lc_ref, channel_of_interest, ii, dt, dimseg ;; __qualifiers); + + variable res = foucalc( struct{ time = channel_of_interest.time, + rate1 = reference.rate, + rate2 = channel_of_interest.rate}, dimseg ; normtype=normtype); + + variable res_f = filterFouQuantForFrequencyRange(res, f_min, f_max ; verbose=verbose); - variable r_filtered = filter_fouquan_for_frequency_range(res, f_min, f_max); - - if (qualifier_exists("raw")) list_append(raw, r_filtered); + %% 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; - variable mom = moment(r_filtered.lag12); - tlags[ii] = mom.ave; - tlagerrs[ii] = mom.sdev/sqrt(mom.num); % standard error on the mean + 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("raw")){ - return struct{ lags = tlags, errs = tlagerrs, raw=raw}; - } else { - return struct{ lags = tlags, errs = tlagerrs}; + + if (qualifier_exists("outfile")==1){ + variable extname = sprintf("BAND_%g_%gHz", f_min, f_max); + variable sdata = struct{ + "ENERGY_LO"=qualifier("elo", Double_Type[length(lc_list)]), + "ENERGY_HI"=qualifier("ehi", Double_Type[length(lc_list)]), + "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 + }; + 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 18ed05bb..27c5b735 100644 --- a/src/data/timing/segment_lc_for_psd.sl +++ b/src/data/timing/segment_lc_for_psd.sl @@ -85,13 +85,13 @@ define segment_lc_for_psd(lc, dt, dimseg) { %{{{ } else { rejected_time += exposure; if (verbose>1){ - vmessage(" ***Skipping lightcurve segment with index %i, length %.2fs",ii,exposure); + vmessage(" ***Skipping lightcurve segment with index %i, length %gs",ii,exposure); } } } if (verbose>0) { - vmessage("\nRejected time: %.2f s (%.3f percent of total data time=%.2fs)", + 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); -- GitLab From b99e47b1806f39428e3b77a81be8014bea57dc53 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Mon, 8 Feb 2021 16:21:29 +0100 Subject: [PATCH 9/9] lag_energy_spectrum: add RMS to outfile --- src/data/timing/lag_energy_spectrum.sl | 86 ++++++++++++++------------ 1 file changed, 47 insertions(+), 39 deletions(-) diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl index 3afd5fce..3f2491c4 100644 --- a/src/data/timing/lag_energy_spectrum.sl +++ b/src/data/timing/lag_energy_spectrum.sl @@ -70,9 +70,15 @@ private define getReferenceLightcurve(lc_ref, cof_lc, index, dt, dimseg){ define lag_energy_spectrum(lc_list, dimseg) { %!%+ %\function{lag_energy_spectrum} -%\synopsis{Calculates the time lag-energy spectrum between -% energy-resolved lightcurves.} -%\usage{} +%\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 @@ -92,17 +98,9 @@ define lag_energy_spectrum(lc_list, dimseg) { % % (Furthermore, the qualifiers of segment_lc_for_psd can be specified) %} -%\description -% * lc_list: List of names of the energy-resolved lightcurves -% * dimseg: Segmentation length in bins, needed to segment the -% lightcurve and do the Fourier calculation -% -% returns: the lag-energy-spectrum as a struct with the fields -% lag, lagerrs %\notes -% * The default energy lag error computation is Nowak et al., 1999, -% ApJ, 510, 874 (Sect. 4), i.e., from coherence without -% bias factor. +% * 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. @@ -110,7 +108,7 @@ define lag_energy_spectrum(lc_list, dimseg) { % * 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(Im(CPD)), mean(Re(CPD)))/(PI*(f_min+f_max)) +% 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). @@ -123,9 +121,10 @@ define lag_energy_spectrum(lc_list, dimseg) { % 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. 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. +% 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 @@ -145,9 +144,12 @@ define lag_energy_spectrum(lc_list, dimseg) { % 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 +% 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. % @@ -179,9 +181,12 @@ define lag_energy_spectrum(lc_list, dimseg) { if (verbose>0) vmessage("Frequency range %g-%g Hz", f_min, f_max); - variable lags = Double_Type[length(lc_list)]; - variable errlag_nowak99 = Double_Type[length(lc_list)]; - variable errlag_uttley14 = Double_Type[length(lc_list)]; + 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){ @@ -197,32 +202,33 @@ define lag_energy_spectrum(lc_list, dimseg) { } else { fp = qualifier("outfile"); } - variable min_freqs = Double_Type[length(lc_list)]; - variable max_freqs = Double_Type[length(lc_list)]; - variable mean_realcpds = Double_Type[length(lc_list)]; - variable mean_realcpd_errs = Double_Type[length(lc_list)]; - variable mean_imagcpds = Double_Type[length(lc_list)]; - variable mean_imagcpd_errs = Double_Type[length(lc_list)]; - variable mean_phases = Double_Type[length(lc_list)]; - variable errlag_properr = Double_Type[length(lc_list)]; - variable errlag_ingram19 = Double_Type[length(lc_list)]; - variable n_segments = Integer_Type[length(lc_list)]; - variable n_frequencies = Integer_Type[length(lc_list)]; - + 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,length(lc_list)-1,1){ + _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); - + 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 @@ -265,8 +271,8 @@ define lag_energy_spectrum(lc_list, dimseg) { if (qualifier_exists("outfile")==1){ variable extname = sprintf("BAND_%g_%gHz", f_min, f_max); variable sdata = struct{ - "ENERGY_LO"=qualifier("elo", Double_Type[length(lc_list)]), - "ENERGY_HI"=qualifier("ehi", Double_Type[length(lc_list)]), + "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, @@ -280,7 +286,9 @@ define lag_energy_spectrum(lc_list, dimseg) { "N_SEGMENTS"=n_segments, "N_FREQUENCIES"=n_frequencies, "MIN_FREQUENCY"=min_freqs, - "MAX_FREQUENCY"=max_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 -- GitLab