Cycle determination (MSAP4-06)#

This notebook provides the test cases described in the MSAP4-06 submodule documentation.

import plato_msap4_demonstrator as msap4

import numpy as np
import matplotlib.pyplot as plt

Preprocessing#

In order to demonstrate the ability of MSAP4-06, we choose KIC~3733735, as its observing time/rotation period ratio is important and will allow us to perform meaningful computation on the \(S_\mathrm{ph}\) time series. We then start by computing this ACF time series, which is normally an IDP input provided by MSAP4-03. We consider the 2.57 days rotation period from Santos et al. (2021).

filename = msap4.get_target_filename (msap4.timeseries, '003733735')
t, s, dt = msap4.load_resource (filename)
pcutoff = (t[-1]-t[0])/2

prot = 2.57
_, t_sph, sph = msap4.compute_sph (t, s, prot,
                                          return_timeseries=True)

fig, ax = plt.subplots (1, 1)

ax.plot (t_sph, sph, color='darkorange', zorder=-1)
ax.scatter (t_sph, sph, color='darkorange', edgecolor='black',
            marker='o', s=40)

ax.set_xlabel ('Time (day)')
ax.set_ylabel (r'$S_\mathrm{ph}$ (ppm)')

fig.tight_layout ()

plt.savefig ('figures/kic3733735_sph_timeseries.png', dpi=300)
../../_images/cycle_determination_6_0.png

Computing the ACF and GLS of the ACF time series#

Now that we have our \(S_\mathrm{ph}\) time series, let’s compute its autocorrelation function and analyse it to extract periodicities above \(P_\mathrm{thresh}\).

dt_sph = np.median (np.diff (t_sph))
p_acf_sph, acf_sph = msap4.compute_acf (sph - np.mean (sph), dt_sph, normalise=True,
                                        use_scipy_correlate=True, smooth=False)
_, _, _, _, pmods_sph_acf, hacf, gacf = msap4.find_period_acf (p_acf_sph, acf_sph, pcutoff=pcutoff)
fig = msap4.plot_acf (p_acf_sph, acf_sph, prot=pmods_sph_acf,
                      xlim=(0,750), filename='figures/kic3733735_sph_acf.png')
../../_images/cycle_determination_9_0.png 
pmods_sph_acf, hacf, gacf

(array([ 89.82634128, 166.82034809, 269.47902383, 346.47303064]),
 array([0.54345038, 0.39418935, 0.23732771, 0.0511166 ]),
 array([0.24586367, 0.21436216, 0.22554245, 0.14562127]))

The second step is to compute the Lomb-Scargle periodogram of our \(S_\mathrm{ph}\) time series.

p_ps, ls, ps_object = msap4.compute_lomb_scargle_sph (t_sph, sph)
(pmods_sph_fourier, e_p,
 E_p, param, h_ps) = msap4.compute_prot_err_gaussian_fit_chi2_distribution (p_ps[p_ps<pcutoff], ls[p_ps<pcutoff],
                                                                            n_profile=5, threshold=0.1, verbose=False)
fig = msap4.plot_ls (p_ps, ls, filename='figures/kic3733735_sph_fourier.png',
                     logscale=False, param_profile=param,
                     ylim=(0, 0.1),
                     xlim=(2*dt_sph, 700))
../../_images/cycle_determination_12_0.png

Building the \(S_\mathrm{ph}\) intermediate data products#

We build here the intermediate data products related to the \(S_\mathrm{ph}\) analysis.

IDP_123_LONGTERM_MODULATION_SPH_FOURIER = msap4.prepare_idp_fourier (param, h_ps, ls.size,
                                                             pcutoff=pcutoff, pthresh=None,
                                                             fapcutoff=1)
IDP_123_LONGTERM_MODULATION_SPH_TIMESERIES = np.c_[pmods_sph_acf,
                                           np.full (pmods_sph_acf.size, -1),
                                           np.full (pmods_sph_acf.size, -1),
                                           hacf, gacf,
                                           np.arange (pmods_sph_acf.size)+1]

pd.DataFrame (data=IDP_123_LONGTERM_MODULATION_SPH_FOURIER)
0 1 2 3 4
0 90.525132 6.492736 7.580067 0.069541 0.651238
1 123.138563 14.257576 18.554155 0.028122 0.999998
2 146.368087 3.438629 3.608162 0.025239 1.000000
3 245.999501 10.181015 11.099773 0.022767 1.000000
4 184.089683 12.047067 13.861266 0.016895 1.000000 
pd.DataFrame (data=IDP_123_LONGTERM_MODULATION_SPH_TIMESERIES)
0 1 2 3 4 5
0 89.826341 -1.0 -1.0 0.543450 0.245864 1.0
1 166.820348 -1.0 -1.0 0.394189 0.214362 2.0
2 269.479024 -1.0 -1.0 0.237328 0.225542 3.0
3 346.473031 -1.0 -1.0 0.051117 0.145621 4.0

Comparing the long term modulations#

Finally, we complete our set with mock (and arbitrary) data to illustrate how long term modulations from different IDP should be compared.

IDP_123_LONGTERM_MODULATION_FOURIER = np.array ([[90, 3, 3, 1, 1e-16],
                                                  [130, 5, 5, 1, 1e-16]])
IDP_123_LONGTERM_MODULATION_TIMESERIES = np.array ([[91, -1, -1, .5, .2, 1],
                                                    [180, -1, -1, .5, .2, 2]])

DP4_123_LONGTERM_MODULATION = msap4.build_long_term_modulation (
                                IDP_123_LONGTERM_MODULATION_FOURIER,
                                IDP_123_LONGTERM_MODULATION_TIMESERIES,
                                IDP_123_LONGTERM_MODULATION_SPH_FOURIER,
                                IDP_123_LONGTERM_MODULATION_SPH_TIMESERIES
                                )

DP4_123_LONGTERM_MODULATION

array([[90.        ,  3.        ,  3.        , 91.        , -1.        ,
        -1.        , 90.52513177,  6.49273646,  7.58006719, 89.82634128,
        -1.        , -1.        ]])

pd.DataFrame (data=DP4_123_LONGTERM_MODULATION)
0 1 2 3 4 5 6 7 8 9 10 11
0 90.0 3.0 3.0 91.0 -1.0 -1.0 90.525132 6.492736 7.580067 89.826341 -1.0 -1.0