{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plasma detection from plasma current signal\n",
    "\n",
    "This is a quick plasma detection code from the plasma current signal.\n",
    "This procedure is run as soon as possible in order to provide the information on plasma existence to other routines.\n",
    "\n",
    "(author: L. Lobko)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import os\n",
    "from scipy import signal as sigproc\n",
    "from scipy import integrate\n",
    "import matplotlib.pyplot as plt"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": "shot_no = 47567",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "def update_db_current_shot(field_name, value):\n",
    "    os.system('export PGPASSWORD=`cat /golem/production/psql_password`;psql -c \"UPDATE operation.discharges SET '+field_name+'='+str(value)+'WHERE shot_no IN(SELECT max(shot_no) FROM operation.discharges)\" -q -U golem golem_database')\n",
    "    os.system('export PGPASSWORD=`cat /golem/production/psql_password`;psql -c \"UPDATE diagnostics.basicdiagnostics SET '+field_name+'='+str(value)+'WHERE shot_no IN(SELECT max(shot_no) FROM diagnostics.basicdiagnostics)\" -q -U golem golem_database')"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "os.makedirs('Results', exist_ok=True)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "def save_scalar(phys_quant, value, format_str='%.3f'):\n",
    "    with open(\"Results/\"+phys_quant, 'w') as f:\n",
    "        f.write(format_str % value)\n",
    "    update_db_current_shot(phys_quant,value)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "ds = np.DataSource('/tmp')  # temporary storage for downloaded files\n",
    "data_rog_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/U_RogCoil.csv'\n",
    "data_ULoop_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/U_Loop.csv'\n",
    "t_cd_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/t_cd_discharge_request'\n",
    "K_RogowskiCoil_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Production/Parameters/SystemParameters/K_RogowskiCoil'\n",
    "L_chamber_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Production/Parameters/SystemParameters/L_chamber'\n",
    "R_chamber_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Production/Parameters/SystemParameters/R_chamber'"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load data"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "def load_data(shot_no, signal_URL):\n",
    "    fname = ds.open(signal_URL.format(shot_no=shot_no)).name\n",
    "    data = np.loadtxt(fname, delimiter=',')\n",
    "    # data[:, 0] = data[:, 0] * 1e3 # from micros to ms\n",
    "    return data"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def load_param(shot_no, param):\n",
    "    data = float(ds.open(param.format(shot_no=shot_no)).read())\n",
    "    # data = data * 1e-3\n",
    "    return data"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "U_rogcoil = load_data(shot_no, data_rog_URL)\n",
    "U_Loop = load_data(shot_no, data_ULoop_URL)\n",
    "t_cd = load_param(shot_no, t_cd_URL)\n",
    "K_RogowskiCoil = load_param(shot_no, K_RogowskiCoil_URL)\n",
    "L_chamber = load_param(shot_no, L_chamber_URL)\n",
    "R_chamber = load_param(shot_no, R_chamber_URL)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "U_rogcoil[:, 1] *= -1",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# fig, ax = plt.subplots()\n",
    "# ax.plot(U_Loop[:, 0], U_Loop[:, 1], 'b-')\n",
    "# plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def offset_remove(data):\n",
    "    x_size, y_size = data.shape\n",
    "    data_for_offset = data[0:int(x_size/100)]\n",
    "    offset = np.mean(data_for_offset[:, 1])\n",
    "    data[:, 1] -= offset\n",
    "    return data"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "U_rogcoil = offset_remove(U_rogcoil)\n",
    "U_Loop = offset_remove(U_Loop)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# fig, ax = plt.subplots()\n",
    "# ax.plot(U_Loop[:, 0], U_Loop[:, 1], 'b-')\n",
    "# plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Integrate signal"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "U_integrated = np.copy(U_rogcoil)\n",
    "U_integrated[:, 1] = (integrate.cumtrapz(U_rogcoil[:, 1], U_rogcoil[:, 0], initial=0))\n",
    "U_integrated[:, 1] *= K_RogowskiCoil"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# fig, ax = plt.subplots()\n",
    "# ax.plot(U_integrated[:, 0], U_integrated[:, 1], 'b-')\n",
    "# plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Calculate chamber current"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# def dIch_dt(t, Ich):\n",
    "#     return (U_l_func(t) - R_chamber * Ich) / L_chamber\n",
    "# \n",
    "# dIch_dt = (U_Loop - R_chamber *)/ L_chamber"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "Ich = np.copy(U_Loop)\n",
    "Ich[:, 1] = U_Loop[:, 1] / R_chamber"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# fig, ax = plt.subplots()\n",
    "# ax.plot(Ich[:, 0], Ich[:, 1], 'b-')\n",
    "# plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "Ip = np.copy(Ich)\n",
    "Ip[:, 1] = U_integrated[:, 1] - Ich[:, 1]"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# fig, ax = plt.subplots()\n",
    "# ax.plot(Ip[:, 0], Ip[:, 1], 'b-')\n",
    "# plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Smooth signal"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "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"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "Ip[:, 1] = smooth(Ip[:, 1], 200)\n",
    "Ip[:, 0] = Ip[:, 0] * 1e3"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# fig, ax = plt.subplots()\n",
    "# ax.plot(Ip[:, 0], Ip[:, 1], 'b-')\n",
    "# plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "def find_peaks(data):\n",
    "    peaks_indexes, _ = sigproc.find_peaks(data[:, 1], prominence=1e-1)\n",
    "    return np.vstack((data[peaks_indexes, 0], data[peaks_indexes, 1])).T"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate plasma boundaries from signal derivation"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "def calc_plasma_boundaries(data, position):\n",
    "    deriv = data.copy()\n",
    "    deriv[:, 1] = np.gradient(data[:, 1])\n",
    "    deriv[:, 1] = smooth(deriv[:, 1], 1000)\n",
    "    if position == 'start':\n",
    "        index = np.where(deriv[:, 1] >= np.max(deriv[:, 1])/5)\n",
    "        deriv = deriv[index]\n",
    "        return deriv[0, 0]\n",
    "    else:\n",
    "        deriv = np.abs(deriv)\n",
    "        # max_time = np.max(deriv[:, 0])\n",
    "        # deriv = deriv[deriv[:, 0] <= (max_time-0.5)]\n",
    "        peaks = find_peaks(deriv)\n",
    "        return (peaks[-1, 0]+0.5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Cut data before t_cd (before is no plasma)"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "Ip_plasma_check = Ip[Ip[:, 0] <= (t_cd/1000+5)]",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "Ip = Ip[Ip[:, 0] > (t_cd/1000)]",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# fig, ax = plt.subplots()\n",
    "# ax.plot(Ip_plasma_check[:, 0], Ip_plasma_check[:, 1], 'b-')\n",
    "# plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "if np.max(Ip_plasma_check[:, 1]) < 100:\n",
    "    print('No plasma in vacuum chamber.')\n",
    "    t_plasma_start = -1.0\n",
    "    t_plasma_end = -1.0\n",
    "else:\n",
    "    t_plasma_start = calc_plasma_boundaries(Ip, 'start')\n",
    "    print('Plasma starts at {:.2f} ms.'.format(t_plasma_start))\n",
    "    t_plasma_end = calc_plasma_boundaries(Ip, 'end')\n",
    "    print('Plasma ends at {:.2f} ms.'.format(t_plasma_end))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "b_plasma = int(t_plasma_start > 0 and t_plasma_end > 0)",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "if b_plasma:\n",
    "    t_plasma_duration = t_plasma_end - t_plasma_start\n",
    "    print('Plasma duration is {:.2f} ms.'.format(t_plasma_duration))\n",
    "else:\n",
    "    t_plasma_duration = -1.0  # convention instead of nan"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "plasma_endpoints = [t_plasma_start, t_plasma_end]",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "fig, axes = plt.subplots()\n",
    "axes.plot(Ip[:, 0], Ip[:, 1]/1000, label ='Ip calculated')\n",
    "for x in plasma_endpoints:\n",
    "    plt.axvline(x=x, color='black', linestyle='--')\n",
    "axes.set(xlabel='$time$ [ms]', ylabel='$I_p$ [kA]')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Save data"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "save_scalar(\"b_plasma\", b_plasma)\n",
    "save_scalar(\"t_plasma_start\", t_plasma_start)\n",
    "save_scalar(\"t_plasma_end\", t_plasma_end)\n",
    "save_scalar(\"t_plasma_duration\", t_plasma_duration)"
   ],
   "outputs": [],
   "execution_count": null
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}