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()
Improve this question
$\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, 2021 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$
    – Mukul Kumar
    Sep 21, 2021 at 11:34

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .