Skip to content

Commit

Permalink
Input of lsft_slow,fast changed to be non stingray specific, added ch…
Browse files Browse the repository at this point in the history
…angelog,
  • Loading branch information
pupperemeritus committed Jun 24, 2023
1 parent 137e22b commit 94e47e0
Show file tree
Hide file tree
Showing 3 changed files with 27 additions and 29 deletions.
2 changes: 2 additions & 0 deletions docs/changes/737.feature.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
- Implemented the Lomb Scargle Fourier Transform
- Using which wrote the corresponding :class:`LombScargleCrossspectrum` and :class:`LombScarglePowerspectrum`
28 changes: 14 additions & 14 deletions stingray/fourier.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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()

Expand All @@ -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,
):
Expand All @@ -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.
Expand All @@ -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))
Expand All @@ -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)))
Expand Down
26 changes: 11 additions & 15 deletions stingray/lombscargle.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit 94e47e0

Please sign in to comment.