{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import constants\n",
    "from scipy import signal as sigproc\n",
    "from scipy.fft import next_fast_len"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shot_no = 36880"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data loading\n",
    "\n",
    "The data is read from ISF files, using the [`isfreader`](https://github.com/justengel/isfreader) library. It is installed bellow automatically if not already installed.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install isfreader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import isfreader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = np.DataSource('/tmp')   # temporary sotarge for downloaded files\n",
    "data_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/Interferometry/DAS_raw_data_dir/ch{ch_id}.isf'\n",
    "scalars_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/Results/{name}'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_scalar(shot_no, name):\n",
    "    return float(ds.open(scalars_URL.format(shot_no=shot_no, name=name)).read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_plasma_start = get_scalar(shot_no, 't_plasma_start')\n",
    "t_plasma_end = get_scalar(shot_no, 't_plasma_end')\n",
    "is_plasma = get_scalar(shot_no, 'is_plasma')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_channel(shot_no, chanel_id, channel_name=None, in_ms=True):\n",
    "    fname = ds.open(data_URL.format(shot_no=shot_no, ch_id=chanel_id)).name\n",
    "    data = isfreader.read_file(fname)\n",
    "    signal = pd.Series(data[:,1], name=channel_name,\n",
    "                       index=pd.Index(data[:,0], name='time'))\n",
    "    if in_ms:\n",
    "        signal.index *= 1e3  # s -> ms\n",
    "    return signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixer = load_channel(shot_no, 1, 'mixer')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_saw = load_channel(shot_no, 3, 'ref_saw')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sampling frequency in kHz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_s = mixer.size / (mixer.index[-1] - mixer.index[0])  # ms -> kHz\n",
    "f_s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral analysis of signals\n",
    "\n",
    "The mixer signals is a base sine wave (the envelope of the mixing) at a frequency close to 500 kHz. The reference saw-tooth frequency sweeping wave has the same base frequency, but with a trail of harmonics forming the sharp saw-tooth shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_spectrum(signal, target_dfreq=10):\n",
    "    nperseg = int(f_s / target_dfreq) \n",
    "    f, psd = sigproc.welch(signal.values, fs=f_s,\n",
    "                           nperseg=nperseg, nfft=next_fast_len(nperseg))\n",
    "    \n",
    "    f_max = f[psd.argmax()]\n",
    "    ax = plt.gca()\n",
    "    ax.plot(f, psd, label=signal.name)\n",
    "    ax.set(xlabel='frequency [kHz]', ylabel='power spectral density [V^2]')\n",
    "    return f_max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_spectrum(ref_saw)\n",
    "f_base = show_spectrum(mixer)\n",
    "plt.axvline(f_base, label=f'base frequency f={f_base:.0f} kHz', color='C3')\n",
    "plt.loglog()\n",
    "plt.legend()\n",
    "f_base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Base band signal analysis\n",
    "\n",
    "The difference between the phase of the reference and mixer signal is proportional to the line-integrated electron density. Therefore, it is necessary to obtain the \"phase\" of both signals.\n",
    "\n",
    "The instantanous phase and amplitude of the base signal can be inferred only for the base bad, i.e. by removing higher and lower frequencies (i.e. extracting the base sine wave from the reference saw-tooth signal)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_band_hwidth = 50 # kHz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_band_filter = sigproc.iirfilter(8, \n",
    "                    [f_base - base_band_hwidth, f_base+base_band_hwidth],\n",
    "                                     fs=f_s,\n",
    "                                    btype='bandpass', output='sos')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def freq_filter(signal, sos_filter):\n",
    "    return pd.Series(name=signal.name+'_filtered', index=signal.index,\n",
    "        data=sigproc.sosfiltfilt(sos_filter, signal.values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixer_filtered = freq_filter(mixer, base_band_filter)\n",
    "ref_saw_filtered = freq_filter(ref_saw, base_band_filter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The instantanous phase and amplitude of this base band signal can be obtained from the [analytic signal](https://en.wikipedia.org/wiki/Analytic_signal) using the Hilbert transform. For convenience of processing and visualization, the phase is unwrapped ($\\pi$ discontinuities are shifted) and it's time derivative, the instantaneous frequency, is calucated as well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inst_phase_amp(signal):\n",
    "    analytic_signal = sigproc.hilbert(signal.values)\n",
    "    inst_phase = pd.Series(np.unwrap(np.angle(analytic_signal)), \n",
    "                      index=signal.index, \n",
    "                     name=signal.name+'_inst_phase')\n",
    "    # remove base band phase evolution trend\n",
    "    inst_freq = pd.Series(np.gradient(inst_phase, inst_phase.index), \n",
    "                          index=inst_phase.index, name=signal.name+'_inst_freq')\n",
    "    inst_amp = pd.Series(np.abs(analytic_signal), index=signal.index, \n",
    "                     name=signal.name+'_inst_amp')\n",
    "    return inst_phase, inst_freq, inst_amp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixer_phase, mixer_freq, mixer_amp = inst_phase_amp(mixer_filtered)\n",
    "ref_saw_phase, ref_saw_freq, ref_saw_amp = inst_phase_amp(ref_saw_filtered)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convenience of visualization it is convenient to subtract the phase evolution of the base nad signal from the phases (which is just a linear trend). This base trend vanishes in the final phase difference calculation anyways as long as the same trend is removed. Therefore, the referense saw-tooth signals is used for the trend estimation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega_base_ref = ref_saw_freq.median()\n",
    "mixer_dfreq = mixer_freq - omega_base_ref\n",
    "mixer_dphase = mixer_dfreq.cumsum() / f_s\n",
    "ref_saw_dphase = (ref_saw_freq - omega_base_ref).cumsum() / f_s\n",
    "omega_base_ref / (2*np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixer_freq.loc[-1:0].median() / (2*np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The signal boundaries may be distorted by the filter, so they are not used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary_Dt = 0.5  # in ms\n",
    "in_boundary = slice(mixer_dfreq.index[0]+boundary_Dt, \n",
    "                    mixer_dfreq.index[-1] - boundary_Dt)\n",
    "in_boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_saw_dphase.loc[in_boundary].plot()\n",
    "mixer_dphase.loc[in_boundary].plot()\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The phase difference between the interferometer signal with respect to the reference signal is proportional to the line-averaged density. An arbitrary initial phase difference is subtracted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dphase = pd.Series(ref_saw_dphase - mixer_dphase, \n",
    "                   index=mixer.index, name='dphase')\n",
    "dphase -= dphase.loc[-1:0].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2, 1, sharex=True, figsize=(8, 8))\n",
    "plt.sca(axs[0])\n",
    "mixer_amp.loc[in_boundary].plot()\n",
    "axs[0].set(ylabel='mixer amplitude')\n",
    "plt.sca(axs[1])\n",
    "dphase.loc[in_boundary].plot()\n",
    "if not is_plasma:\n",
    "    plt.savefig('icon-fig.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimation of the electron density\n",
    "\n",
    "The ordinary wave (O-mode) with a carrier frequency $\\omega$ traveling through a collisionless plasma with the plasma frequency $\\omega_{p} = \\sqrt{\\frac{n_e e^2}{\\epsilon_0 m_e}}$ has a refractive index $$N_O=\\sqrt{1-\\frac{\\omega_p^2}{\\omega^2}}$$ \n",
    "Under the assumption that the carrier wave frequency is much larger than the plasma frequency $\\omega>>\\omega_p$ this formula can be expanded into a Taylor series as $$N_O\\approx 1-\\frac{\\omega_p^2}{2\\omega^2}$$\n",
    "A wave traveling through a medium with a refractive index $N(l)$ accumulates a total phase shift $\\varphi = \\frac{\\omega}{c} \\int N(l) \\mathrm{d}l$. Therefore, in comparison to a wave traveling in vacuum (or clear air) with $N\\approx 1$, the  wave traveling through the plasma over a chord with length $L$ accumulates a relative phase shift of $$\\Delta \\varphi = \\frac{e^2}{2\\omega c\\epsilon_0 m_e}\\int\\limits_L n_e(l) \\mathrm{d}l$$\n",
    "Therefore, it is possible to estimate the line-averaged density $\\bar n_e = \\frac{1}{L} \\int\\limits_L n_e(l) \\mathrm{d}l$ from the detected phase shift between the reference and mixer signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega_carrier = 2*np.pi*71e9  # 71 GHz microwave generator\n",
    "a = 0.085   # limiter radius\n",
    "L = 2*a  # central plasma chord estimate\n",
    "prop_const = constants.elementary_charge**2 / (\n",
    "    2*omega_carrier* constants.speed_of_light*constants.epsilon_0*constants.m_e)\n",
    "if is_plasma:\n",
    "    ne_lav = dphase / prop_const / L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only save every 100-th point as there is not much useful signal above 2 MHz\n",
    "save_sl = slice(0, t_plasma_end+5, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if is_plasma:\n",
    "    ne_lav.loc[save_sl].to_csv('ne_lav.csv', header=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if is_plasma:\n",
    "    fig = plt.figure(figsize=(6, 4.5), dpi=300)\n",
    "    (ne_lav.loc[save_sl] / 1e19).plot()\n",
    "    for t in (t_plasma_start, t_plasma_end):\n",
    "        plt.axvline(t, color='k', linestyle='--')\n",
    "    plt.gca().set(xlabel='time [ms]', ylim=(0,None),\n",
    "                  ylabel='$\\\\bar n_e$ [10$^{19}$ m$^{-3}$]')\n",
    "    plt.savefig('icon-fig.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if is_plasma:\n",
    "    ne_lav_plasma = ne_lav.loc[t_plasma_start:t_plasma_end]\n",
    "    for m in ('mean', 'max'):\n",
    "        v = getattr(ne_lav_plasma, m)()\n",
    "        with open(f'ne_lav_{m}', 'w') as f:\n",
    "            f.write(f'{v:.3e}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}