diff --git a/src/data/timing/colacal.sl b/src/data/timing/colacal.sl index a865f53362bb9e623fb032d8160c2ac5767dc678..53437a58f2a941db49ce37a5ac90ae5d789ed0d2 100644 --- a/src/data/timing/colacal.sl +++ b/src/data/timing/colacal.sl @@ -23,7 +23,6 @@ define colacal(freq, cpd, noicpd, lopsd, hipsd, noilopsd, noihipsd, siglopsd, si % 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 %!%- { @@ -36,9 +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); - variable rawcof_w_bias = (sqr(abs(cpd)) - noicpd) / (lopsd * hipsd); % Uttley+14 Eq. 11, note that this divides by psd and not sigpsd - 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}; + lag=lag, errlag=errlag, sigcpd=sigcpd}; } diff --git a/src/data/timing/foucalc.sl b/src/data/timing/foucalc.sl index 1613ec0a6dae53c6e538180ae384d4ffaac8bfb9..a1a4588e81c610e9c4f84fe1a54586260ebafda1 100644 --- a/src/data/timing/foucalc.sl +++ b/src/data/timing/foucalc.sl @@ -14,8 +14,10 @@ define foucalc() %\qualifier{deadtime}{ [\code{=1e-5}]: detector deadtime in seconds} %\qualifier{fmin}{: minimum frequency used for RMS calculation} %\qualifier{fmax}{: maximum frequency used for RMS calculation} -%\qualifier{RMS}{: reference to a variable to store the RMS of each light curve -% in the [fmin, fmax] band. Can only be used with normtype="Miyamoto".} +%\qualifier{RMS}{: DEPRECATED, use rms qualifier instead.} +%\qualifier{rms}{: reference to a variable to store the signal, noise RMS and error +% of each light curve in the [fmin, fmax] band. Can only be used with +% normtype="Miyamoto".} %\qualifier{avgrate}{: reference to a variable to store the average rate of each light curve} %\qualifier{compact}{: compact the output structure, only keep the most important quantities} %\qualifier{noCPD}{: do not calculate cross power densities and derived quantities} @@ -26,7 +28,8 @@ 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. +% Also specifying "rate" instead of "rate1" is possible if no CPD +% should be computed. % % \code{dimseg} is the segment size used for the FFTs, which should therefore be a power of 2. % @@ -76,20 +79,27 @@ define foucalc() variable numseg = int(dimlc/dimseg); variable tseg = bintime * dimseg; variable timelc = bintime * dimlc; - variable numchan = length(get_struct_field_names(lc[0])) - 1; % lc[0] = lc, if lc is no array + variable str_fields = get_struct_field_names(lc[0]); % lc[0]=lc, if lc is no array but struct + variable numchan = length(str_fields[where(is_substr(str_fields, "rate"))]); % Lightcurve parameters variable avgrate = Double_Type[numchan]; - variable rms = Double_Type[numchan]; - + variable rms_deprecated = Double_Type[numchan]; % used for deprecated RMS qualifier + variable fieldnames = ["freq", "numavgall"]; % field names of output structure + variable rms_fieldnames = String_Type[0]; + variable i, j; - _for i (1, numchan, 1) + _for i (1, numchan, 1){ fieldnames = [fieldnames, ["rawpsd", "errpsd", "normpsd", "errnormpsd", "noipsd", "sigpsd", "noinormpsd", "signormpsd", "effnoinormpsd"]+string(i)]; + rms_fieldnames = [rms_fieldnames, ["sigrms", "errrms", "noirms"]+string(i)]; + } + _for i (1, numchan, 1) _for j (i+1, numchan, 1) 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 + variable rmsquant = @Struct_Type(rms_fieldnames); % Calculating Fourier Frequency Array fouquant.numavgall = ones(dimseg/2) * numseg; @@ -114,7 +124,7 @@ define foucalc() vmessage("warning (foucalc): Deadtime keyword not set. Using default for RXTE: 1e-05 s"); } - % RMS Qualifiers + % RMS Qualifiers (see, e.g., Uttley et al., A&A Review 22, 72, 2014, Sect. 2.2.3) variable fmin = qualifier("fmin", min(fouquant.freq)); variable fmax = qualifier("fmax", max(fouquant.freq)); @@ -157,14 +167,12 @@ define foucalc() % Calculating observational noise with deadtime influence (noipsd, noinormpsd) variable noipsd, noinormpsd; -% (, noipsd) = psdcorr_zhang(avgrate[i-1], tseg, dimseg;; struct_combine(__qualifiers(), struct { avgbkg=avgbkg[i-1], unnormalized })); -% (, noinormpsd) = psdcorr_zhang(avgrate[i-1], tseg, dimseg;; __qualifiers()); (, noipsd, noinormpsd) = psdcorr_zhang(avgrate[i-1], tseg, dimseg;; struct_combine(__qualifiers(), struct { avgbkg=avgbkg[i-1] })); set_struct_field(fouquant, "noipsd"+istr, noipsd); set_struct_field(fouquant, "noinormpsd"+istr, noinormpsd); % effective noise level - set_struct_field(fouquant, "effnoinormpsd"+istr, noinormpsd/sqrt(numseg)); % fouquant.numavgall + set_struct_field(fouquant, "effnoinormpsd"+istr, noinormpsd/sqrt(numseg)); % calculate normalized signal psd variable signormpsd = tempnormpsd - noinormpsd; set_struct_field(fouquant, "signormpsd"+istr, signormpsd); @@ -175,8 +183,15 @@ define foucalc() variable idxmin = min(idx); variable idxmax = max(idx); variable df = fouquant.freq[[idxmin+1:idxmax]] - fouquant.freq[[idxmin:idxmax-1]]; - rms[i-1] = sqrt( sum(signormpsd[[idxmin:idxmax-1]] * df) ); - if(verbose) vmessage("RMS(%g-%g Hz) = %.1f%%", fmin, fmax, rms[i-1]*100); + rms_deprecated[i-1] = sqrt( sum(signormpsd[[idxmin:idxmax-1]] * df) ); + + variable sigrms = sqrt( sum(signormpsd[[idxmin:idxmax-1]] * df) ); + variable noirms = sqrt( sum(noinormpsd[[idxmin:idxmax-1]] * df) ); + variable errrms = sqrt((2.*sigrms^2*noirms^2+noirms^4)/(2.*numseg*length(idx)*sigrms^2)); + set_struct_field(rmsquant, "sigrms"+istr, sigrms); + set_struct_field(rmsquant, "errrms"+istr, errrms); + set_struct_field(rmsquant, "noirms"+istr, noirms); + if(verbose) vmessage("RMS(%g-%g Hz) = (%g +/- %g) %%", fmin, fmax, sigrms*100, errrms*100); } else if(verbose) message(""); % to finish the earlier printf without "\n" @@ -200,11 +215,10 @@ define foucalc() 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); + (cpd, errcpd) = cross_power_density(array_struct_field(lc, "rate"+lostr), + array_struct_field(lc, "rate"+histr), + 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)); @@ -227,19 +241,30 @@ define foucalc() 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); cola= NULL; cpd = NULL; errcpd=NULL; } } + variable rmsref = qualifier("rms"); + if(normtype=="Miyamoto" and rmsref!=NULL) + { + if(typeof(rmsref) == Ref_Type) + (@rmsref) = rmsquant; + else + vmessage("warning (%s): rms qualifier has to be a reference", _function_name()); + } + + % Deprecated RMS qualifier (kept for backwards-compatibility) variable RMSref = qualifier("RMS"); if(normtype=="Miyamoto" and RMSref!=NULL) { - if(typeof(RMSref)==Ref_Type) - (@RMSref) = rms; - else + if(typeof(RMSref)==Ref_Type){ + vmessage("warning (%s): RMS qualifier is deprecated, please change to rms.", _function_name()); + (@RMSref) = rms_deprecated; + } else { vmessage("warning (%s): RMS qualifier has to be a reference", _function_name()); + } } variable avgrateref = qualifier("avgrate"); diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl index 79fef22a5023eb3566127149c46918e11ba2d101..851a94a6fa49f2aa4fcfaf0e9f25c56395a7ef5e 100644 --- a/src/data/timing/lag_energy_spectrum.sl +++ b/src/data/timing/lag_energy_spectrum.sl @@ -1,14 +1,20 @@ % -*- mode: slang; mode: fold -*- -private define filterFouQuantForFrequencyRange(fouquan, f_min, f_max) { +private define filterFouQuantForFrequencyRange(fouquant, fmin, fmax) { %{{{ - variable filt = where(f_min < fouquan.freq <= f_max); - if (qualifier("verbose", 0) > 1){ - vmessage(" ***(lag_energy_spectrum): Minimal frequency bin %d at %g Hz", filt[0], fouquan.freq[filt[0]]); - vmessage(" ***(lag_energy_spectrum): Maximal frequency bin %d at %g Hz", filt[-1], fouquan.freq[filt[-1]]); - vmessage(" ***(lag_energy_spectrum): Number of frequency bins: %d", length(filt)); - } - return struct_filter(fouquan, filt ; copy); + variable idx = where(fmin <= fouquant.freq <= fmax); + variable idxmin = min(idx); + variable idxmax = max(idx); + variable df = fouquant.freq[[idxmin+1:idxmax]] - fouquant.freq[[idxmin:idxmax-1]]; + + if (length(idx) == 0) vmessage("ERROR (lag_energy_spectrum): No frequency bins between %g-%g in Fourier products!", fmin, fmax); + if (qualifier("verbose", 0) > 1){ + vmessage(" ***(lag_energy_spectrum): Minimal frequency bin %d at %g Hz", idxmin, fouquant.freq[idxmin]); + vmessage(" ***(lag_energy_spectrum): Maximal frequency bin %d at %g Hz", idxmax, fouquant.freq[idxmax]); + vmessage(" ***(lag_energy_spectrum): Number of frequency bins: %d", length(idx)); + } + variable ret = struct_filter(fouquant, idx; copy); + return struct_combine(ret, struct{df = df, idxmin = idxmin, idxmax = idxmax, n_freqs = length(idx)}); } %}}} private define verifyLightcurveTime(lc1, lc2){ @@ -75,11 +81,11 @@ private define getReferenceLightcurve(lc_ref, cof_lc, index, dt, dimseg){ return struct{time=lc_ref.time, rate=rate}; } %}}} -private define write_cpd_to_ascii(realcpd, realcpd_err, imagcpd, imagcpd_err, f_min, f_max){ +private define write_cpd_to_ascii(realcpd, realcpd_err, imagcpd, imagcpd_err, fmin, fmax){ %{{{ % ASCII files of the cross-spectrum are needed for fitting with RELTRANS - variable fptr_imag = fopen(sprintf("freq_cross_spec_imag_%f_%fHz.ascii", f_min, f_max), "w"); - variable fptr_real = fopen(sprintf("freq_cross_spec_real_%f_%fHz.ascii", f_min, f_max), "w"); + variable fptr_imag = fopen(sprintf("freq_cross_spec_imag_%f_%fHz.ascii", fmin, fmax), "w"); + variable fptr_real = fopen(sprintf("freq_cross_spec_real_%f_%fHz.ascii", fmin, fmax), "w"); variable ii; _for ii (0, length(realcpd)-1, 1){ () = fputs(sprintf("%i\t%f\t%f\n", ii+1, imagcpd[ii], imagcpd_err[ii]), fptr_imag); @@ -89,6 +95,84 @@ private define write_cpd_to_ascii(realcpd, realcpd_err, imagcpd, imagcpd_err, f_ () = fclose(fptr_real); } %}}} +private define verify_frequency_binning(dat){ + %{{{ + %% Check that the frequency binning of every energy-resolved + %% lightcurve was the same for the individual frequency bins. If + %% this is not the case this means that the energy-resolved + %% lightcurves were extracted with different binning/time + %% resolution. + variable ii; + variable N_energies = array_shape(dat.lag)[0]; + _for ii (1, N_energies-1, 1){ + if (_eqs(dat.n_segments[ii,*], dat.n_segments[0,*]) == 0 or + _eqs(dat.n_frequencies[ii,*], dat.n_frequencies[0,*]) == 0 or + _eqs(dat.f_min_dft[ii,*], dat.f_min_dft[0,*]) == 0 or + _eqs(dat.f_max_dft[ii,*], dat.f_max_dft[0,*]) == 0){ + vmessage("WARNING (%s): The DFT binning of the energy-resolved lightcurves differs! Check the time resolutios of your individual lightcurves!", + _function_name()); + } + } +} %}}} + +private define write_to_outfile(dat, f_lo, f_hi, keywords, lc_list){ + %{{{ + variable fp = fits_open_file(qualifier("outfile"), "c"); + + variable N_energies, N_les; + (N_energies, N_les) = array_shape(dat.lag)[0], array_shape(dat.lag)[1]; + + variable jj; + _for jj (0, N_les-1, 1){ + % write_cpd_to_ascii(dat.mean_realcpds[*,jj], dat.mean_realcpd_errs[*,jj], + % dat.mean_imagcpds[*,jj], dat.mean_imagcpd_errs[*,jj], + % f_lo[jj], f_hi[jj]); + + variable extname = sprintf("BAND_%f_%fHz", f_lo[jj], f_hi[jj]); + variable sdata = struct{ + "ENERGY_LO" = qualifier("elo", Double_Type[N_energies]), + "ENERGY_HI" = qualifier("ehi", Double_Type[N_energies]), + "LAG" = dat.lag[*,jj], + "LAG_ERR" = dat.errlag[*,jj], + % "LAG_ERR_PROPERR" = dat.errlag_properr[*,jj], + % "LAG_ERR_INGRAM19" = dat.errlag_ingram19[*,jj], + "MEAN_IMAG_CPD" = dat.mean_imagcpds[*,jj], + "MEAN_IMAG_CPD_ERR" = dat.mean_imagcpd_errs[*,jj], + "MEAN_REAL_CPD" = dat.mean_realcpds[*,jj], + "MEAN_REAL_CPD_ERR" = dat.mean_realcpd_errs[*,jj], + "MEAN_PHASE" = dat.mean_phases[*,jj], + "RMS" = dat.rms[*,jj], + "RMS_ERR" = dat.rms_err[*,jj], + }; + + variable keys = struct_combine(struct{ + F_MIN=f_lo[jj], F_MAX=f_hi[jj], N_SEGMENTS=dat.n_segments[0,jj], + N_FREQUENCIES=dat.n_frequencies[0,jj], F_MIN_DFT=dat.f_min_dft[0,jj], + F_MAX_DFT=dat.f_max_dft[0,jj]}, keywords); + variable hist = struct{ + history = [sprintf("Wrote lag-energy spectrum in frequency range %f-%fHz", f_lo[jj], f_hi[jj])], + 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); + } + + %% Write power spectrum into "PSD" extension + variable colnames = String_Type[N_energies+1]; + colnames[0] = "freq"; + variable ii; + _for ii (0, N_energies-1, 1) colnames[ii+1] = sprintf("psd%i", ii); + variable psddata = @Struct_Type(colnames); + set_struct_field(psddata, "freq", dat.freq[0]); % frequency grids are all the same (because LCs have same binning) + _for ii (0, N_energies-1, 1) set_struct_field(psddata, colnames[ii+1], dat.psd[ii]); + fits_write_binary_table(fp, "PSD", psddata); + + fits_close_file(fp); + () = system(sprintf("fthedit %s\[0] keyword=N_LES operation=add value=%i comment=\"Number of frequency-resolved lag-energy spectra\"", + qualifier("outfile"), N_les)); + vmessage("Wrote lag-energy spectrum data to %s", qualifier("outfile")); +} %}}} + define lag_energy_spectrum(lc_list, dimseg) { %!%+ @@ -110,21 +194,20 @@ define lag_energy_spectrum(lc_list, dimseg) { % dt - Time resolution of the lightcurve in seconds % [default: TIMEDEL keyword] % deadtime - Detector deadtime in seconds [default: 0.0s] -% 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 +% f_lo - Array_Type. For frequency-resolved lag-energy +% spectra. Default: lowest sampled frequency [an array +% containing only fmin=1/(dt*dimseg)]. Unit: Hz +% f_hi - See f_lo. 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) +% outfile - String_Type. If set, write various timing products +% into a FITS file. % 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. +% 874 (Sect. 4). % * 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. @@ -132,7 +215,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(Imag(CPD)), mean(Real(CPD)))/(PI*(f_min+f_max)) +% atan2(mean(Imag(CPD)), mean(Real(CPD)))/(PI*(fmin+fmax)) % * 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). @@ -141,192 +224,165 @@ define lag_energy_spectrum(lc_list, dimseg) { % 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: +% 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 extension, +% called BAND__Hz, 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 +% 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*(fmin+fmax)) +% LAG_ERR - Error computed as in Nowak et al., 1999, ApJ, 510, 874 Eq. 16 +% LAG_ERR_PROPERR - Error computed from Gaussian error propagation +% 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. +% real/imaginary segment-averaged cross-spectrum +% RMS/RMS_ERR - Root Mean Squared variability of the +% channel-of-interest lightcurve +% +% Following keywords are written into each extension: % +% N_SEGMENTS - Number of segments averaged over (foucalc.numavgall) +% N_FREQUENCIES - Number of frequencies averaged over (as in +% "mean of the CPD") +% F_MIN/F_MAX - Minimum and maximum frequency as given by +% f_lo, f_hi qualifiers +% F_MIN_DFT/F_MAX_DFT - Minimal/Maximal frequency of the frequency-filtered +% Fourier products - this can be different from +% F_MIN/F_MAX due to the discretization of the +% *Discrete* Fourier Transform +% % 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) +% (f_lo, f_hi) = 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); +% variable les = lag_energy_spectrum(lc_names, dimseg +% ; f_lo = f_lo, f_hi=f_hi, +% verbose=2, outfile="test.fits", elo=elo, ehi=ehi); +% _for ii (0, n_freqs-1, 1) { +% ohplot_with_err(elo, ehi, les.lag[*,ii], les.errlag[*,ii]); % } -% 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 f_lo = qualifier("f_lo", [1/(dt*dimseg)]); + variable f_hi = qualifier("f_hi", [1/(2*dt)]); % Nyquist frequency variable verbose = qualifier("verbose", 0); variable deadtime = qualifier("deadtime", 0.0); variable normtype = "Miyamoto"; if (dt == NULL) return "ERROR: Couldn't receive the time resolution from TIMEDEL, need dt qualifier!"; - if (verbose > 0) vmessage("Calculating lag-energy spectrum in 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 N_les = length(f_lo); + variable dat = struct{ + lag = Double_Type[N_energies, N_les], + errlag = Double_Type[N_energies, N_les], + rms = Double_Type[N_energies, N_les], + rms_err = Double_Type[N_energies, N_les], + min_freqs = Double_Type[N_energies, N_les], + max_freqs = Double_Type[N_energies, N_les], + mean_realcpds = Double_Type[N_energies, N_les], + mean_realcpd_errs = Double_Type[N_energies, N_les], + mean_imagcpds = Double_Type[N_energies, N_les], + mean_imagcpd_errs = Double_Type[N_energies, N_les], + mean_phases = Double_Type[N_energies, N_les], + errlag_properr = Double_Type[N_energies, N_les], + % errlag_ingram19 = Double_Type[N_energies, N_les], + freq = Array_Type[N_energies], + psd = Array_Type[N_energies], + n_segments = Integer_Type[N_energies, N_les], + n_frequencies = Integer_Type[N_energies, N_les], + f_min_dft = Double_Type[N_energies, N_les], + f_max_dft = Double_Type[N_energies, N_les] + }; + variable lc_ref; if (qualifier_exists("lc_ref") == 1){ - lc_ref = getLightcurve(qualifier("lc_ref"), dt, dimseg ;; __qualifiers); + 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){ - if (verbose > 1) vmessage(" ***(%s): Calculate time lags for %s", _function_name(), lc_list[ii]); + %% Loop through energy-resolved lightcurves (channel-of-interests) + variable ii, jj; + _for ii (0, N_energies-1, 1){ + + if (verbose > 0) vmessage(" ***(%s): Calculate time lags for %s", _function_name(), lc_list[ii]); variable channel_of_interest = getLightcurve (lc_list[ii], dt, dimseg ;; __qualifiers); - variable reference = getReferenceLightcurve (lc_ref, channel_of_interest, 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, deadtime=deadtime); - (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; + variable fouquant = foucalc(struct{time = channel_of_interest.time, + rate1 = reference.rate, + rate2 = channel_of_interest.rate}, + dimseg; normtype=normtype, deadtime=deadtime); + dat.freq[ii] = fouquant.freq; + dat.psd[ii] = fouquant.signormpsd2; - 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); + %% Loop through frequency ranges + _for jj (0, N_les-1, 1){ + if (verbose > 1) vmessage(" ***(%s): Calculating lag-energy spectrum in frequency range %g-%g Hz", + _function_name(), f_lo[jj], f_hi[jj]); + variable fouquant_f = filterFouQuantForFrequencyRange(fouquant, f_lo[jj], f_hi[jj]; verbose=verbose); + + % RMS calculation (same as in foucalc but we want it for every + % energy lightcurve and frequency range - running foucalc + % everytime drastically increases the runtime due to the DFT!) + variable sigrms = sqrt( sum(fouquant.signormpsd2[[fouquant_f.idxmin:fouquant_f.idxmax-1]] * fouquant_f.df) ); + variable noirms = sqrt( sum(fouquant.noinormpsd2[[fouquant_f.idxmin:fouquant_f.idxmax-1]] * fouquant_f.df) ); + variable errrms = sqrt((2.*sigrms^2*noirms^2+noirms^4)/(2.*fouquant_f.numavgall[0]*fouquant_f.n_freqs*sigrms^2)); + + dat.rms[ii,jj] = sigrms; + dat.rms_err[ii,jj] = errrms; + + %% Coherence and lag calculation on averaged cross-spectrum + variable cola = colacal(0.5*(f_lo[jj]+f_hi[jj]), % mid frequency + mean(fouquant_f.realcpd12)+mean(fouquant_f.imagcpd12)*1i, + mean(fouquant_f.noicpd12), % not needed for lag+error + mean(fouquant_f.rawpsd1), mean(fouquant_f.rawpsd2), + mean(fouquant_f.noipsd1), mean(fouquant_f.noipsd2), + mean(fouquant_f.sigpsd1), mean(fouquant_f.sigpsd2), % not needed for lag+err + fouquant_f.numavgall[0]*length(fouquant_f.freq)); % K*M with M=number of frequencies averaged over and K=number of segments + (dat.lag[ii,jj], dat.errlag[ii,jj]) = cola.lag, cola.errlag; + + variable m_ReCPD = mean(fouquant_f.realcpd12); + variable m_ReCPD_err = mean(fouquant_f.errrealcpd12); + variable m_ImCPD = mean(fouquant_f.imagcpd12); + variable m_ImCPD_err = mean(fouquant_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))* + dat.errlag_properr[ii,jj]=1/(PI*(f_lo[jj]+f_hi[jj]))* 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 - ); + (m_ImCPD/(m_ImCPD^2+m_ReCPD^2))^2*m_ReCPD_err^2); - %% Implement error calculation from Ingram 2019, Eq. 19 + %% Error calculation from Ingram 2019, Eq. 19: To be implemented - 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); + dat.mean_realcpds[ii,jj] = m_ReCPD; + dat.mean_realcpd_errs[ii,jj] = m_ReCPD_err; + dat.mean_imagcpds[ii,jj] = m_ImCPD; + dat.mean_imagcpd_errs[ii,jj] = m_ImCPD_err; + dat.mean_phases[ii,jj] = m_phase; + dat.f_min_dft[ii,jj] = min(fouquant_f.freq); + dat.f_max_dft[ii,jj] = max(fouquant_f.freq); + dat.n_segments[ii,jj] = fouquant_f.numavgall[0]; + dat.n_frequencies[ii,jj] = length(fouquant_f.freq); } } - if (qualifier_exists("outfile") == 1){ - % write_cpd_to_ascii(mean_realcpds, mean_realcpd_errs, mean_imagcpds, mean_imagcpd_errs, f_min, f_max); - - variable extname = sprintf("BAND_%f_%fHz", 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_UTTLEY14"=errlag_uttley14, - "LAG_ERR_NOWAK99"=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); + verify_frequency_binning(dat); + + if (qualifier_exists("outfile")){ + variable keywords = struct{DIMSEG=dimseg, TIMEDEL=dt, NORMTYPE=normtype}; + write_to_outfile(dat, f_lo, f_hi, keywords, lc_list;; __qualifiers()); } - return struct{ lag = lags, errlag = errlag_nowak99 }; + return dat; }