{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3e048b8e",
   "metadata": {},
   "source": [
    "# Timing methods\n",
    "\n",
    "This notebook demonstrates some of the \"quick-look\" timing analyses available in `phoptic`, as well as how to extend the available timing analyses using [`stingray`](https://docs.stingray.science/en/stable/).\n",
    "\n",
    "## Generating and Reducing Data\n",
    "\n",
    "Before we can explore the different timing analyses provided in `phoptic`, we need some light curves to analyse. For this example, I'll generate and reduce some gappy observations (see the [reduction tutorial](reduction.ipynb) for more details on data reduction):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1dfb3b9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:25.792451Z",
     "iopub.status.busy": "2026-01-23T11:21:25.792360Z",
     "iopub.status.idle": "2026-01-23T11:22:24.596025Z",
     "shell.execute_reply": "2026-01-23T11:22:24.595706Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "os.environ['OMP_NUM_THREADS'] = '1'\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "import phoptic\n",
    "\n",
    "\n",
    "out_dir = Path('out')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3717730a",
   "metadata": {},
   "outputs": [],
   "source": [
    "phoptic.generate_gappy_observations(\n",
    "    out_directory=out_dir / 'data',  # path to the directory where the generated data will be saved\n",
    "    circular_aperture=False,  # disable circular aperture shadow\n",
    "    n_images=500,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "818383d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "reducer = phoptic.Reducer(\n",
    "    out_directory=out_dir / 'reduced',\n",
    "    data_directory=out_dir / 'data',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0cb355a",
   "metadata": {},
   "outputs": [],
   "source": [
    "reducer.create_catalogs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aac2caed",
   "metadata": {},
   "outputs": [],
   "source": [
    "phot = phoptic.OptimalPhotometer()\n",
    "reducer.photometry(phot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6259a29",
   "metadata": {},
   "outputs": [],
   "source": [
    "dphot = phoptic.DifferentialPhotometer(\n",
    "    out_directory=out_dir / 'reduced',\n",
    ")\n",
    "\n",
    "unnormed_analyser = dphot.get_relative_light_curve(\n",
    "    key='1:g',\n",
    "    target=2,\n",
    "    comparisons=[1, 3, 6],\n",
    "    phot_label='optimal',\n",
    "    prefix='example',\n",
    "    match_other_cameras=True,\n",
    "    show_diagnostics=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e341082",
   "metadata": {},
   "source": [
    "We now have an `Analyzer` instance that provides a number of quick-look timing analyses routines. Under-the-hood, these routines are mostly wrappers around `astropy` routines.\n",
    "\n",
    "By default, light curves are not normalised when performing differential photometry. To normalise our light curves to a mean of 1, we can define a new `Analyzer` instance with `norm='mean'`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e9f35ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "analyser = phoptic.Analyzer(\n",
    "    out_directory=out_dir / 'reduced',\n",
    "    light_curves=unnormed_analyser.light_curves,\n",
    "    norm='mean',\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79137bea",
   "metadata": {},
   "source": [
    "### `phoptic.Analyzer` routines\n",
    "\n",
    "Let's first take a look at our `Analyzer` light curves using the `plot()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbd3c0c6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:24.597191Z",
     "iopub.status.busy": "2026-01-23T11:22:24.597093Z",
     "iopub.status.idle": "2026-01-23T11:22:24.879315Z",
     "shell.execute_reply": "2026-01-23T11:22:24.879068Z"
    }
   },
   "outputs": [],
   "source": [
    "analyser.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a0641e5",
   "metadata": {},
   "source": [
    "As we can see, we have three variable, gappy light curves, all with a mean flux of 1. The filter for each light curve is shown in the top right. We can also see that this figure was saved to `timing_analysis_tutorial/reduced/plots/source_2_optimal_light_curves.pdf`.\n",
    "\n",
    "We can inspect our light curve data via the `light_curves` attribute. Within an `Analyzer`, light curves are stored as [`astropy.timeseries.TimeSeries`](https://docs.astropy.org/en/stable/timeseries/index.html) instances:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46bdc24a",
   "metadata": {},
   "source": [
    "We can also rebin our light curves to a lower time resolution using the `rebin()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65d175bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "\n",
    "binned_analyser = analyser.rebin(time_bin_size=10 * u.s)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "153ba679",
   "metadata": {},
   "source": [
    "This method takes a `time_bin_size` parameter, which must be a `astropy.units.Quantity`, and returns a new `Analyzer` instance containing the binned light curves. Under-the-hood, this uses `astropy`'s `aggregate_downsample()` function, which is slow for long timeseries.\n",
    "\n",
    "Let's take a look at the binned light curves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1c33c50",
   "metadata": {},
   "outputs": [],
   "source": [
    "binned_analyser.plot(save=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16a69531",
   "metadata": {},
   "source": [
    "`phoptic` handles the error propagation during binning and avoids padding gaps in the original light curves with zeros. In this case, the bin size was too large, and the variability of the light curves has been lost. As such, we'll carry on using the un-binned light curves.\n",
    "\n",
    "Let's now perform a few quick-look timing analyses on these light curves.\n",
    "\n",
    "### Lomb-Scargle Periodograms\n",
    "\n",
    "Often times, you might want to compute the Lomb-Scargle periodograms of your light curves to check for periododic signals. `astropy` implements two flavours of the Lomb-Scargle periodogram: [`LombScargle`](https://docs.astropy.org/en/stable/timeseries/lombscargle.html) and [`LombScargleMultiband`](https://docs.astropy.org/en/stable/timeseries/lombscarglemb.html). For convenience, `phoptic`'s `Analyzer` class implements wrappers for both routines that handle the boilerplate, allowing you to get straight to the analysis.\n",
    "\n",
    "#### Lomb-Scargle\n",
    "\n",
    "Let's first take a look at the standard Lomb-Scargle periodogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "347b8a5c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:24.901797Z",
     "iopub.status.busy": "2026-01-23T11:22:24.901714Z",
     "iopub.status.idle": "2026-01-23T11:22:25.099235Z",
     "shell.execute_reply": "2026-01-23T11:22:25.099023Z"
    }
   },
   "outputs": [],
   "source": [
    "lsps = analyser.lomb_scargle(scale='semilogx')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee343dce",
   "metadata": {},
   "source": [
    "We can see that each light curve has a dedicated periodogram, and the corresponding filter is shown in the top right of each panel. In this case, all three periodograms show a sharp peak between 0.1 and 0.2 Hz.\n",
    "\n",
    "By default, the `autofrequency()` method of `astropy`'s `LombScargle` class is used to construct the frequency grid of the periodogram. Alternatively, you may pass your own frequency grid. To avoid confusion with units, custom frequency grids must be passed as an `astropy` [`Quantity`](https://docs.astropy.org/en/stable/units/quantity.html). For example, let's recompute these periodograms in the 0.1-0.2 Hz frequency range. However, we do not want to overwrite our previous plot, so we'll also pass `save=False` to prevent the new plot from being saved:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50778d7d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:25.100095Z",
     "iopub.status.busy": "2026-01-23T11:22:25.100014Z",
     "iopub.status.idle": "2026-01-23T11:22:25.174779Z",
     "shell.execute_reply": "2026-01-23T11:22:25.174512Z"
    }
   },
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "import numpy as np\n",
    "\n",
    "frequency = np.linspace(0.1, 0.2, 100) * u.Hz  # units of Hz\n",
    "\n",
    "lsps = analyser.lomb_scargle(frequency=frequency, save=False, scale='semilogx')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bccc0f9d",
   "metadata": {},
   "source": [
    "Now we can see that all three bands suggest a signal with a frequency of between 0.12 and 0.14 Hz.\n",
    "\n",
    "Like many of the timing analysis methods of `Analyzer`, `lomb_scargle()` returns a dictionary of results that may be used to conduct further analyses. In this case, the keys list the filters and the values list the corresponding [`astropy.timeseries.LombScargle`](https://docs.astropy.org/en/stable/api/astropy.timeseries.LombScargle.html#astropy.timeseries.LombScargle.from_timeseries) instances: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2983f95f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:25.175636Z",
     "iopub.status.busy": "2026-01-23T11:22:25.175548Z",
     "iopub.status.idle": "2026-01-23T11:22:25.177851Z",
     "shell.execute_reply": "2026-01-23T11:22:25.177668Z"
    }
   },
   "outputs": [],
   "source": [
    "lsps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88c7f75b",
   "metadata": {},
   "source": [
    "Let's use the results returned by the `lomb_scargle()` method to precisely infer the peak frequencies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a9caa29",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:25.178623Z",
     "iopub.status.busy": "2026-01-23T11:22:25.178540Z",
     "iopub.status.idle": "2026-01-23T11:22:25.183284Z",
     "shell.execute_reply": "2026-01-23T11:22:25.183022Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "for fltr, lsp in lsps.items():\n",
    "    power = lsp.power(frequency=frequency)\n",
    "    peak_index = np.argmax(power)\n",
    "    print(f'{fltr} Peak frequency = {frequency[peak_index]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f8aeed0",
   "metadata": {},
   "source": [
    "Now we can see that all three bands show a periodicity at $\\sim$ 0.135 Hz which, as we can see at the top of this notebook, is the true frequency of variability for our source of interest.\n",
    "\n",
    "#### Multiband Lomb-Scargle\n",
    "\n",
    "In some cases, it may be preferable to compute a single Lomb-Scargle periodogram for all light curves; this can be done using `astropy`'s [`LombScargleMultiband`](https://docs.astropy.org/en/stable/timeseries/lombscarglemb.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afb32a02",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:25.184079Z",
     "iopub.status.busy": "2026-01-23T11:22:25.184005Z",
     "iopub.status.idle": "2026-01-23T11:22:25.646343Z",
     "shell.execute_reply": "2026-01-23T11:22:25.646128Z"
    }
   },
   "outputs": [],
   "source": [
    "periodograms = analyser.multiband_lomb_scargle(scale='semilogx')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f31895ce",
   "metadata": {},
   "source": [
    "As with `lomb_scargle()`, custom frequency grids can also be passed to `multiband_lomb_scargle()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "722658d9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:25.647186Z",
     "iopub.status.busy": "2026-01-23T11:22:25.647105Z",
     "iopub.status.idle": "2026-01-23T11:22:25.738817Z",
     "shell.execute_reply": "2026-01-23T11:22:25.738592Z"
    }
   },
   "outputs": [],
   "source": [
    "periodograms = analyser.multiband_lomb_scargle(frequency=frequency, save=False, scale='semilogx')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a19103e0",
   "metadata": {},
   "source": [
    "### Phase Folding/Binning\n",
    "\n",
    "Now that we have identified a candidate signal, we may want fold our light curves to see what it looks like. `astropy`'s `TimeSeries` class implements a `fold()` method, and `phoptic`'s `Analyzer` class provides a convenience wrapper for this method with some added functionality (binning). To fold a light curve, we need a period on which to fold; to avoid confusion with units, `phoptic` requires that the period being passed is a `Quantity`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "961e7042",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:25.739800Z",
     "iopub.status.busy": "2026-01-23T11:22:25.739716Z",
     "iopub.status.idle": "2026-01-23T11:22:26.144849Z",
     "shell.execute_reply": "2026-01-23T11:22:26.144624Z"
    }
   },
   "outputs": [],
   "source": [
    "f_signal = 0.135 * u.Hz  # periodogram peak frequency\n",
    "p_signal = (1 / f_signal).to(u.s)  # fold period in seconds\n",
    "\n",
    "folded_lcs = analyser.fold(\n",
    "    period=p_signal,\n",
    "    save=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d89a69f8",
   "metadata": {},
   "source": [
    "The `fold()` method of `phoptic.Analyzer` returns a `TimeSeries` whose `\"time\"` column quantifies the folded phase:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b9b90a9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.145792Z",
     "iopub.status.busy": "2026-01-23T11:22:26.145706Z",
     "iopub.status.idle": "2026-01-23T11:22:26.164094Z",
     "shell.execute_reply": "2026-01-23T11:22:26.163883Z"
    }
   },
   "outputs": [],
   "source": [
    "folded_lcs.pprint()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f046be8f",
   "metadata": {},
   "source": [
    "Additionally, the `fold()` method can be used to phase bin your light curves by passing the desired number of bins to to the `nbins` parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e2d91aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "folded_lcs = analyser.fold(\n",
    "    period=p_signal,\n",
    "    nbins=10,\n",
    "    save=False,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6426c335",
   "metadata": {},
   "source": [
    "Phase binning more clearly reveals the pulse shape compared to simply folding. From the above plots, it's clear that there is a lag in the pulsations between the different bands. Currently, `phoptic` does not support quantify lags natively. However, `phoptic`'s `Analyzer` class can be used to export your light curves to [`stingray`](https://docs.stingray.science/en/stable/), which provides a much more complete timing analysis framework - including ways to compute lags. Let's take a look at this now.\n",
    "\n",
    "## `stingray`\n",
    "\n",
    "`stingray` is a feature-rich spectral timing Python package developed primarily for use with X-ray data. `stingray` implements a number of generic timing analyis routines that can be applied to arbitrary time series. However, since `stingray` was developed with X-ray data in mind, it makes certain assumptions about the data. These assumptions can cause issues for non-X-ray data, as we'll see soon.\n",
    "\n",
    "To export your light curves to `stingray`, call the `export_light_curves_to_stingray()` method of `phoptic.Analyzer`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83afac1b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.164913Z",
     "iopub.status.busy": "2026-01-23T11:22:26.164835Z",
     "iopub.status.idle": "2026-01-23T11:22:26.282510Z",
     "shell.execute_reply": "2026-01-23T11:22:26.282252Z"
    }
   },
   "outputs": [],
   "source": [
    "stingray_lcs = analyser.export_light_curves_to_stingray()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c87f4afc",
   "metadata": {},
   "source": [
    "We can see that a warning about error distributions has been triggered. We'll take a look at this in more detail later. For now, let's ignore it.\n",
    "\n",
    "`export_light_curves_to_stingray()` returns a dictionary of {filter: [`stingray.Lightcurve`](https://docs.stingray.science/en/stable/notebooks/Lightcurve/Lightcurve%20tutorial.html)} pairs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bb5f377",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.283510Z",
     "iopub.status.busy": "2026-01-23T11:22:26.283423Z",
     "iopub.status.idle": "2026-01-23T11:22:26.285427Z",
     "shell.execute_reply": "2026-01-23T11:22:26.285245Z"
    }
   },
   "outputs": [],
   "source": [
    "stingray_lcs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87cd0e8b",
   "metadata": {},
   "source": [
    "Let's take a closer look at one of these light curves using the `plot()` method of `stingray.Lightcurve`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3497aa3e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.286187Z",
     "iopub.status.busy": "2026-01-23T11:22:26.286109Z",
     "iopub.status.idle": "2026-01-23T11:22:26.333142Z",
     "shell.execute_reply": "2026-01-23T11:22:26.332913Z"
    }
   },
   "outputs": [],
   "source": [
    "ax = stingray_lcs['1:g'].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42a6d13b",
   "metadata": {},
   "source": [
    "As we can see, the light curve times are converted into seconds, with $t=0$ corresponding to the `t_ref` attribute of the `Analyzer` instance. Additionally, `phoptic` has automatically inferred GTIs for each light curve, which are required for compatibility with many of `stingray`'s methods. The shaded regions in the above plot represent the time between GTIs. Let's use these `stingray.Lightcurve` instances to do some more advanced analyses.\n",
    "\n",
    "### Lags\n",
    "\n",
    "#### Cross-correlation\n",
    "\n",
    "Having recognised a lag between our light curves, there are several ways we can quantify it. For example, we can use the cross-correlations between light curves. This is a relatively simple approach, but it has the benefit of working with arbitrary time series. We can compute cross-correlations in `stingray` using [`CrossCorrelation`](https://docs.stingray.science/en/stable/notebooks/CrossCorrelation/cross_correlation_notebook.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2024f430",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.334009Z",
     "iopub.status.busy": "2026-01-23T11:22:26.333914Z",
     "iopub.status.idle": "2026-01-23T11:22:26.374464Z",
     "shell.execute_reply": "2026-01-23T11:22:26.374258Z"
    }
   },
   "outputs": [],
   "source": [
    "from stingray.crosscorrelation import CrossCorrelation\n",
    "\n",
    "cross_corr = CrossCorrelation(\n",
    "    lc1=stingray_lcs['1:g'],\n",
    "    lc2=stingray_lcs['3:i'],\n",
    "    )\n",
    "\n",
    "ax = cross_corr.plot(\n",
    "    labels=['Lag [s]', 'Correlation'],\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb2fadbd",
   "metadata": {},
   "source": [
    "We can see that the cross-correlation spectrum peaks at a non-zero lag. We can get the exact time lag from the `time_shift` attribute of the `CrossCorrelation` instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ff28660",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.375268Z",
     "iopub.status.busy": "2026-01-23T11:22:26.375183Z",
     "iopub.status.idle": "2026-01-23T11:22:26.376860Z",
     "shell.execute_reply": "2026-01-23T11:22:26.376683Z"
    }
   },
   "outputs": [],
   "source": [
    "print(f'Time lag: {cross_corr.time_shift:.2f}s')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bc03897",
   "metadata": {},
   "source": [
    "We find a negative lag of about 4 s, meaning $i$ is lagging $g$ (as expected).\n",
    "\n",
    "When we generated this data, we were told that the $i$-band variability has a phase lag of 3.142 (i.e., $1 \\pi$) radians, while the $g$-band variability has a phase lag of 0 radians. We can convert this time lag into a phase lag by dividing it by the period of the variability:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5b58506",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.377592Z",
     "iopub.status.busy": "2026-01-23T11:22:26.377516Z",
     "iopub.status.idle": "2026-01-23T11:22:26.379108Z",
     "shell.execute_reply": "2026-01-23T11:22:26.378943Z"
    }
   },
   "outputs": [],
   "source": [
    "print(f'Phase lag: {2 * np.pi * cross_corr.time_shift / p_signal.to_value(u.s):.2f} radians')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc4e239c",
   "metadata": {},
   "source": [
    "The phase lag is a little higher than expected, but this is due to the time resolution of our light curves.\n",
    "\n",
    "Cross-correlating generally works well, provided the light curves are dominated by the signal of interest. In many cases, however, light curves show variability on many time-scales, and we may want to isolate the lags of signals at specific frequencies. In such cases, it may be preferable to compute the cross-spectrum.\n",
    "\n",
    "#### Cross-spectra\n",
    "\n",
    "##### A Cautionary Tale\n",
    "\n",
    "To compute cross-spectra, we can use `stingray`'s `AveragedCrossspectrum` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc7360e5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.379816Z",
     "iopub.status.busy": "2026-01-23T11:22:26.379740Z",
     "iopub.status.idle": "2026-01-23T11:22:26.452904Z",
     "shell.execute_reply": "2026-01-23T11:22:26.452658Z"
    }
   },
   "outputs": [],
   "source": [
    "from stingray import AveragedCrossspectrum\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "segment_size = 30 * u.s\n",
    "\n",
    "avg_cs = AveragedCrossspectrum.from_lightcurve(\n",
    "    lc1=stingray_lcs['1:g'],\n",
    "    lc2=stingray_lcs['3:i'],\n",
    "    norm='frac',\n",
    "    segment_size=segment_size.to_value(u.s),  # stingray assumes all time values are in seconds\n",
    ")\n",
    "\n",
    "avg_cs_amplitude = np.abs(avg_cs.power)\n",
    "\n",
    "fig, ax = plt.subplots(tight_layout=True)\n",
    "\n",
    "ax.step(\n",
    "    avg_cs.freq,  # convert freq from cyc/d to Hz\n",
    "    avg_cs_amplitude,\n",
    "    where='mid',\n",
    "    c='k',\n",
    ")\n",
    "\n",
    "ax.set_xlabel('Frequency [Hz]')\n",
    "ax.set_ylabel('Power [(frac. rms)$^2$ / Hz]')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f1fc122",
   "metadata": {},
   "source": [
    "As we can see, the cross spectrum amplitude looks similar to the Lomb-Scargle periodograms we saw earlier, which is what we would expect. So far so good. Let's try to compute the coherence between these two light curves and check it's in the expected range ($0 \\leq \\text{coherence} \\leq 1$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fcabe0b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.453961Z",
     "iopub.status.busy": "2026-01-23T11:22:26.453817Z",
     "iopub.status.idle": "2026-01-23T11:22:26.456209Z",
     "shell.execute_reply": "2026-01-23T11:22:26.455957Z"
    }
   },
   "outputs": [],
   "source": [
    "coh, coh_err = avg_cs.coherence()\n",
    "\n",
    "print(np.logical_and(coh >= 0, coh <= 1).all())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d985924",
   "metadata": {},
   "source": [
    "Uh oh, we have unphysical coherence values! Let's take a look at some of them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22f3456c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.457147Z",
     "iopub.status.busy": "2026-01-23T11:22:26.456998Z",
     "iopub.status.idle": "2026-01-23T11:22:26.458977Z",
     "shell.execute_reply": "2026-01-23T11:22:26.458771Z"
    }
   },
   "outputs": [],
   "source": [
    "for i in range(5):\n",
    "    print(f'Coherence: {coh[i]}, error: {coh_err[i]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f9d4bb2",
   "metadata": {},
   "source": [
    "Our coherence estimates are in the thousands, and our error estimates in the *negative thousands*. Something is very wrong...\n",
    "\n",
    "The issue is due to `stingray` 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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22bc3840",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.459864Z",
     "iopub.status.busy": "2026-01-23T11:22:26.459741Z",
     "iopub.status.idle": "2026-01-23T11:22:26.540569Z",
     "shell.execute_reply": "2026-01-23T11:22:26.540325Z"
    }
   },
   "outputs": [],
   "source": [
    "plags, plags_err = avg_cs.phase_lag()\n",
    "\n",
    "fig, ax = plt.subplots(tight_layout=True)\n",
    "\n",
    "ax.errorbar(\n",
    "    avg_cs.freq,\n",
    "    plags,\n",
    "    plags_err,\n",
    "    fmt='kx',\n",
    ")\n",
    "\n",
    "ax.axvline(\n",
    "    0.135,\n",
    "    c='red',\n",
    "    label='Lag frequency'\n",
    "    )\n",
    "ax.axhline(\n",
    "    np.pi,\n",
    "    c='blue',\n",
    "    label='Expected lag'\n",
    "    )\n",
    "ax.legend()\n",
    "\n",
    "ax.set_xlabel('Frequency [Hz]')\n",
    "ax.set_ylabel('Phase lag [rad]')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "367c021e",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53c5cb83",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:22:26.541419Z",
     "iopub.status.busy": "2026-01-23T11:22:26.541333Z",
     "iopub.status.idle": "2026-01-23T11:22:26.543047Z",
     "shell.execute_reply": "2026-01-23T11:22:26.542873Z"
    }
   },
   "outputs": [],
   "source": [
    "print(plags_err)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4256290a",
   "metadata": {},
   "source": [
    "Our errors are all `NaN`s! This is another consequence of our light curves not being Poissonian. For cross spectral analyses, `stingray` may therefore not be suitable. For power spectra, however, `stingray` can be used without issue provided a suitable normalisation is specified (e.g., `norm=\"frac\"`).\n",
    "\n",
    "That concludes the timing methods tutorial for `phoptic`! The timing methods of `Analyzer` demonstrated here are intended as \"quick-look\" analyses, while more complete analyses are made easy thanks to `stingray`. Note, however, that `stingray` may not always work as intended when applied to relative light curves."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "opticam7",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}