{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plasma Position\n",
    "The mission of this script is to process data from Mirnov coils, determine plasma position and thus help with plasma stabilization. Mirnov coils are type of magnetic diagnostics and are used to the plasma position determination on GOLEM. Four Mirnov coils are placed inside the tokamak 93 mm from the center of the chamber. The effective area of each coil is $A_{eff}=3.8\\cdot10^{-3}$ m. Coils $mc_9$ and $mc_1$ are used to determination the horizontal plasma position and coils $mc_5$ and $mc_{13}$ to determination vertical plasma position.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table>\n",
    "  <tr>\n",
    "    <td valign=\"top\"><img src=\"icon-fig.png\" width=80%></td>\n",
    "    <td valign=\"top\"><img src=\"plasma_position.png\" width=100%></td>\n",
    "  </tr>\n",
    "</table>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Procedure (<a href=20_Position_DK.ipynb>This notebook to download</a>)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "\n",
    "#remove files from the previous discharge\n",
    "def remove_old(name):\n",
    "    if os.path.exists(name):\n",
    "        os.remove(name)\n",
    "        \n",
    "names = ['analysis.html','icon-fig.png', 'plasma_position.png']\n",
    "for name in names:\n",
    "    remove_old(name)\n",
    "    \n",
    "if os.path.exists('Results/'):\n",
    "    file=glob.glob('Results/*')\n",
    "    for f in file:\n",
    "        os.remove(f)\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy import integrate, signal, interpolate\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "import holoviews as hv\n",
    "hv.extension('bokeh')\n",
    "import hvplot.pandas\n",
    "import requests\n",
    "\n",
    "from IPython.display import Markdown\n",
    "\n",
    "#execute time\n",
    "from timeit import default_timer as timer\n",
    "start_execute = timer()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BDdata_URL = \"http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/Basic/{identifier}.csv\"\n",
    "data_URL = \"http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/LimiterMirnovCoils/Default/DAS_raw_data_dir/NIdata_6133.lvm\" #NI Standart (DAS) new address\n",
    "\n",
    "\n",
    "shot_no = 43655 # to be replaced by the actual discharge number\n",
    "\n",
    "vacuum_shot =  # to be replaced by the discharge command line paramater \"vacuum_shot\"\n",
    "\n",
    "\n",
    "ds = np.DataSource(destpath='')  "
   ]
  },
  {
   "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', '13', '14', '15', '16']\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Data integration and $B_t$ elimination\n",
    "\n",
    "The coils measure voltage induced by changes in poloidal magnetic field. In order to obtain measured quantity (poloidal magnetic field), it is necessary to integrate the measured voltage and multiply by constant: $B(t)=\\frac{1}{A_{eff}}\\int_{0}^{t} U_{sig}(\\tau)d\\tau$\n",
    "\n",
    "\n",
    "Ideally axis of the coils is perpendicular to the toroidal magnetic field, but in fact they are slightly deflected and measure toroidal magnetic field too. To determine the position of the plasma, we must eliminate this additional signal. For this we use vacuum discharge with the same parameters."
   ]
  },
  {
   "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",
    "       \n",
    "    if vacuum_shot == False: \n",
    "        signal_start = mc.index[0]\n",
    "        length = len(mc)\n",
    "        Bt = read_signal(shot_no, 'BtCoil_integrated', BDdata_URL).loc[signal_start:signal_start+length*1e-06]\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": [
    "### Plasma life time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    \n",
    "plasma_start = float(requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/DetectPlasma/Results/t_plasma_start').text)\n",
    "plasma_end = float(requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/DetectPlasma/Results/t_plasma_end').text)\n",
    "\n",
    "\n",
    "print ('Plasma start =', round(plasma_start, 3), 'ms')\n",
    "print ('Plasma end =', round(plasma_end, 3), 'ms')\n",
    "# print (CD_orientation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plasma Position\n",
    "To calculate displacement of plasma column, we can use approximation of the straight conductor. If we measure the poloidal field on the two opposite sides of the column, its displacement can be expressed as: $\\Delta=\\frac{B_{\\Theta=0}-B_{\\Theta=\\pi}}{B_{\\Theta=0}+B_{\\Theta=\\pi}}\\cdot b$ "
   ]
  },
  {
   "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",
    "    \n",
    "#     return r.loc[plasma_start:]\n",
    "    return r.loc[plasma_start:plasma_end]\n",
    "#     return r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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]', title = 'Horizontal plasma position #{}'.format(shot_no))\n",
    "ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)\n"
   ]
  },
  {
   "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:]\n",
    "    return z.loc[plasma_start:plasma_end]\n",
    "#     return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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:]\n",
    "    return a.loc[plasma_start:plasma_end]\n",
    "#     return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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]', title = 'Plasma column radius #{}'.format(shot_no))"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "#### Data"
   ]
  },
  {
   "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",
    "\n",
    "df_processed = pd.concat([r_cut.rename('r'), z_cut.rename('z'), a_cut.rename('a')], axis= 'columns')\n",
    "df_processed\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "savedata = 'plasma_position.csv'\n",
    "df_processed.to_csv(savedata)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Markdown(\"[Plasma position data - $r, z, a$ ](./{})\".format(savedata))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hline = hv.HLine(0)\n",
    "hline.opts(color='k',line_dash='dashed', alpha = 0.4,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),  hv.opts.Curve('a', xlabel='time [ms]'))\n",
    "\n",
    "plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_processed.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_end=round(z_cut.iat[-25],2)\n",
    "z_start=round(z_cut.iat[25],2)\n",
    "r_end=round(r_cut.iat[-25],2)\n",
    "r_start=round(r_cut.iat[25],2)\n",
    "\n",
    "name=('z_start', 'z_end', 'r_start', 'r_end','plasma_start', 'plasma_end')\n",
    "values=(z_start, z_end, r_start, r_end, start, end)\n",
    "data={'Values':[z_start, z_end, r_start, r_end, start, end]}\n",
    "results=pd.DataFrame(data, index=name)\n",
    "    \n",
    "dictionary='Results/'\n",
    "j=0\n",
    "for i in values :\n",
    "    with open(dictionary+name[j], 'w') as f:\n",
    "        f.write(str(i))\n",
    "    j+=1\n",
    "    \n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Icon fig"
   ]
  },
  {
   "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 not '-1' in stab or 'Not Found' in stab and exStab==False:\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]', \n",
    "               title = 'Horizontal, vertical plasma position and radius #{}'.format(shot_no))\n",
    "    \n",
    "else:\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]', \n",
    "               title = 'Horizontal, vertical plasma position and radius #{}'.format(shot_no))\n",
    "    \n",
    "plt.savefig('icon-fig')   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plasma position during discharge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tms = np.arange(np.round(start),np.round(end)-1,1)\n",
    "\n",
    "circle1 = plt.Circle((0, 0), 100, color='black', fill = False, lw = 5)\n",
    "circle2 = plt.Circle((0, 0), 85, color='black', fill = False, lw = 3)\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (10,10)) \n",
    "ax.axvline(x=0, color = 'black', alpha = 0.5)\n",
    "ax.axhline(y=0, color = 'black', alpha = 0.5)\n",
    "\n",
    "ax.set_xlim(-100,100)\n",
    "ax.set_ylim(-100,100)\n",
    "ax.add_patch(circle1)\n",
    "ax.add_patch(circle2)\n",
    "plt.annotate('vessel',xy=(27-100, 170-100), xytext=(2-100, 190-100),\n",
    "             arrowprops=dict(facecolor='black', shrink=0.05), fontsize = 12)\n",
    "plt.annotate('limiter',xy=(175-100, 140-100), xytext=(170-100, 190-100),\n",
    "             arrowprops=dict(facecolor='black', shrink=0.05), fontsize = 12)\n",
    "\n",
    "color_idx = np.linspace(0, 1, tms.size)\n",
    "a = 0\n",
    "for i in tms:\n",
    "    mask_pos = ((df_processed.index>i) & (df_processed.index<i+1))\n",
    "    a_mean = df_processed['a'][mask_pos].mean()\n",
    "    r_mean = df_processed['r'][mask_pos].mean()\n",
    "    z_mean = df_processed['z'][mask_pos].mean()   \n",
    "    circle3 = plt.Circle((r_mean, z_mean), a_mean, fill = False, \n",
    "                         color=plt.cm.jet(color_idx[a]),label = str(i) +'-' + str(i+1) + ' ms')\n",
    "    ax.add_patch(circle3)\n",
    "    a+=1\n",
    "    \n",
    "    \n",
    "plt.grid(True)\n",
    "plt.title('Estimated plasma position during discharge', fontweight = 'bold', fontsize = 20)\n",
    "leg=plt.legend(bbox_to_anchor = (1,1))\n",
    "for text in leg.get_texts():\n",
    "     plt.setp(text, fontsize = 16)\n",
    "ax.set_xlabel('r [mm]',fontweight = 'bold', fontsize = 15)\n",
    "ax.set_ylabel('z [mm]', fontweight = 'bold', fontsize = 15)\n",
    "plt.xticks( fontweight = 'bold', fontsize = 14)\n",
    "plt.yticks( fontweight = 'bold', fontsize = 14)\n",
    "\n",
    "plt.savefig('plasma_position.png', bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if not '-1' in stab and not 'Not Found' in stab:\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",
    "    \n",
    "    Irog*=1/5e-3\n",
    "    Irog=Irog.loc[start:end]\n",
    "    data = pd.concat([df_processed, Irog], axis='columns')\n",
    "    data.to_csv('plasma_position.csv')\n",
    "    data['RogQuadr'].plot()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": 2
}