Timing methods

This notebook demonstrates some of the “quick-look” timing analyses available in opticam.

Preface

stingray is designed primarily for analysing X-ray data, and so is perhaps not ideal for OPTICAM data products. That said, stingray provides a number of useful routines that can (mostly) be applied to arbitrary time series, provided some care is taken.

Generating and Reducing Data

Before we can explore the different timing analyses provided in opticam, we need some light curves to analyse. For this example, I’ll quickly generate and reduce some gappy observations (see the Reduction Tutorial for more details on data reduction):

import opticam

opticam.generate_gappy_observations(
    out_dir='timing_analysis_tutorial/data',
    n_images=500,  # generate lots of observations to create long light curves
    binning_scale=8,  # use large binning scale to reduce disk space
    circular_aperture=False,  # disable circular aperture shadow for simplicity
)

reducer = opticam.Reducer(
    data_directory='timing_analysis_tutorial/data',
    out_directory='timing_analysis_tutorial/reduced',
    show_plots=True,
    verbose=True,
    remove_cosmic_rays=False,
    )
reducer.create_catalogs()

phot = opticam.OptimalPhotometer()
reducer.photometry(phot)

[OPTICAM] variable source is at (65, 57)
[OPTICAM] variability frequency: 0.135 Hz
[OPTICAM] variability phase lags:
    [OPTICAM] g-band: 0.000 radians
    [OPTICAM] r-band: 1.571 radians
    [OPTICAM] i-band: 3.142 radians

Generating observations: 100%|██████████|[00:21<00:00]

[OPTICAM] timing_analysis_tutorial/reduced not found, attempting to create ...
[OPTICAM] timing_analysis_tutorial/reduced created.

[OPTICAM] Scanning data directory: 100%|██████████|[00:01<00:00]
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101g200000108o.fits.gz and 240101g200000374o.fits.gz (48.995 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101g200000449o.fits.gz and 240101g200000010o.fits.gz (19.998 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101g200000469o.fits.gz and 240101g200000233o.fits.gz (58.994 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101g200000369o.fits.gz and 240101g200000415o.fits.gz (10.999 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101g200000319o.fits.gz and 240101g200000464o.fits.gz (24.998 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101i200000118o.fits.gz and 240101i200000169o.fits.gz (48.995 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101i200000045o.fits.gz and 240101i200000378o.fits.gz (19.998 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101i200000347o.fits.gz and 240101i200000102o.fits.gz (58.994 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101i200000004o.fits.gz and 240101i200000167o.fits.gz (10.999 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101i200000125o.fits.gz and 240101i200000036o.fits.gz (24.998 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101r200000375o.fits.gz and 240101r200000344o.fits.gz (48.995 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101r200000321o.fits.gz and 240101r200000048o.fits.gz (19.998 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101r200000489o.fits.gz and 240101r200000446o.fits.gz (58.994 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101r200000149o.fits.gz and 240101r200000456o.fits.gz (10.999 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/opticam/utils/data_checks.py:205: UserWarning: [OPTICAM] Large time gap detected between 240101r200000166o.fits.gz and 240101r200000194o.fits.gz (24.998 s compared to the median time difference of 1.000 s). This may cause alignment issues. If so, consider moving all files after this gap to a separate directory.
  warnings.warn(string)

[OPTICAM] Binning: 8x8
[OPTICAM] Filters: g-band, r-band, i-band
[OPTICAM] 341 g-band images.
[OPTICAM] 341 r-band images.
[OPTICAM] 341 i-band images.

[OPTICAM] Creating source catalogs

[OPTICAM] Aligning g-band images: 100%|██████████|[00:02<00:00]

[OPTICAM] Done.
[OPTICAM] 341 image(s) aligned.
[OPTICAM] 0 image(s) could not be aligned.

[OPTICAM] Aligning r-band images: 100%|██████████|[00:02<00:00]

[OPTICAM] Done.
[OPTICAM] 341 image(s) aligned.
[OPTICAM] 0 image(s) could not be aligned.

[OPTICAM] Aligning i-band images: 100%|██████████|[00:02<00:00]

[OPTICAM] Done.
[OPTICAM] 341 image(s) aligned.
[OPTICAM] 0 image(s) could not be aligned.

[OPTICAM] Photometry results will be saved to lcs/optimal in timing_analysis_tutorial/reduced.

[OPTICAM] Performing photometry on g-band images: 100%|██████████|[00:02<00:00]
[OPTICAM] Performing photometry on r-band images: 100%|██████████|[00:02<00:00]
[OPTICAM] Performing photometry on i-band images: 100%|██████████|[00:02<00:00]

dphot = opticam.DifferentialPhotometer(
    out_directory='timing_analysis_tutorial/reduced',
    show_plots=False,
)

analyser = dphot.get_relative_light_curve(
    fltr='g-band',
    target=6,
    comparisons=5,
    phot_label='optimal',
    match_other_cameras=True,
    show_diagnostics=False,
)

[OPTICAM] Filters: g-band, r-band, i-band
[OPTICAM] g-band target ID 6 was matched to r-band target ID 6
[OPTICAM] g-band comparison ID 5 was matched to r-band comparison ID 5
[OPTICAM] g-band target ID 6 was matched to i-band target ID 6
[OPTICAM] g-band comparison ID 5 was matched to i-band comparison ID 5

/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/stingray/utils.py:486: UserWarning: SIMON says: Stingray only uses poisson err_dist at the moment. All analysis in the light curve will assume Poisson errors. Sorry for the inconvenience.
  warnings.warn("SIMON says: {0}".format(message), **kwargs)

We now have an Analyzer instance that provides a number of quick-look timing analyses routines. Under-the-hood, these routines use the stingray Python package, so if a particular stingray routine is not implemented in opticam, a light curve from the Analyzer instance can be passed directly to stingray instead (more on this later). By using stingray, the results returned by opticam.Analyzer methods are intended to be compatible with existing analysis pipelines, or provide a feature-rich framework in which to develop analysis pipelines.

Let’s first take a look at our Analyzer light curves. Since these observations contain gaps, we can also visualise the Good Time Intervals by passing show_gtis=True to plot_light_curves():

fig = analyser.plot_light_curves(
    show_gtis=True,
)

As we can see, we have three clearly variable, gappy light curves. Let’s now perform a few quick-look timing analyses on these light curves.

Periodograms

Often times, users might want to compute the periodograms of their light curves. stingray implements both the Lomb-Scargle periodogram (LSP) and the classical periodogram (FFT). Both of these are also implemented in opticam for convenience.

Since opticam computes relative light curves, while stingray assumes that light curves are in units of counts, any Analyzer methods that use stingray routines will default to norm="frac". The fractional RMS normalization is valid for relative light curves, while the absolute RMS normalization (norm="abs"), will give the same results as the fractional RMS normalisation provided the light curve is normalised to a mean flux of 1. In general, try to stick with norm="frac" since it’s scale-invariant (more on this later).

Let’s first take a look at the Lomb-Scargle periodogram:

lsps = analyser.compute_lomb_scargle_periodograms()

We have three Lomb-Scarlge periodograms, all of which show a strong peak at a frequency of about 0.14 Hz.

Like many of the timing analysis methods of Analyzer, compute_lomb_scargle_periodograms() returns a dictionary of results. In this case, the keys list the filters and the values list the corresponding stingray.lombscargle.LombScarglePowerspectrum instances:

lsps

{'g-band': <stingray.lombscargle.LombScarglePowerspectrum at 0x73179933c6e0>,
 'r-band': <stingray.lombscargle.LombScarglePowerspectrum at 0x7317993c7b10>,
 'i-band': <stingray.lombscargle.LombScarglePowerspectrum at 0x7317993c7ed0>}

Let’s find the peak frequency for each periodogram:

import numpy as np

for fltr, lsp in lsps.items():
    print('Peak frequency:', lsp.freq[np.argmax(lsp.power)], 'Hz')

Peak frequency: 0.13528395820116063 Hz
Peak frequency: 0.13528395820116063 Hz
Peak frequency: 0.13528395820116063 Hz

As we can see, all three bands show a periodicity at \(\sim\) 0.135 Hz (which is the true frequency of variability for our source of interest). We also see a bit of aliasing due to the gaps in our observations.

Let’s compare the LSP to the classical periodogram:

periodograms = analyser.compute_power_spectra()

[OPTICAM] Unable to compute periodogram for g-band light curve due to gaps. Consider using either the compute_lomb_scargle_periodograms() or compute_averaged_power_spectra() methods instead.
[OPTICAM] Unable to compute periodogram for r-band light curve due to gaps. Consider using either the compute_lomb_scargle_periodograms() or compute_averaged_power_spectra() methods instead.
[OPTICAM] Unable to compute periodogram for i-band light curve due to gaps. Consider using either the compute_lomb_scargle_periodograms() or compute_averaged_power_spectra() methods instead.

As we can see, computing the classical periodogram fails because the data contain gaps. Let’s try computing the averaged periodogram instead. When computing an averaged periodogram, a segment size must be specified; to avoid issues with units, opticam requires this segment size by an astropy.units.quantity.Quantity instance:

import astropy.units as u

averaged_periodograms = analyser.compute_averaged_power_spectra(
    segment_size=50 * u.s,
)

[OPTICAM] 5 g-band segments averaged.
[OPTICAM] 5 r-band segments averaged.
[OPTICAM] 5 i-band segments averaged.

Computing averaged periodograms works, since the gaps can be filtered out, but we don’t really have enough data in this case to get meaningful results (note each averaged periodogram is comprised of just five segments). When averaging Fourier products, the general consensus is to average at least 30 segments to get trustworthy results. However, generating enough data to get 30 complete 50 s segments would take up a lot of disk space, so I’ll leave this as an exercise for the user.

Phase Folding

It is common to phase fold light curves on candidate frequencies or periods to visualise the pulse profile. This is implemented in opticam for convenience via the phase_fold_light_curves() method of Analyzer. First, however, we must define our candidate period. To avoid issues with units, candidate periods must be an astropy.units.quantity.Quantity instance:

f_signal = 0.135 * u.Hz  # periodogram peak frequency
p_signal = 1 / f_signal

phase_folds = analyser.phase_fold_light_curves(
    period=p_signal,
)

As we can see, each band shows a sinusoidal periodicity at 0.135 Hz. We can also phase bin to make the pulse profile clearer:

phase_bins = analyser.phase_bin_light_curves(
    period=p_signal,
    n_bins=20,
)

Now we can see the pulse profiles more clearly. It is also clear from the above plots that there is a lag in the pulsations between the different bands. Let’s quantify this.

Lags

One of the great things about OPTICAM is that it provides simultaneous three-colour observations. As such, it is possible to search for wavelength-dependent lags; opticam’s Analyzer class makes this easy via the compute_cross_correlations() method:

corrs = analyser.compute_cross_correlations()

The compute_cross_correlations() method computes the cross-correlations between each pair of light curves in the Analyzer instance. As we can see from the above plot, the \(r\)-band lags the \(g\)-band by 2 s (since the time shift is negative, the second filter lags the first filter), and the \(i\)-band lags the \(r\)-band by 2 s as well. The true lag between these bands is a quarter period; since the period is \(1 / 0.135 = 7.4\) s, this is 1.85 s. Within the time resolution of our light curves (1 s), we can therefore correctly recover the lag.

Using stingray Directly

A Cautionary Tale

Since Analyzer stores light curves as stingray.Lightcurve instances, we can pass them directly to a number of stingray timing analysis routines. However, stingray assumes light curves are in units of counts with Poisson errors, while our light curves are relative (unitless) with Gaussian errors. As such, care must be taken in regards to units (e.g., in normalisations), and some routines either won’t work as expected or won’t work at all. As an instructive example, let’s say we want to compute the cross spectrum between our \(g\)- and \(r\)-band light curves:

from stingray import AveragedCrossspectrum
from matplotlib import pyplot as plt

avg_cs = AveragedCrossspectrum.from_lightcurve(
    lc1=analyser.light_curves['g-band'],
    lc2=analyser.light_curves['r-band'],
    norm='frac',
    segment_size=(50 * u.s).to_value(u.d),  # stingray assumes all time values are in seconds but OPTICAM uses Barycentric MJD
)

avg_cs_amplitude = np.abs(avg_cs.power)

fig, ax = plt.subplots(tight_layout=True)

ax.step(
    avg_cs.freq / 86400,  # convert freq from cyc/d to Hz
    avg_cs_amplitude,
    where='mid',
    c='k',
)

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Power [(frac. rms)$^2$ / Hz]')

plt.show()

5it [00:00, 5283.83it/s]
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/stingray/fourier.py:1134: UserWarning: n_ave is below 30. Please note that the error bars on the quantities derived from the cross spectrum are only reliable for a large number of averaged powers.
  warnings.warn(
/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/stingray/fourier.py:1161: RuntimeWarning: invalid value encountered in sqrt
  dphi = np.sqrt((1 - gsq) / (2 * gsq * n_ave))

That looks how we’d expect. So far so good. Let’s try to compute the coherence between these two light curves:

fig, ax = plt.subplots(tight_layout=True)

coh, coh_err = avg_cs.coherence()

ax.errorbar(
    avg_cs.freq / 86400,
    coh,
    coh_err,
    c='black',
)

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Coherence')

plt.show()

/home/zac/miniforge3/envs/opticam/lib/python3.13/site-packages/stingray/utils.py:486: UserWarning: SIMON says: Number of segments used in averaging is significantly low. The result might not follow the expected statistical distributions.
  warnings.warn("SIMON says: {0}".format(message), **kwargs)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[13], line 5
      1 fig, ax = plt.subplots(tight_layout=True)
      3 coh, coh_err = avg_cs.coherence()
----> 5 ax.errorbar(
      6     avg_cs.freq / 86400,
      7     coh,
      8     coh_err,
      9     c='black',
     10 )
     12 ax.set_xlabel('Frequency [Hz]')
     13 ax.set_ylabel('Coherence')

File ~/miniforge3/envs/opticam/lib/python3.13/site-packages/matplotlib/__init__.py:1473, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1470 @functools.wraps(func)
   1471 def inner(ax, *args, data=None, **kwargs):
   1472     if data is None:
-> 1473         return func(
   1474             ax,
   1475             *map(sanitize_sequence, args),
   1476             **{k: sanitize_sequence(v) for k, v in kwargs.items()})
   1478     bound = new_sig.bind(ax, *args, **kwargs)
   1479     auto_label = (bound.arguments.get(label_namer)
   1480                   or bound.kwargs.get(label_namer))

File ~/miniforge3/envs/opticam/lib/python3.13/site-packages/matplotlib/axes/_axes.py:3743, in Axes.errorbar(self, x, y, yerr, xerr, fmt, ecolor, elinewidth, capsize, barsabove, lolims, uplims, xlolims, xuplims, errorevery, capthick, **kwargs)
   3740 res = np.zeros(err.shape, dtype=bool)  # Default in case of nan
   3741 if np.any(np.less(err, -err, out=res, where=(err == err))):
   3742     # like err<0, but also works for timedelta and nan.
-> 3743     raise ValueError(
   3744         f"'{dep_axis}err' must not contain negative values")
   3745 # This is like
   3746 #     elow, ehigh = np.broadcast_to(...)
   3747 #     return dep - elow * ~lolims, dep + ehigh * ~uplims
   3748 # except that broadcast_to would strip units.
   3749 low, high = dep + np.vstack([-(1 - lolims), 1 - uplims]) * err

ValueError: 'yerr' must not contain negative values

We get negative coherence errors, which is obviously unphysical, and we can’t even plot the it. This is due to how stingray is computing the cross spectrum (assuming Poisson statistics). Since we know the phase lag between these two light curves, let’s see what happens if we try to compute it from the averaged cross spectrum:

plags, plags_err = avg_cs.phase_lag()

fig, ax = plt.subplots(tight_layout=True)

ax.errorbar(
    avg_cs.freq / 86400,
    plags,
    plags_err,
    fmt='kx',
)

ax.axvline(
    0.135,
    c='red',
    label='Lag frequency'
    )
ax.axhline(
    np.pi / 2,
    c='blue',
    label='Expected lag'
    )
ax.legend()

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Phase lag [$\\pi$]')

plt.show()

We get the correct phase lag at the correct frequency. However, we don’t seem to have any error bars, and so it’s impossible to quantify the significance of this lag. Let’s see what our error bars are:

print(plags_err)

[nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan
 nan nan nan nan nan nan]

Our errors are all NaNs! This is another consequency of our light curves not having Poisson errors.

A Hacky Fix

To work around stingray’s requirement for Poisson errors, we can trick it into thinking our relative light curves are actually in units of counts by scaling them by a suitably large value (e.g., 10,000). We’ll also want to convert the time values from Barycentric MJD to seconds:

from stingray import Lightcurve

lc = analyser.light_curves['g-band']
t_ref = np.min(lc.time)
t = 86400 * (lc.time - t_ref)  # convert time from MJD to s
gti = 86400 * (lc.gti - t_ref)  # convert GTIs from MJD to s
f = lc.counts * 10000  # scale fluxes so stingray thinks they're counts
g_lc = Lightcurve(t, f, gti=gti)  # do not pass errors or define err_dist

lc = analyser.light_curves['r-band']
t = 86400 * (lc.time - t_ref)  # convert time from MJD to s
gti = 86400 * (lc.gti - t_ref)  # convert GTIs from MJD to s
f = lc.counts * 10000  # scale fluxes so stingray thinks they're counts
r_lc = Lightcurve(t, f, gti=gti)  # do not pass errors or define err_dist

avg_cs = AveragedCrossspectrum.from_lightcurve(
    g_lc,
    r_lc,
    segment_size=50,  # segment size in seconds
    norm='frac',
    )

5it [00:00, 6434.96it/s]

The fractional RMS normalisation (norm='frac') is scale-invariant, so it should give us the same averaged periodogram as before (provided we use the same segment length). Let’s verify this:

from stingray import AveragedPowerspectrum

avg_ps = AveragedPowerspectrum.from_lightcurve(
    g_lc,
    segment_size=50,
    norm='frac',
    )

print(np.allclose(avg_ps.power, averaged_periodograms['g-band'].power))

5it [00:00, 14453.15it/s]

True

As we can see, our scaled \(g\)-band light curve gives the same averaged power spectrum as our initial relative light curve. Let’s now try to compute the averaged cross spectrum, coherence, and phase lags using our scaled light curves:

fig, ax = plt.subplots(tight_layout=True)

ax.step(
    avg_cs.freq,
    avg_cs_amplitude,
    where='mid',
    c='k',
)

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Power [(frac. rms)$^2$ / Hz]')

plt.show()

As with the averaged power spectrum, the averaged cross spectrum looks the same as before (since we’re using norm="frac"). Again, so far so good. Let’s see if we can compute the coherence this time:

fig, ax = plt.subplots(tight_layout=True)

coh, coh_err = avg_cs.coherence()

ax.errorbar(
    avg_cs.freq,
    coh,
    coh_err,
    c='black',
)

ax.axvline(
    0.135,
    c='red',
)

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Coherence')

plt.show()

We have a coherence plot! We can see that the coherence is generally quite noisy, as we would expect for noisy light curves. At the frequency of variability, however, the coherence peaks. This is expected because we know both light curves show sinusoidal variability at this frequency.

Let’s see how the phase lags look now:

plags, plags_err = avg_cs.phase_lag()

fig, ax = plt.subplots(tight_layout=True)

ax.errorbar(
    avg_cs.freq,
    plags,
    plags_err,
    fmt='kx',
)

ax.axvline(
    0.135,
    c='red',
    label='Lag frequency'
    )
ax.axhline(
    np.pi / 2,
    c='blue',
    label='Expected lag'
    )
ax.axhline(
    0,
    c='black',
    ls='--'
    )
ax.legend()

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Phase lag')

plt.show()

We have errors on our phase lags! Let’s quantify the significance of the lag:

fig, ax = plt.subplots(tight_layout=True)

ax.step(
    avg_cs.freq,
    plags / plags_err,
    where='mid',
    c='black',
)

ax.fill_between(
    ax.set_xlim(),
    -5, 5,
    color='grey',
    alpha=.5,
    edgecolor='none',
)

ax.axvline(
    0.135,
    c='red',
    label='Lag frequency'
    )

ax.legend()

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Phase lag significance')

plt.show()

Now we can see that the phase lag is highly significant (\(> 30 \sigma\))!

That concludes the timing methods tutorial for opticam! The timing methods of Analyzer demonstrated here are intended as “quick-look” analyses, while more complete analyses are made easy thanks to stingray. Note that special care (and hacky workarounds) may be needed to get the most out of stingray.