{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IMPORT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import constants\n",
    "from scipy import signal as sigproc\n",
    "from scipy.fft import next_fast_len\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shot_no = 40485 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data loading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "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 storage for downloaded files\n",
    "data_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/Interferometry/LukLob/DAS_raw_data_dir/ch{ch_id}.isf'\n",
    "scalars_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/DetectPlasma/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):\n",
    "    fname = ds.open(data_URL.format(shot_no=shot_no, ch_id=chanel_id)).name\n",
    "    data = isfreader.read_file(fname)\n",
    "    data[:, 0] = data[:, 0] * 1e3\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixer = load_channel(shot_no, 1)\n",
    "ref_saw = load_channel(shot_no, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = mixer.shape\n",
    "f_s = x / (mixer[-1, 0] - mixer[0, 0])  # ms -> kHz\n",
    "f_s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral analysis of signals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mixer signals is a base sine wave (the envelope of the mixing) at a frequency close to 500 kHz. The reference\n",
    "saw-tooth frequency sweeping wave has the same base frequency, but with a trail of harmonics forming the sharp\n",
    "saw-tooth shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_spectrum(signal, target_dfreq=10):\n",
    "    nperseg = int(f_s / target_dfreq)\n",
    "    f, psd = sigproc.welch(signal[:, 1], fs=f_s, nperseg=nperseg, nfft=next_fast_len(nperseg))\n",
    "    return f, psd\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_saw_f, ref_saw_psd = calculate_spectrum(ref_saw)\n",
    "mixer_f, mixer_psd = calculate_spectrum(mixer)\n",
    "f_base = mixer_f[mixer_psd.argmax()]\n",
    "f_base_ref_saw = ref_saw_f[ref_saw_psd.argmax()]\n",
    "f_base  # [kHz]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_base_ref_saw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.set(xlabel='frequency [kHz]', ylabel='power spectral density [V^2]')\n",
    "ax.plot(ref_saw_f, ref_saw_psd, label='ref_saw')\n",
    "ax.plot(mixer_f, mixer_psd, label='mixer')\n",
    "plt.axvline(f_base, label=f'base frequency f={f_base:.0f} kHz', color='C3')\n",
    "ax.loglog()\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Base band signal analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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 instantaneous phase and amplitude of the base signal can be inferred only for the base band, 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\n",
    "base_band_filter = sigproc.iirfilter(8, [f_base - base_band_hwidth, f_base + base_band_hwidth], fs=f_s, btype='bandpass', output='sos')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def freq_filter(signal, sos_filter):\n",
    "    signal[:, 1] = sigproc.sosfiltfilt(sos_filter, signal[:, 1])\n",
    "    return signal"
   ]
  },
  {
   "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": [
    "cut 0.1 ms from the beginning and from the end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixer_filtered = mixer_filtered[(mixer_filtered[:, 0] < (mixer_filtered[-1, 0] - 0.1)) & ((mixer_filtered[0, 0] + 0.1) < mixer_filtered[:, 0])]\n",
    "ref_saw_filtered = ref_saw_filtered[(ref_saw_filtered[:, 0] < (ref_saw_filtered[-1, 0] - 0.1)) & ((ref_saw_filtered[0, 0] + 0.1) < ref_saw_filtered[:, 0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_peaks(data):\n",
    "    peaks_indexes, _ = sigproc.find_peaks(data[:, 1])\n",
    "    return np.vstack((data[peaks_indexes, 0], data[peaks_indexes, 1])).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initial_phase_shift(peaks, peaks_ref, mixer_filtered):\n",
    "    phase_mean = peaks[0, 0] - peaks_ref[0, 0]\n",
    "    mixer_filtered[:, 0] -= phase_mean\n",
    "    peaks[:, 0] -= phase_mean\n",
    "    return peaks, mixer_filtered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def differences(signal):\n",
    "    data = signal.copy()\n",
    "    x, y = data.shape\n",
    "    for i in range(1, int(x) - 1):\n",
    "        data[i, 0] = signal[i, 0] - signal[i - 1, 0]\n",
    "        data[i, 1] = signal[i + 1, 0] - signal[i, 0]\n",
    "\n",
    "    return np.hstack((signal, data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cut_edges(diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref):\n",
    "    diff_peaks_mixer = diff_peaks_mixer[(diff_peaks_mixer[0, 0] < diff_peaks_mixer[:, 0]) & (diff_peaks_mixer[:, 0] < diff_peaks_mixer[-1, 0])]\n",
    "    diff_peaks_ref = diff_peaks_ref[(diff_peaks_ref[0, 0] < diff_peaks_ref[:, 0]) & (diff_peaks_ref[:, 0] < diff_peaks_ref[-1, 0])]\n",
    "    peaks = peaks[(peaks[0, 0] < peaks[:, 0]) & (peaks[:, 0] < peaks[-1, 0])]\n",
    "    peaks_ref = peaks_ref[(peaks_ref[0, 0] < peaks_ref[:, 0]) & (peaks_ref[:, 0] < peaks_ref[-1, 0])]\n",
    "    return diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_basic_diff(diff_peaks_mixer, diff_peaks_ref):\n",
    "    basic_diff = np.mean(diff_peaks_ref[:, 2])\n",
    "    diff_peaks_mixer[:, 2] -= basic_diff\n",
    "    diff_peaks_ref[:, 2] -= basic_diff\n",
    "    return diff_peaks_mixer, diff_peaks_ref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sort_differencies(diff_peaks_mixer):\n",
    "    data = diff_peaks_mixer.copy()\n",
    "    x, y = data.shape\n",
    "    for i in range(0, x):\n",
    "        if data[i, 2] < 0:\n",
    "            data[i, 2] *= 0.9\n",
    "    data[:, 2] = abs(data[:, 2])\n",
    "\n",
    "    return data[data[:, 2].argsort()[::-1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def smooth(y, box_pts):\n",
    "    box = np.ones(box_pts)/box_pts\n",
    "    y_smooth = np.convolve(y, box, mode='same')\n",
    "    return y_smooth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def without_correction(mixer_filtered, ref_saw_filtered):\n",
    "    mixer_data = mixer_filtered.copy()\n",
    "    ref_saw_data = ref_saw_filtered.copy()\n",
    "    peaks = find_peaks(mixer_data)\n",
    "    peaks_ref = find_peaks(ref_saw_data)\n",
    "    peaks, mixer_data = initial_phase_shift(peaks, peaks_ref, mixer_data)\n",
    "    diff_peaks_mixer = differences(peaks)\n",
    "    diff_peaks_ref = differences(peaks_ref)\n",
    "    diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref = cut_edges(diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref)\n",
    "    diff_peaks_mixer, diff_peaks_ref = remove_basic_diff(diff_peaks_mixer, diff_peaks_ref)\n",
    "\n",
    "    return mixer_data, peaks, peaks_ref, diff_peaks_mixer, diff_peaks_ref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def basic_cycle(mixer_filtered, ref_saw_filtered):\n",
    "    peaks = find_peaks(mixer_filtered)\n",
    "    peaks_ref = find_peaks(ref_saw_filtered)\n",
    "    peaks, mixer_filtered = initial_phase_shift(peaks, peaks_ref, mixer_filtered)\n",
    "    diff_peaks_mixer = differences(peaks)\n",
    "    diff_peaks_ref = differences(peaks_ref)\n",
    "    diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref = cut_edges(diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref)\n",
    "    diff_peaks_mixer, diff_peaks_ref = remove_basic_diff(diff_peaks_mixer, diff_peaks_ref)\n",
    "    sorted_diff_mixer = sort_differencies(diff_peaks_mixer)\n",
    "    return peaks, peaks_ref, mixer_filtered, diff_peaks_mixer, diff_peaks_ref, sorted_diff_mixer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_nearest(array, value):\n",
    "    array = array.copy()\n",
    "    x, y = array.shape\n",
    "    array_data = np.zeros((x, 1))\n",
    "    for i in range(0, x):\n",
    "        array_data[i, 0] = array[i, 0]\n",
    "    idx = (np.abs(array_data[:, 0] - value)).argmin()\n",
    "    return array_data[idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_dphase_unchanged(peaks, peaks_ref):\n",
    "    x_peaks, y_peaks = peaks.shape\n",
    "    x_ref_peaks, y_ref_peaks = peaks_ref.shape\n",
    "\n",
    "    dphase = np.ones((min(x_peaks, x_ref_peaks), 2))\n",
    "    for i in range(0, int(len(dphase))):\n",
    "        dphase[i, 0] = peaks[i, 0] # TODO: zmenit na peaks_ref?\n",
    "        dphase[i, 1] = peaks[i, 0] - peaks_ref[i, 0]\n",
    "\n",
    "    return dphase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_lost_phase(peaks, peaks_ref, diff_peaks_mixer, diff_peaks_ref):\n",
    "    dphase = calc_dphase_unchanged(peaks, peaks_ref)\n",
    "    time_interval = 0.1 # ms\n",
    "    indexes = np.argwhere(dphase[:, 0] > (dphase[-1, 0] - time_interval))\n",
    "    data = dphase[indexes[:, 0]]\n",
    "    return 2*math.pi*f_base*np.average(data[:, 1]), dphase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def repair(mixer_filtered, ref_saw_filtered, lost_phase, dphase):\n",
    "    peaks, peaks_ref, mixer_filtered, diff_peaks_mixer, diff_peaks_ref, sorted_diff_mixer = basic_cycle(mixer_filtered, ref_saw_filtered)\n",
    "    # make repairing cycle\n",
    "    number_of_bad_peaks = round(lost_phase/(2*math.pi))\n",
    "    remainder = lost_phase - 2*math.pi*number_of_bad_peaks\n",
    "    x_peaks_mixer, y_peaks_mixer = peaks.shape\n",
    "    x_peaks_ref, y_peaks_ref = peaks_ref.shape\n",
    "    deriv = dphase.copy()\n",
    "    deriv[:, 1] = np.gradient(dphase[:, 1])\n",
    "    deriv[:, 1] = smooth((deriv[:, 1]), 10)\n",
    "    deriv_sort = deriv.copy()\n",
    "    # deriv_sort[:, 1] = abs(deriv_sort[:, 1])\n",
    "    deriv_sort = deriv_sort[deriv_sort[:, 1].argsort()[::-1]]\n",
    "    bad_peaks_indexes = np.empty((0, 1))\n",
    "    k = 0\n",
    "    l = 0\n",
    "    while k < number_of_bad_peaks:\n",
    "        index = np.argwhere(diff_peaks_mixer[:, 0] == deriv_sort[l, 0])\n",
    "        if (l > 0) and (abs((index[0, 0] - find_nearest(bad_peaks_indexes, index[0, 0]))) < 10):\n",
    "            l += 1\n",
    "        else:\n",
    "            bad_peaks_indexes = np.vstack((bad_peaks_indexes, index[0, 0]))\n",
    "            peaks_ref = np.delete(peaks_ref, index, 0)\n",
    "            k += 1\n",
    "            l += 1\n",
    "\n",
    "    if number_of_bad_peaks == 0:\n",
    "        for i in range(x_peaks_mixer):\n",
    "            if remainder > math.pi:\n",
    "                index = np.argwhere(diff_peaks_mixer[:, 0] == deriv_sort[i, 0])\n",
    "                peaks_ref = np.delete(peaks_ref, index, 0)\n",
    "                bad_peaks_indexes = np.vstack((bad_peaks_indexes, index[0, 0]))\n",
    "\n",
    "    return peaks, peaks_ref, bad_peaks_indexes, diff_peaks_mixer, deriv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_dphase(peaks, peaks_ref, repaired_mixer_peaks, repaired_ref_peaks):\n",
    "    x_peaks, y_peaks = peaks.shape\n",
    "    x_ref_peaks, y_ref_peaks = peaks_ref.shape\n",
    "\n",
    "    dphase = np.ones((min(x_peaks, x_ref_peaks), 2))\n",
    "    dphase[:, 1] *= 2*math.pi*f_base\n",
    "    for i in range(0, int(len(dphase))):\n",
    "        dphase[i, 0] = peaks_ref[i, 0]\n",
    "        dphase[i, 1] = peaks[i, 0] - peaks_ref[i, 0]\n",
    "\n",
    "    x_peaks, y_peaks = repaired_mixer_peaks.shape\n",
    "    x_ref_peaks, y_ref_peaks = repaired_ref_peaks.shape\n",
    "\n",
    "    dphase1 = np.ones((min(x_peaks, x_ref_peaks), 2))\n",
    "    for i in range(0, int(len(dphase1))):\n",
    "        dphase1[i, 0] = repaired_ref_peaks[i, 0]\n",
    "        dphase1[i, 1] = repaired_mixer_peaks[i, 0] - repaired_ref_peaks[i, 0]\n",
    "\n",
    "    dphase1[:, 1] *= 2*math.pi*f_base\n",
    "    return dphase, dphase1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixer_filtered_no_corr, peaks_no_corr, peaks_ref_no_corr, diff_peaks_mixer_no_corr, diff_peaks_ref_no_corr = without_correction(mixer_filtered, ref_saw_filtered)\n",
    "lost_phase, dphase_zero = calc_lost_phase(peaks_no_corr, peaks_ref_no_corr, diff_peaks_mixer_no_corr, diff_peaks_ref_no_corr)\n",
    "repaired_mixer_peaks, repaired_ref_peaks, bad_peaks_indexes, diff_peaks_mixer, deriv = repair(mixer_filtered, ref_saw_filtered, lost_phase, dphase_zero)\n",
    "dphase, dphase1 = calc_dphase(peaks_no_corr, peaks_ref_no_corr, repaired_mixer_peaks, repaired_ref_peaks)\n",
    "dphase1[:, 1] = smooth(dphase1[:, 1], 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimation of the electron density"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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 / (2 * omega_carrier * constants.speed_of_light * constants.epsilon_0 * constants.m_e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ne_lav = dphase1.copy()\n",
    "ne_lav[:, 1] = ne_lav[:, 1] * (1 / (prop_const*L))\n",
    "ne_lav = ne_lav[(ne_lav[:, 0] >= 0) & (ne_lav[:, 0] <= (t_plasma_end + 5))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FINAL FIGURE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4.5), dpi=300)\n",
    "ax.plot(ne_lav[:, 0], ne_lav[:, 1]/1e18)\n",
    "for t in (t_plasma_start, t_plasma_end):\n",
    "    plt.axvline(t, color='k', linestyle='--')\n",
    "ax.set(xlabel='time [ms]', ylim=(0, None), ylabel='$\\\\bar n_e$ [10$^{18}$ m$^{-3}$]')\n",
    "plt.savefig('icon-fig.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SAVE DATA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt('ne_lav.csv', ne_lav, delimiter=\",\")\n",
    "ne_lav_plasma = ne_lav[((t_plasma_start) < ne_lav[:, 0]) & (ne_lav[:, 0] < (t_plasma_end))]\n",
    "for m in ('mean', 'max'):\n",
    "    v = getattr(ne_lav_plasma[:, 1], m)()\n",
    "    with open(f'ne_lav_{m}', 'w') as f: f.write(f'{v:.3e}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}