1
$\begingroup$

I am using the astropy Lomb-Scargle periodogram here, https://docs.astropy.org/en/stable/timeseries/lombscargle.html to get the best period fit using its model() method for simulated data given in their examples. I am then plotting this fit over the phase-folded data using the pyasl.foldAt() method but the period fits are not within the range of 0 - 1 for periods above 1d?

When the periodicity within the simulated data is set to 1d, the best period fit works perfectly as shown below; Best_period_fit_1d

But if the period within the simulated data is set to 1.2d, the best period fit is not contained within the 0 - 1 phase range even though the t_fit array given to the ls.model() method is from 0 - 1, plot shown below; Best_period_fit_1.2d

Below is the code I used.

%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
from PyAstronomy import pyasl
from astropy.timeseries import LombScargle

# Quick benchmark of this function on a simple simulated and noisy periodic data!

period = 1.2 # periodicity within the simulated data
rand = np.random.default_rng(42)
t = 100 * rand.random(100)
dy = 0.1 * (1 + rand.random(100))
y = np.sin(2 * np.pi * t * (1/period)) + dy * rand.standard_normal(100)

frequency, power = LombScargle(t, y, dy).autopower()
print('Period at max. power: {}'.format(1/frequency[np.argmax(power)]))

best_frequency = frequency[np.argmax(power)]
t_fit = np.linspace(0.0, 1.0, 1000)
ls = LombScargle(t, y, dy)
y_fit = ls.model(t_fit, best_frequency)

phases = pyasl.foldAt(t, period=1/frequency[np.argmax(power)])

plt.errorbar(phases, y, yerr=dy, fmt='.b')
plt.plot(t_fit, y_fit)
plt.ylabel('Measurement values')
plt.xlabel('Time (d)')
plt.show()
$\endgroup$
2
  • $\begingroup$ Have you checked whether the frequency corresponding to 1 day is contained in the peridogram? Maybe autopower is failing to include that frequency. You could try to force it to include it by using the options "minimum_frequency" and "maximum_frequency" $\endgroup$
    – Prallax
    Sep 21 at 11:14
  • $\begingroup$ @Prallax, Yes, the periodogram does include the 1d frequency in the periodogram and finds it as its peak, i.e. period with max. power. The issue is that when I use ls.model() to calculate y_fit values for a period of 1.2 days with t_fit ranging from 0-1, it instead just returns the y_fit values actually ranging from 0 - 1.2 days and only returns those up to 1d. If I were to change the t_fit to (0.0, 1.2), it shows the complete phase of 1.2 days but its not fit to the data at all $\endgroup$ Sep 21 at 11:34
1
$\begingroup$
y_fit = ls.model(t_fit, best_frequency)

I figured it out. The input t_fit in the code line above is taken as 0.0 to 1.0 DAYS instead of PHASES so the y_fit values are the values of a 1.2d orbit but only up to 1.0 days. If the t_fit range is set instead as (0.0 - 1/best_frequency), we get the full phase of the period orbit. To then plot this full orbit ranging from 0 - 1 over the phase-folded data, the t_fit can be divided by 1/best_frequency and it then fits the data correctly.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.