{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Generating signals\n",
    "\n",
    "Example on how to generate continuous gravitational-wave signals."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "0. [Introduction](#introduction)\n",
    "0. [Generating an ensemble of signals](#ensemble)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import pyfstat\n",
    "from pyfstat.utils import get_sft_as_arrays\n",
    "\n",
    "# Local module to simplify plotting\n",
    "import tutorial_utils\n",
    "\n",
    "logger = pyfstat.set_up_logger(label=\"1_generating_signals\", log_level=\"INFO\")\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "## Introduction<a name=\"the-sft-format\"></a>\n",
    "Continuous gravitational-wave signals (CWs) are long-lasting forms of gravitational radiation.\n",
    "They are usually characterized using two sets of parameters, namely the *amplitude* parameters\n",
    "and the *Doppler* parameters, respectively referred to as $\\mathcal{A}$ and $\\lambda$.\n",
    "\n",
    "For the typical case of a rapidly-spinning neutron star, the amplitude parameter contain the nominal \n",
    "CW amplitude $h_0$, the (cosine of the) source inclination angle with respect to the line of sight $\\cos\\iota$,\n",
    "the polarization angle $\\psi$ and the initial phase of the wave $\\psi$.\n",
    "Depending on the emission mechanism, $h_0$ can be further described using further physical quantities such as\n",
    "the source's frequency, ellipticity, or distance to the detector.\n",
    "\n",
    "Doppler parameters describe the evolution of the gravitational-wave frequency both due to physical processes\n",
    "undergoing at the source and the motion of the interferometric detector with respect to the \n",
    "Solar system barycenter (i.e. the Sun). In this tutorial, we will limit ourselves to gravitational wave\n",
    "frequency $f_0$, spindown $f_1$ and sky position $\\hat{n}$, which we will parametrize using the right ascension\n",
    "$\\alpha$ and declination $\\delta$ angles.\n",
    "\n",
    "A detailed explanation of these parameters can be found in\n",
    "[this technical document](https://dcc.ligo.org/LIGO-T0900149/public).\n",
    "\n",
    "`pyfstat.Writer` allows these variables to be given as inputs to produce SFTs containing a CW signal and\n",
    "(optionally) Gaussian noise. Signal parameters are both available as attributes in the `pyfstat.Writer` instance\n",
    "as as a `.cff` file in `outdir`. \n",
    "Alternatively, as exemplified in [PyFstat_example_injecting_into_noise_sfts.py](../other_examples),\n",
    "the `noiseSFTs` can be used to provide a pre-generated set of SFTs as background noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "writer_kwargs = {\n",
    "    \"label\": \"SingleDetectorGaussianNoise\",\n",
    "    \"outdir\": \"PyFstat_example_data\",\n",
    "    \"tstart\": 1238166018,\n",
    "    \"duration\": 365 * 86400,\n",
    "    \"detectors\": \"H1\",\n",
    "    \"sqrtSX\": 1e-23,\n",
    "    \"Tsft\": 1800,\n",
    "    \"SFTWindowType\": \"tukey\",\n",
    "    \"SFTWindowParam\": 0.01,\n",
    "}\n",
    "\n",
    "signal_parameters = {\n",
    "    \"F0\": 100.0,\n",
    "    \"F1\": -1e-9,\n",
    "    \"Alpha\": 0.0,\n",
    "    \"Delta\": 0.0,\n",
    "    \"h0\": 1e-22,\n",
    "    \"cosi\": 1,\n",
    "    \"psi\": 0.0,\n",
    "    \"phi\": 0.0,\n",
    "    \"tref\": writer_kwargs[\"tstart\"],\n",
    "}\n",
    "\n",
    "writer = pyfstat.Writer(**writer_kwargs, **signal_parameters)\n",
    "writer.make_data()\n",
    "frequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "tutorial_utils.plot_real_imag_spectrograms(\n",
    "    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n",
    ")\n",
    "fig, ax = tutorial_utils.plot_amplitude_phase_spectrograms(\n",
    "    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "These plots show the typical features of CW signals.\n",
    "\n",
    "Short-scale amplitude modulations are due to the rotation of the Earth, \n",
    "as the detector has different sensitivities depending on its orientation.\n",
    "\n",
    "The frequency is modulated in a quasi-periodic fashion due to the translation of the interferometric\n",
    "detector around the Sun. The downwards trend is due to the presence of a negative source spindown due to the\n",
    "emission of energy as gravitational waves.\n",
    "\n",
    "The instantaneous frequency of a CW at the detector can be expressed as\n",
    "$$\n",
    "f(t) = \\left[f_0 + f_1 \\left(t - t_{\\text{ref}}\\right)\\right]\\left[1 + \\frac{\\vec{v} \\cdot \\hat{n}}{c}\\right],\n",
    "$$\n",
    "where $\\vec{v}/c$ is the velocity of the detector normalized by the speed of light in vacuum \n",
    "and $t_{\\text{ref}}$ is a fiducial reference time at which $f_0$ and $f_1$ are specified.\n",
    "\n",
    "Detector velocities can be retrieved from SFT files using `pyfstat.DetectorStates`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "states = pyfstat.DetectorStates().get_multi_detector_states_from_sfts(\n",
    "    writer.sftfilepath, central_frequency=writer.F0, time_offset=0\n",
    ")\n",
    "\n",
    "ts = np.array([data.tGPS.gpsSeconds for data in states.data[0].data])\n",
    "velocities = np.vstack([data.vDetector for data in states.data[0].data]).T\n",
    "\n",
    "n = np.array(\n",
    "    [\n",
    "        [\n",
    "            np.cos(writer.Alpha) * np.cos(writer.Delta),  # Cartesian X\n",
    "            np.sin(writer.Alpha) * np.cos(writer.Delta),  # Cartesian Y\n",
    "            np.sin(writer.Delta),  # Cartesian Z\n",
    "        ]\n",
    "    ]\n",
    ")\n",
    "\n",
    "f_inst = ((writer.F0 + (ts - writer.tref) * writer.F1) * (1 + np.dot(n, velocities)))[0]\n",
    "\n",
    "for a in ax:\n",
    "    for off in [-2, 2]:\n",
    "        a.plot(\n",
    "            (ts - ts[0]) / 1800,\n",
    "            f_inst + off / writer.Tsft,\n",
    "            color=\"white\",\n",
    "            ls=\"--\",\n",
    "            label=\"Instantaneous frequency\" if off == 2 else \"\",\n",
    "        )\n",
    "    a.legend()\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "The presence of visible CW features in the spectrogram can be quantified using the so-called\n",
    "[sensitivity depth](https://arxiv.org/abs/1808.02459), which is defined as the ratio between the noise\n",
    "and CW amplitude\n",
    "$$\n",
    "\\mathcal{D} = \\frac{\\sqrt{S_{\\text{n}}}}{h_0}\\;.\n",
    "$$\n",
    "\n",
    "Current searches for CW signals from unknown sources are able to detect a sensitivity depth ranging between\n",
    "$10\\;\\text{Hz}^{-1/2}$ and $50\\;\\text{Hz}^{-1/2}$. Visual features, however, tend to disappear around a depth\n",
    "of $20\\;\\text{Hz}^{-1/2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CW signal at a depth of 20\n",
    "signal_parameters[\"h0\"] = writer_kwargs[\"sqrtSX\"] / 20.0\n",
    "\n",
    "writer = pyfstat.Writer(**writer_kwargs, **signal_parameters)\n",
    "\n",
    "# Create SFTs\n",
    "writer.make_data()\n",
    "\n",
    "frequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "tutorial_utils.plot_real_imag_spectrograms(\n",
    "    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n",
    ")\n",
    "tutorial_utils.plot_amplitude_phase_spectrograms(\n",
    "    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "Alternatively, CW signals can be characterized in terms of their (squared) signal-to-noise ratio (SNR) $\\rho^2$.\n",
    "As opposed to sensitivity depth, which focuses on instantaneous features of the CW, $\\rho^2$ is an integrated\n",
    "quantity along the full duration of the observing run. \n",
    "The dependency of the *optimal* SNR (i.e. assuming the Doppler parameters match those of the signal) depends in\n",
    "a non-trivial way on the amplitude parameters and sky position of the source (due to the anisotropy of the\n",
    "detector response function) [[Eq. (77)](https://dcc.ligo.org/LIGO-T0900149/public)]; \n",
    "the average optimal SNR for a uniform distribution of sources across the sky with an isotropic polarization angle, however, \n",
    "can be readily expressed as\n",
    "$$\n",
    "\\langle \\rho^2 \\rangle_{\\vec{n}, \\psi} = \n",
    "\\frac{1}{20}\n",
    "\\frac{h_0^2 \\left(\\cos^4\\iota + 6 \\cos^2\\iota + 1\\right)}\n",
    "{S_{\\textrm{n}}}\n",
    "T_{\\textrm{obs}},\n",
    "$$\n",
    "where $T_{\\textrm{obs}}$ represents full duration of the data stream.\n",
    "\n",
    "The optimal SNR for a specific template can be computed using `pyfstat.SignalToNoiseRatio` assuming Gaussian noise with a given\n",
    "`sqrtSX` value or from a specific set of SFTs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "snr = pyfstat.SignalToNoiseRatio.from_sfts(F0=writer.F0, sftfilepath=writer.sftfilepath)\n",
    "squared_snr = snr.compute_snr2(\n",
    "    Alpha=writer.Alpha,\n",
    "    Delta=writer.Delta,\n",
    "    psi=writer.psi,\n",
    "    phi=writer.phi,\n",
    "    h0=writer.h0,\n",
    "    cosi=writer.cosi,\n",
    ")\n",
    "print(f\"SNR: {np.sqrt(squared_snr)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "## Generating an ensemble of signals<a name=\"ensemble\"></a>\n",
    "\n",
    "Characterizing a CW search or follow-up method usually involves analyzing its response to an ensemble of signals\n",
    "from a population of interest.\n",
    "For the case of all-sky searches, this population is usually consistent with a uniform distribution of sources\n",
    "across the sky with isotropic orientation.\n",
    "\n",
    "The `pyfstat.InjectionParametersGenerator` provides methods to draw parameters \n",
    "from generic distributions in a suitable format, ready to be fed into `pyfstat.Writer`.\n",
    "Uniform sampling across the sky is baked in the child class `pyfstat.AllSkyInjectionParametersGenerator`,\n",
    "and isotropic priors on amplitude parameters are available in \n",
    "`pyfstat.isotropic_amplitude_distribution`.\n",
    "\n",
    "All-sky searches tend to perform injection campaigns at fixed values of `h0` (or, equivalently, $\\mathcal{D}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate 10 signals with parameters drawn from a specific population\n",
    "num_signals = 10\n",
    "\n",
    "writer_kwargs = {\n",
    "    \"tstart\": 1238166018,\n",
    "    \"duration\": 365 * 86400,\n",
    "    \"detectors\": \"H1\",\n",
    "    \"sqrtSX\": 1e-23,\n",
    "    \"Tsft\": 1800,\n",
    "    \"SFTWindowType\": \"tukey\",\n",
    "    \"SFTWindowParam\": 0.01,\n",
    "}\n",
    "\n",
    "signal_parameters_generator = pyfstat.AllSkyInjectionParametersGenerator(\n",
    "    priors={\n",
    "        \"F0\": {\"stats.uniform\": {\"loc\": 100.0, \"scale\": 0.1}},\n",
    "        \"F1\": -1e-10,\n",
    "        \"F2\": 0,\n",
    "        \"h0\": writer_kwargs[\"sqrtSX\"] / 10,  # Fix amplitude at depth 10.\n",
    "        **pyfstat.isotropic_amplitude_distribution,\n",
    "        \"tref\": writer_kwargs[\"tstart\"],\n",
    "    },\n",
    ")\n",
    "\n",
    "for ind in range(num_signals):\n",
    "    params = signal_parameters_generator.draw()\n",
    "    writer_kwargs[\"outdir\"] = f\"PyFstat_example_data_ensemble/Signal_{ind}\"\n",
    "    writer_kwargs[\"label\"] = f\"Signal{ind}\"\n",
    "\n",
    "    writer = pyfstat.Writer(**writer_kwargs, **params)\n",
    "    writer.make_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "598f833d7ffd87ccd52f02f676ab7108c9e37681cbcc5e582981e8e544dc7b95"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}