{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plasma column position\n",
    "\n",
    "Adapted version of InnerStabilisation.ipynb by Daniela Kropáčková, Honza Stöckel and Vojtěch Svoboda.\n",
    "\n",
    "This notebook loads Mirnov coil data and calculates the plasma column position from them.\n",
    "\n",
    "### Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "#delete files from the previous discharge\n",
    "if os.path.exists('icon-fig.png'):\n",
    "    os.remove('icon-fig.png')\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import holoviews as hv\n",
    "hv.extension('bokeh')\n",
    "import hvplot.pandas\n",
    "from scipy import integrate, signal, interpolate\n",
    "import pandas as pd\n",
    "import requests\n",
    "from IPython.display import Markdown"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set discharge number\n",
    "\n",
    "To calculate the plasma position, $B_t$ contribution must be subtracted from the Mirnov coil signal. This contribution can be acquired by two ways: from a dedicated vacuum shot with the same toroidal magnetic field, or it can be estimated from the standard diagnostics. If a suitable vacuum shot is available, its number is stored in the variable `vacuum_shot`. If it isn't available, the `vacuum_shot` variable shall be set to `False` and the $B_t$ contribution shall be estimated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_URL = \"http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/LimiterMirnovCoils/{identifier}.csv\"  #Mirnov coils and quadrupole\n",
    "BDdata_URL = \"http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/{identifier}.csv\" #BD = basic diagnostic\n",
    "# data_URL = \"http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/022Daniela0InnerStabilization/{identifier}.csv\"  \n",
    "data_URL = \"http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/LimiterMirnovCoils/DAS_raw_data_dir/Nidatap_6133.lvm\" #NI Standart (DAS)\n",
    "\n",
    "from urllib.request import urlopen\n",
    "\n",
    "def get_file(url, shot, silent=False):\n",
    "    URL = 'http://golem.fjfi.cvut.cz/shots/%i/%s' % (shot, url)\n",
    "    if not silent:\n",
    "        print(URL)\n",
    "    f = urlopen(URL)\n",
    "    try:\n",
    "        return np.loadtxt(f)\n",
    "    except ValueError: # the data couldn't be converted to a row of numbers\n",
    "        print('The data could not be read! Using default values.')\n",
    "        return np.array([35022])\n",
    "\n",
    "shot_no = 35805\n",
    "#shot_no = 35021\n",
    "vacuum_shot = int(get_file('/Diagnostics/LimiterMirnovCoils/Parameters/vacuum_shot', 0))\n",
    "#vacuum_shot = 35022 #number of the vacuum shot or 'False'\n",
    "print('Vacuum shot: %i' % vacuum_shot)\n",
    "\n",
    "ds = np.DataSource(destpath='')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data access"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def open_remote(shot_no, identifier, url_template=data_URL):\n",
    "    return ds.open(url_template.format(shot_no=shot_no, identifier=identifier))\n",
    "\n",
    "def read_signal(shot_no, identifier, url = data_URL, data_type='csv'): \n",
    "    file = open_remote(shot_no, identifier, url)\n",
    "    if data_type == 'lvm':\n",
    "        channels = ['Time', 'mc1', 'mc5', 'mc9', 'mc13', '5', 'RogQuadr', '7', '8', '9', '10', '11', '12']\n",
    "        return pd.read_table(file, sep = '\\t', names=channels)\n",
    "    else:\n",
    "        return pd.read_csv(file, names=['Time', identifier],\n",
    "                     index_col = 'Time', squeeze=True)\n",
    "\n",
    "def get_basic_diag_data(shot):\n",
    "    url = ( 'http://golem.fjfi.cvut.cz/shots/{}/Diagnostics/BasicDiagnostics/'\n",
    "            'basig_diagnostics_processed.csv'.format(shot) )\n",
    "    print(url)\n",
    "    return pd.read_csv(url, header=0, index_col='Time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Mirnov coil signal processing: integration and $B_t$ elimination\n",
    "\n",
    "Load the data of a given Mirnov coil, integrate it and subtract the $B_t$ contribution from them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def elimination (shot_no, identifier, vacuum_shot = False):\n",
    "    #load data \n",
    "    mc = (read_signal(shot_no, identifier, data_URL, 'lvm'))\n",
    "    mc = mc.set_index('Time')\n",
    "    mc = mc.replace([np.inf, -np.inf, np.nan], value = 0)\n",
    "    mc = mc[identifier]\n",
    "    \n",
    "    konst = 1/(3.8e-03)\n",
    "\n",
    "    if vacuum_shot == False: \n",
    "        signal_start = mc.index[0]\n",
    "        length = len(mc)\n",
    "        basic_diag = get_basic_diag_data(shot_no)\n",
    "        Bt = basic_diag['Bt']\n",
    "        if len(Bt)>len(mc):\n",
    "            Bt = Bt.iloc[:length]            \n",
    "        if len(Bt)<len(mc):\n",
    "            mc = mc.iloc[:len(Bt)]\n",
    "        \n",
    "        if identifier == 'mc1':\n",
    "            k=300\n",
    "        elif identifier == 'mc5':\n",
    "            k= 14\n",
    "        elif identifier == 'mc9':\n",
    "            k = 31\n",
    "        elif identifier == 'mc13':\n",
    "            k = -100 \n",
    "        \n",
    "        mc_vacuum = Bt/k\n",
    "    else:\n",
    "        mc_vacuum = (read_signal(vacuum_shot, identifier, data_URL, 'lvm'))\n",
    "        mc_vacuum = mc_vacuum.set_index('Time')\n",
    "        mc_vacuum = mc_vacuum[identifier]\n",
    "        mc_vacuum -= mc_vacuum.loc[:0.9e-3].mean()    #remove offset\n",
    "        mc_vacuum = mc_vacuum.replace([np.inf, -np.inf], value = 0)\n",
    "        \n",
    "        mc_vacuum = pd.Series(integrate.cumtrapz(mc_vacuum, x=mc_vacuum.index, initial=0) * konst,\n",
    "                    index=mc_vacuum.index*1000, name= identifier)    #integration\n",
    "\n",
    "    mc -= mc.loc[:0.9e-3].mean()  #remove offset\n",
    "       \n",
    "    mc = pd.Series(integrate.cumtrapz(mc, x=mc.index, initial=0) * konst,\n",
    "                    index=mc.index*1000, name= identifier)    #integration\n",
    "    \n",
    "    #Bt elimination\n",
    "    mc_vacuum = np.array(mc_vacuum) \n",
    "    mc_elim = mc - mc_vacuum\n",
    "    \n",
    "    return mc_elim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find plasma start and end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plasma_start = float(get_file('Diagnostics/BasicDiagnostics/Results/t_plasma_start', shot_no))\n",
    "plasma_end = float(get_file('Diagnostics/BasicDiagnostics/Results/t_plasma_end', shot_no))\n",
    "\n",
    "print ('Plasma start =', round(plasma_start, 3), 'ms')\n",
    "print ('Plasma end =', round(plasma_end, 3), 'ms')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Horizontal plasma position $\\Delta r$ calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def horpol(shot_no, vacuum_shot=False):\n",
    "    mc1 = elimination(shot_no, 'mc1', vacuum_shot)\n",
    "    mc9 = elimination (shot_no, 'mc9', vacuum_shot)\n",
    "    \n",
    "    b = 93\n",
    "    \n",
    "    r = ((mc1-mc9)/(mc1+mc9))*b\n",
    "    r = r.replace([np.nan], value = 0)\n",
    "    return r.loc[plasma_start:plasma_end]\n",
    "\n",
    "r = horpol(shot_no, vacuum_shot)\n",
    "ax = r.plot()\n",
    "ax.set(ylim=(-85,85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$\\Delta$r [mm]',\n",
    "       title = 'Horizontal plasma position #{}'.format(shot_no))\n",
    "ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vertical plasma position $\\Delta z$ calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vertpol(shot_no, vacuum_shot = False):\n",
    "    mc5 = elimination(shot_no, 'mc5', vacuum_shot)\n",
    "    mc13 = elimination (shot_no, 'mc13', vacuum_shot)\n",
    "    \n",
    "    b = 93\n",
    "    \n",
    "    z = ((mc5-mc13)/(mc5+mc13))*b\n",
    "    z = z.replace([np.nan], value = 0)\n",
    "    return z.loc[plasma_start:plasma_end]\n",
    "\n",
    "z = vertpol (shot_no, vacuum_shot)\n",
    "ax = z.plot()\n",
    "ax.set(ylim=(-85, 85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$\\Delta$z [mm]', title = 'Vertical plasma position #{}'.format(shot_no))\n",
    "ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plasma column radius $a$ calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plasma_radius(shot_no, vacuum_shot=False):\n",
    "    r = horpol(shot_no, vacuum_shot) \n",
    "    z = vertpol(shot_no, vacuum_shot) \n",
    "    \n",
    "    a0 = 85\n",
    "    a = a0 - np.sqrt((r**2)+(z**2)) \n",
    "    a = a.replace([np.nan], value = 0)\n",
    "    return a.loc[plasma_start:plasma_end]\n",
    "\n",
    "\n",
    "a = plasma_radius(shot_no,vacuum_shot)\n",
    "ax = a.plot()\n",
    "ax.set(ylim=(0,85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$a$ [mm]',\n",
    "        title = 'Plasma column radius #{}'.format(shot_no))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plasma_time = []\n",
    "t = 0\n",
    "for i in a:\n",
    "    if i>85 or i <0:\n",
    "        a = a.replace(i, value = 0)\n",
    "    else:\n",
    "\n",
    "        plasma_time.append(a.index[t])\n",
    "\n",
    "    t+=1\n",
    "start = plasma_time[0]-1e-03 \n",
    "end = plasma_time[-1]-1e-03 \n",
    "print('start =', round(start, 3), 'ms')\n",
    "print('end =', round(end, 3), 'ms')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_cut = r.loc[start:end]\n",
    "a_cut = a.loc[start:end]\n",
    "z_cut = z.loc[start:end]\n",
    "df_processed = pd.concat(\n",
    "    [r_cut.rename('r'), z_cut.rename('z'), a_cut.rename('a')], axis= 'columns')\n",
    "df_processed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "savedata = 'plasma_position.csv'\n",
    "df_processed.to_csv(savedata)\n",
    "\n",
    "Markdown(\"[Plasma position data - $r, z, a$ ](./{})\".format(savedata))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hline = hv.HLine(0)\n",
    "hline.opts(\n",
    "    color='k', \n",
    "    line_dash='dashed',\n",
    "    alpha = 0.4,\n",
    "    line_width=1.0)\n",
    "\n",
    "layout = hv.Layout([df_processed[v].hvplot.line(\n",
    "    xlabel='', ylabel=l,ylim=(-85,85), xlim=(start,end),legend=False, title='', grid=True, group_label=v)\n",
    "                    for (v, l) in [('r', ' r [mm]'), ('z', 'z [mm]'), ('a', 'a [mm]')] ])*hline\n",
    "\n",
    "plot=layout.cols(1).opts(hv.opts.Curve(width=600, height=200),  \n",
    "                    hv.opts.Curve('a', xlabel='time [ms]'))\n",
    "plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "responsestab = requests.get(\"http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/PowSup4InnerStab\"%shot_no)\n",
    "stab = responsestab.text\n",
    "\n",
    "if '-1' in stab or 'Not Found' in stab:\n",
    "    fig, axs = plt.subplots(3, 1, sharex=True, dpi=200)\n",
    "    r.plot(grid=True, ax=axs[0])\n",
    "    z.plot(grid=True, ax=axs[1])\n",
    "    a.plot(grid=True, ax=axs[2])\n",
    "    axs[2].set(ylim=(0,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$a$ [mm]')\n",
    "    axs[1].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\\Delta$z [mm]')\n",
    "    axs[0].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\\Delta$r [mm]', title = 'Horizontal, vertical plasma position and radius #{}'.format(shot_no))\n",
    "\n",
    "else:\n",
    "    fig, axs = plt.subplots(4, 1, sharex=True, dpi=200)\n",
    "    dataNI = (read_signal(shot_no, 'Rog', data_URL, 'lvm')) #data from NI\n",
    "    dataNI = dataNI.set_index('Time')\n",
    "    dataNI = dataNI.set_index(dataNI.index*1e3)\n",
    "    Irog = dataNI['RogQuadr']\n",
    "#     Irog=abs(Irog)\n",
    "    Irog*=1/5e-3\n",
    "    r.plot(grid=True, ax=axs[0])\n",
    "    z.plot(grid=True, ax=axs[1])\n",
    "    a.plot(grid=True, ax=axs[2])\n",
    "    Irog.plot(grid=True, ax=axs[3])\n",
    "    axs[3].set(xlim=(start,end), ylabel = 'I [A]')\n",
    "    axs[2].set(ylim=(0,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$a$ [mm]')\n",
    "    axs[1].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\\Delta$z [mm]')\n",
    "    axs[0].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\\Delta$r [mm]', title = 'Horizontal, vertical plasma position and radius #{}'.format(shot_no))\n",
    "\n",
    "\n",
    "plt.savefig('icon-fig')   "
   ]
  }
 ],
 "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": 2
}