From 94e47e049e42c4ef82fd578ff810957d0c9199a3 Mon Sep 17 00:00:00 2001 From: pupperemeritus Date: Sat, 10 Jun 2023 19:40:25 +0530 Subject: [PATCH] Input of lsft_slow,fast changed to be non stingray specific, added changelog, --- docs/changes/737.feature.rst | 2 ++ stingray/fourier.py | 28 ++++++++++++++-------------- stingray/lombscargle.py | 26 +++++++++++--------------- 3 files changed, 27 insertions(+), 29 deletions(-) create mode 100644 docs/changes/737.feature.rst diff --git a/docs/changes/737.feature.rst b/docs/changes/737.feature.rst new file mode 100644 index 000000000..534f4c495 --- /dev/null +++ b/docs/changes/737.feature.rst @@ -0,0 +1,2 @@ +- Implemented the Lomb Scargle Fourier Transform +- Using which wrote the corresponding :class:`LombScargleCrossspectrum` and :class:`LombScarglePowerspectrum` \ No newline at end of file diff --git a/stingray/fourier.py b/stingray/fourier.py index f39ccd426..bb6b3574f 100644 --- a/stingray/fourier.py +++ b/stingray/fourier.py @@ -4,6 +4,7 @@ from typing import Optional import numpy as np +import numpy.typing as npt from astropy.table import Table from astropy.timeseries.periodograms.lombscargle.implementations.utils import trig_sum @@ -1928,11 +1929,8 @@ def avg_cs_from_events( return results -def lsft_fast(lc: Lightcurve, w0, df, Nf): +def lsft_fast(y, t, dy, w0, df, Nf): # Work in Progress - y = lc.counts - t = lc.time - dy = lc.counts_err weights = dy**-2.0 weights /= weights.sum() @@ -1941,8 +1939,9 @@ def lsft_fast(lc: Lightcurve, w0, df, Nf): def lsft_slow( - lc: Lightcurve, - ww: np.ndarray, + y: npt.ArrayLike, + t: npt.ArrayLike, + ww: npt.ArrayLike, sign: Optional[int] = 1, fullspec: Optional[bool] = False, ): @@ -1951,8 +1950,11 @@ def lsft_slow( Parameters ---------- - lc : :class:`stingray.lightcurve.Lightcurve` - A light curve object. + y : a `:class:numpy.array` of floats + Observations to be transformed. + + y : `:class:numpy.array` of floats + Times of the observations freqs : numpy.ndarray An array of frequencies at which the transform is sampled. @@ -1974,11 +1976,9 @@ def lsft_slow( const1 = 1 / np.sqrt(2) const2 = const1 * sign - xx = lc.counts - sum_xx = np.sum(xx) - t = lc.time + sum_xx = np.sum(y) - num_xt = len(xx) + num_xt = len(y) num_ww = len(ww) ft_real = ft_imag = np.zeros((num_ww)) @@ -1997,8 +1997,8 @@ def lsft_slow( watan = np.arctan2(ssum, csum) wtau = 0.5 * watan - sumr = np.sum(np.multiply(xx, np.cos(wrun * t - wtau))) - sumi = np.sum(np.multiply(xx, np.sin(wrun * t - wtau))) + sumr = np.sum(np.multiply(y, np.cos(wrun * t - wtau))) + sumi = np.sum(np.multiply(y, np.sin(wrun * t - wtau))) scos2 = np.sum((np.power(np.cos(wrun * t - wtau), 2))) ssin2 = np.sum((np.power(np.sin(wrun * t - wtau), 2))) diff --git a/stingray/lombscargle.py b/stingray/lombscargle.py index eb0ab6e4b..b3ed92d91 100644 --- a/stingray/lombscargle.py +++ b/stingray/lombscargle.py @@ -2,17 +2,15 @@ from typing import Optional, Union import numpy as np - +import numpy.typing as npt from astropy.timeseries.periodograms import LombScargle -from stingray.crossspectrum import Crossspectrum -from stingray.exceptions import StingrayError -from stingray.utils import simon - -from .events import EventList # Convert to relative nomenclature -from .lightcurve import Lightcurve # Convert to relative nomenclature - -from .fourier import lsft_slow, lsft_fast +from .crossspectrum import Crossspectrum +from .events import EventList +from .exceptions import StingrayError +from .fourier import lsft_fast, lsft_slow +from .lightcurve import Lightcurve +from .utils import simon # method default will change to fast after implementation of the fast algorithm @@ -263,8 +261,8 @@ def _make_crossspectrum(self, lc1, lc2, fullspec=False, method="fast"): if lc1.n != lc2.n: raise StingrayError("Lightcurves do not have the same number of bins per segment.") - # if not np.isclose(lc1.dt, lc2.dt, rtol=0.1 * lc1.dt / lc1.tseg): - # raise StingrayError("Lightcurves do not have the same time binning dt.") + if not np.isclose(lc1.dt, lc2.dt, rtol=0.1 * lc1.dt / lc1.tseg): + raise StingrayError("Lightcurves do not have the same time binning dt.") lc1.dt = lc2.dt self.dt = lc1.dt @@ -290,14 +288,12 @@ def _make_crossspectrum(self, lc1, lc2, fullspec=False, method="fast"): ) if self.__class__.__name__ == "LombScarglePowerspectrum": - print("ps") self.power_err = self.unnorm_power_err = self.power / np.sqrt(self.m) elif self.__class__.__name__ == "LombScargleCrossspectrum": simon( "Errorbars on cross spectra are not thoroughly tested." "Please report any inconsistencies." ) - print("ls") self.unnorm_power_err = np.sqrt(2) / np.sqrt(self.m) self.unnorm_power_err /= np.divide(2, np.sqrt(np.abs(self.nphots1 * self.nphots2))) self.unnorm_power_err += np.zeros_like(self.unnorm_power) @@ -389,8 +385,8 @@ def _ls_cross(self, lc1, lc2, ww=None, fullspec=False, method="slow"): ww = ww2 if method == "slow": - lsft1 = lsft_slow(lc1, ww, sign=1, fullspec=fullspec) - lsft2 = lsft_slow(lc2, ww, sign=-1, fullspec=fullspec) + lsft1 = lsft_slow(lc1.counts, lc1.time, ww, sign=1, fullspec=fullspec) + lsft2 = lsft_slow(lc2.counts, lc2.time, ww, sign=-1, fullspec=fullspec) # elif method == "fast": # lsft1 = lsft_fast(lc1, ww, fullspec=fullspec) cross = np.multiply(lsft1, lsft2)