{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mock traces example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "from radiocalibrationtoolkit import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "piko = 1e-12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some global plot settings\n",
    "plt.rcParams[\"axes.labelweight\"] = \"bold\"\n",
    "plt.rcParams[\"font.weight\"] = \"bold\"\n",
    "plt.rcParams[\"font.size\"] = 16\n",
    "plt.rcParams[\"legend.fontsize\"] = 12\n",
    "\n",
    "plt.rcParams[\"xtick.major.width\"] = 2\n",
    "plt.rcParams[\"ytick.major.width\"] = 2\n",
    "\n",
    "plt.rcParams[\"xtick.major.size\"] = 5\n",
    "plt.rcParams[\"ytick.major.size\"] = 5\n",
    "\n",
    "plt.rcParams[\"xtick.labelsize\"] = 14\n",
    "plt.rcParams[\"ytick.labelsize\"] = 14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read HW response\n",
    "hw_file_path = \"./antenna_setup_files/HardwareProfileList_realistic.xml\"\n",
    "hw_dict = read_hw_file(hw_file_path, interp_args={\"fill_value\": \"extrapolate\"})\n",
    "\n",
    "hw_reponse_1 = hw_dict[\"RResponse\"][\"LNA\"]\n",
    "hw_reponse_2 = hw_dict[\"RResponse\"][\"digitizer\"]\n",
    "hw_reponse_3 = hw_dict[\"RResponse\"][\"cable_fromLNA2digitizer\"]\n",
    "hw_reponse_4 = hw_dict[\"RResponse\"][\"impedance_matching_EW\"]\n",
    "\n",
    "\n",
    "# merge all hw responses to one function\n",
    "def hw_response_func(x):\n",
    "    return dB2PowerAmp(\n",
    "        hw_reponse_1(x) + hw_reponse_2(x) + hw_reponse_3(x) + hw_reponse_4(x)\n",
    "    )\n",
    "\n",
    "\n",
    "# impedance function\n",
    "impedance_func = hw_dict[\"IImpedance\"][\n",
    "    \"antenna_EW\"\n",
    "]\n",
    "\n",
    "# read sidereal voltage square spectral density\n",
    "sidereal_voltage2_density_DF = pd.read_csv(\n",
    "    \"./voltage2_density/voltage2_density_Salla_EW_GSM16.csv\",\n",
    "    index_col=0,\n",
    ")\n",
    "sidereal_voltage2_density_DF.columns = sidereal_voltage2_density_DF.columns.astype(\n",
    "    float\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "frequencies = np.linspace(0,125,126)\n",
    "ax.plot(frequencies, hw_reponse_1(frequencies), label='LNA')\n",
    "ax.plot(frequencies, hw_reponse_2(frequencies), label='digitizer')\n",
    "ax.plot(frequencies, hw_reponse_3(frequencies), label='cable')\n",
    "ax.plot(frequencies, hw_reponse_4(frequencies), label='imp. match')\n",
    "ax.set_xlabel('frequency [MHz]')\n",
    "ax.set_ylabel('dB')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "frequencies = np.linspace(0,125,126)\n",
    "ax.plot(frequencies, impedance_func(frequencies), label='impedance')\n",
    "ax.set_xlabel('frequency [MHz]')\n",
    "ax.set_ylabel('$\\Omega$')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power obtained from the spectra calculated by applying the Discrete Fourier Transform (DFT) to time traces is calculated as follows:\n",
    "\n",
    " \\begin{equation}\n",
    "P_s = 2\\frac{1}{T} \\sum_{k=f}^{f+\\delta f} \\frac{|X(k)|^2}{R(f)} \\Delta f\n",
    "\\end{equation}\n",
    "\n",
    "where $\\Delta t$ is the sampling time, $T$ the time trace length, $N$ the number of samples, $f_s$ sampling frequency and $R$ the impedance,\n",
    "\n",
    "where\n",
    "\n",
    " \\begin{equation}\n",
    "X(k) = X(k)_{DFT} \\Delta t = \\frac{X(k)_{DFT}}{fs}\n",
    "\\end{equation}\n",
    "\n",
    "where $X(k)_{DFT}$ is \"one-sided\" absolute spectrum obtained by applying discrete Fourier transform to the time trace. \n",
    "\n",
    "On the other hand, the power deposited to antenna by the galactic radio emission is calculated as:\n",
    "\n",
    " \\begin{equation}\n",
    "P_{sky}(t,f) = \\frac{k_{b}}{c^{2}} \\sum_{f} f^{2} \\int_{\\Omega} T_{sky}(t,f,\\theta,\\phi) \\frac{|H(f,\\theta,\\phi)|^{2}Z_{0}}{R(f)} \\Delta f d\\Omega  \n",
    "\\end{equation} \n",
    "\n",
    "where $Z_0$ is the vacuum impedance, $R_r$ is the antenna impedance, $k_b$ the Boltzmann constant and $c$ is the speed of light.\n",
    "\n",
    "Denoting power spectral density as \n",
    "\n",
    " \\begin{equation}\n",
    "S_{P}(t,f) = \\frac{k_{b}}{c^{2}} f^{2} \\int_{\\Omega} T_{sky}(t,f,\\theta,\\phi) \\frac{|H(f,\\theta,\\phi)|^{2}Z_{0}}{R(f)} d\\Omega  \n",
    "\\end{equation} \n",
    "\n",
    "we can rewrite the equation as\n",
    "\n",
    " \\begin{equation}\n",
    "P_{sky}(t,f) = \\sum_{f} S_{P}(t,f)  \\Delta f\n",
    "\\end{equation} \n",
    "\n",
    "Further, the voltage square density  $S_{\\nu}(t,f)$ we defined as  \n",
    "\n",
    " \\begin{equation}\n",
    "S_{\\nu}(t,f) = \\frac{S_{P}(t,f)}{R(f)}\n",
    "\\end{equation} \n",
    "\n",
    "or as:\n",
    "\n",
    " \\begin{equation}\n",
    "S_{\\nu}(t,f) = U_{\\nu}^2 \n",
    "\\end{equation} \n",
    "\n",
    "Equating the two definitions of power yields:\n",
    "\n",
    " \\begin{equation}\n",
    "2\\frac{1}{T} \\sum_{k=f}^{f+\\delta f} \\frac{|X(k)|^2}{R(f)} \\Delta f = \n",
    "\\frac{k_{b}}{c^{2}} \\sum_{f} f^{2} \\int_{\\Omega} T_{sky}(t,f,\\theta,\\phi) \\frac{|H(f,\\theta,\\phi)|^{2}Z_{0}}{R_{r}} \\Delta f d\\Omega  \n",
    "\\end{equation}\n",
    " \n",
    "\n",
    " \\begin{equation}\n",
    "2\\frac{1}{T} \\frac{|X(k)|^2}{R(f)}  = \n",
    "\\frac{k_{b}}{c^{2}} f^{2} \\int_{\\Omega} T_{sky}(t,f,\\theta,\\phi) \\frac{|H(f,\\theta,\\phi)|^{2}Z_{0}}{R(f)}  d\\Omega  \n",
    "\\end{equation}\n",
    "\n",
    " \\begin{equation}\n",
    "2\\frac{1}{T} |X(k)|^2 = \n",
    "\\frac{k_{b}}{c^{2}} f^{2} Z_{0}  \\int_{\\Omega} T_{sky}(t,f,\\theta,\\phi) |H(f,\\theta,\\phi)|^{2} d\\Omega  \n",
    "\\end{equation}\n",
    "\n",
    " \\begin{equation}\n",
    "2\\frac{1}{T} |X(k)|^2 = S_{\\nu}(t,f)  \n",
    "\\end{equation}\n",
    "\n",
    " \\begin{equation}\n",
    "|X(k)|^2 = \\frac{1}{2} T  S_{\\nu}(t,f)\n",
    "\\end{equation}\n",
    "\n",
    " \\begin{equation}\n",
    "|X(k)_{DFT}|^2 \\Delta t^2 = \n",
    "\\frac{1}{2} T  S_{\\nu}(t,f)\n",
    "\\end{equation}\n",
    "\n",
    " \\begin{equation}\n",
    "|X(k)_{DFT}|^2  = \n",
    "\\frac{T}{\\Delta t^2} \\frac{1}{2}  S_{\\nu}(t,f)  \n",
    "\\end{equation}\n",
    "\n",
    " \\begin{equation}\n",
    "|X(k)_{DFT}|^2  = \n",
    "\\frac{N f_{s}^2}{ f_s} \\frac{1}{2}   S_{\\nu}(t,f) \n",
    "\\end{equation}\n",
    "\n",
    " \\begin{equation}\n",
    "|X(k)_{DFT}|^2  =  \\frac{N f_{s}}{2}   S_{\\nu}(t,f)  \n",
    "\\end{equation}\n",
    "\n",
    " \\begin{equation}\n",
    "|X(k)_{DFT}|  = \\sqrt{\\frac{N f_{s} }{2}}   \\sqrt{S_{\\nu}(t,f)}  \n",
    "\\end{equation}\n",
    "\n",
    "A small reminder of the system parameter relations:\n",
    " \n",
    " \\begin{equation}\n",
    "\\Delta f = \\frac{fs}{N} = \\frac{1}{N\\Delta t} = \\frac{1}{T} \n",
    "\\end{equation}\n",
    "\n",
    " \\begin{equation}\n",
    "fs = \\frac{1}{\\Delta t}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thermal noise:\n",
    "\n",
    " \\begin{equation}\n",
    "v^2 = 4 k_B T R\n",
    "\\end{equation}\n",
    "\n",
    "where $k_B$ is the Boltzmann constant, $T$ is the temperature and $R$ is the impedance.\n",
    "\n",
    "source: https://en.wikipedia.org/wiki/Johnson%E2%80%93Nyquist_noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mock_trace_generator = Mock_trace_generator(\n",
    "    sidereal_voltage2_density_DF=sidereal_voltage2_density_DF,\n",
    "    hw_response_func=hw_response_func,\n",
    "    impedance_func=impedance_func,\n",
    "    voltage2ADC=2048,\n",
    "    time_trace_size=2048,\n",
    "    sampling_frequency_MHz=250,\n",
    ")\n",
    "freq_MHz_bins = mock_trace_generator.get_frequency_bins()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mock_traces_DF = mock_trace_generator.generate_mock_trace(5, nbi = {\"67\": 1}, nbi_err=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(mock_traces_DF.columns.values[2:], mock_traces_DF.iloc[1,2:])\n",
    "ax.set_xlabel('time [$\\mu$s]')\n",
    "ax.set_ylabel('amplitude [ADC]')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "additional_noise = 5e-4*piko\n",
    "\n",
    "debug_spectra_dict = mock_trace_generator.generate_mock_trace(\n",
    "    1, lst=15, temp_celsius=4, nbi={\"67\": 1}, nbi_err=0.3, return_debug_dict=True, additional_noise=additional_noise\n",
    ")[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "debug_spectra_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = debug_spectra_dict[\"freq_MHz_bins\"]\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_title(\"Signal at antenna\")\n",
    "\n",
    "ax.plot(\n",
    "    x,\n",
    "    debug_spectra_dict[\"sidereal_voltage2_density\"],\n",
    "    lw=3,\n",
    "    label=\"galactic signal\",\n",
    "    color=\"b\",\n",
    ")\n",
    "ax.plot(\n",
    "    x,\n",
    "    debug_spectra_dict[\"thermal_noise_voltage2_density\"],\n",
    "    lw=3,\n",
    "    label=\"thermal noise\",\n",
    "    color=\"green\",\n",
    ")\n",
    "\n",
    "ax.plot(\n",
    "    x,\n",
    "    debug_spectra_dict[\"sidereal_voltage2_density\"]\n",
    "    + debug_spectra_dict[\"thermal_noise_voltage2_density\"],\n",
    "    label=\"galactic signal + thermal noise\",\n",
    "    lw=3,\n",
    "    color=\"r\",\n",
    ")\n",
    "\n",
    "ax.plot(\n",
    "    x,\n",
    "    debug_spectra_dict[\"total_voltage2_density_scattered\"],\n",
    "    alpha=0.4,\n",
    "    label=\"galactic signal + thermal noise scattered\",\n",
    "    color=\"r\",\n",
    ")\n",
    "\n",
    "ax.set_xlabel(\"frequency [MHz]\")\n",
    "ax.set_ylabel(\"amplitude [V/$\\sqrt{Hz}]$\")\n",
    "ax.legend()\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_title(\"Signal in hardware\")\n",
    "\n",
    "volts2adc = 2048\n",
    "ax.plot(\n",
    "    x,\n",
    "    debug_spectra_dict[\"spectrum_voltage_density_in_HW_with_additional_noise_with_NBI\"] * volts2adc,\n",
    "    alpha=0.5,\n",
    "    label=\"spectrum + extra noise + NBI\",\n",
    ")\n",
    "ax.plot(\n",
    "    x,\n",
    "    debug_spectra_dict[\"spectrum_voltage_density_in_HW_with_additional_noise\"] * volts2adc,\n",
    "    alpha=0.5,\n",
    "    label=\"spectrum + extra noise\",\n",
    ")\n",
    "ax.plot(\n",
    "    x,\n",
    "    debug_spectra_dict[\"spectrum_voltage_density_in_HW\"] * volts2adc,\n",
    "    alpha=1,\n",
    "    label=\"spectrum\",\n",
    ")\n",
    "\n",
    "ax.set_xlabel(\"frequency [MHz]\")\n",
    "ax.set_ylabel(\"amplitude [ADC]\")\n",
    "ax.set_ylim(0, 360)\n",
    "ax.legend()"
   ]
  }
 ],
 "metadata": {
  "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.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}