{
 "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. It also calculates the plasma column position from the fast camera data and compares the two diagnostics.\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('InnerStabilization.html'):\n",
    "    os.remove('InnerStabilization.html')\n",
    "    \n",
    "if os.path.exists('video.mp4'):\n",
    "    os.remove('video.mp4')\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",
    "\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()"
   ]
  },
  {
   "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",
    "shot_no = 35362 # to be replaced by the actual discharge number\n",
    "#shot_no = 35021\n",
    "vacuum_shot = 35353 # to be replaced by the discharge command line paramater \"vacuum_shot\"\n",
    "#vacuum_shot = 35022 #number of the vacuum shot or 'False'\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",
    "\n",
    "Find plasma start and end from the basic diagnostics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loop_voltage = read_signal(shot_no, 'U_Loop', BDdata_URL)\n",
    "\n",
    "dIpch = read_signal(shot_no, 'U_RogCoil', BDdata_URL)\n",
    "\n",
    "dIpch -= dIpch.loc[:0.9e-3].mean()\n",
    "\n",
    "Ipch = pd.Series(integrate.cumtrapz(dIpch, x=dIpch.index, initial=0) * (-5.3*1e06),\n",
    "                index=dIpch.index, name='Ipch')\n",
    "\n",
    "U_l_func = interpolate.interp1d(loop_voltage.index, loop_voltage)  \n",
    "def dIch_dt(t, Ich):\n",
    "    return (U_l_func(t) - 0.0097 * Ich) / (1.2e-6/2)\n",
    "t_span = loop_voltage.index[[0, -1]]\n",
    "solution = integrate.solve_ivp(dIch_dt, t_span, [0], t_eval=loop_voltage.index, )\n",
    "Ich = pd.Series(solution.y[0], index=loop_voltage.index, name='Ich')\n",
    "Ip = Ipch - Ich\n",
    "Ip.name = 'Ip'\n",
    "\n",
    "Ip_detect = Ip.loc[0.0025:]\n",
    "\n",
    "dt = (Ip_detect.index[-1] - Ip_detect.index[0]) / (Ip_detect.index.size) \n",
    "\n",
    "window_length = int(0.5e-3/dt)  \n",
    "if window_length % 2 == 0:  \n",
    "    window_length += 1\n",
    "dIp = pd.Series(signal.savgol_filter(Ip_detect, window_length, polyorder=3, deriv=1, delta=dt),\n",
    "                name='dIp', index=Ip_detect.index) / 1e6 \n",
    "\n",
    "threshold = 0.05\n",
    "\n",
    "CD = requests.get(\"http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/CD_orientation\" % shot_no)\n",
    "CD_orientation = CD.text\n",
    "\n",
    "if \"ACW\" in CD_orientation:\n",
    "    plasma_start = dIp[dIp < dIp.min()*threshold].index[0]*1e3 \n",
    "    plasma_end = dIp.idxmax()*1e3 \n",
    "else: \n",
    "    plasma_start = dIp[dIp > dIp.max()*threshold].index[0]*1e3 \n",
    "    plasma_end = dIp.idxmin()*1e3 + 0.5    \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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mc1 = elimination(shot_no, 'mc1', vacuum_shot)"
   ]
  },
  {
   "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]',\n",
    "       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": "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\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "\n",
    "plt.savefig('icon-fig')   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fast cameras\n",
    "\n",
    "The plasma vertical position can also be calculated from the fast camera signal as the centre of the radiating zone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "from io import BytesIO\n",
    "from IPython.display import Image as DispImage\n",
    "\n",
    "FastCamera = True\n",
    "\n",
    "\n",
    "try:\n",
    "    url = \"http://golem.fjfi.cvut.cz/shots/%i/Diagnostics/FastExCameras/plasma_film.png\"%(shot_no)\n",
    "    image = requests.get(url)\n",
    "    img = Image.open(BytesIO(image.content)).convert('P') #'L','1'\n",
    "\n",
    "    data = np.asanyarray(img)\n",
    "\n",
    "    r = []\n",
    "    for i in range(data.shape[1]):\n",
    "        a=0\n",
    "        b=0\n",
    "        for j in range(336):\n",
    "            a += data[j,i]*j\n",
    "            b += data[j,i]\n",
    "        r.append((a/b)*(170/336))\n",
    "    Tcd = 3 #delay\n",
    "    camera_time = np.linspace(start+abs(start-Tcd), end, len(r))\n",
    "    r_camera = pd.Series(r, index = camera_time)-85\n",
    "    r_mirnov = horpol(shot_no, vacuum_shot).loc[start:end]\n",
    "#     print(start, end, len(r))\n",
    "        \n",
    "    fig = plt.figure()\n",
    "    ax = r_camera.plot(label = 'Fast Camera')\n",
    "    ax = r_mirnov.plot(label = 'Mirnov coils')\n",
    "    plt.legend()\n",
    "    ax.set(ylim=(-85, 85), title='#%i'%shot_no,ylabel='$\\Delta$r visible radiation')\n",
    "    ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)\n",
    "    fig.savefig('FastCameras_deltaR')\n",
    "except OSError:\n",
    "    FastCamera = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if FastCamera == True:\n",
    "    plt.imshow(img)"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}