{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e6cb0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle\n",
    "import numpy as np\n",
    "from lmfit.models import LinearModel\n",
    "from ipywidgets import interact, widgets\n",
    "\n",
    "def eckstein(alpha, Y0, alpha0, b, c, f):\n",
    "    return Y0*(np.cos((alpha/alpha0*np.pi/2)**c))**(-1*f)*np.exp(b*(1-1/(np.cos((alpha/alpha0*np.pi/2)**c))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e530d89",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in all the evaluated results from SDTrimSP-3D simulations\n",
    "df = pd.read_csv('Data/Porous_Sputter_Yields.csv')\n",
    "rough = pd.read_csv('Data/Rough_Sputter_Yields.csv')\n",
    "\n",
    "# Eckstein parameters from one-dimensional SDTrimSP simulations\n",
    "eckstein_params = {\n",
    "    500:\n",
    "        {'Ar':[0.62119884, 1.39466856, 1.00532345, 0.73832892, 2.0396859 ],\n",
    "         'D':[1.18845934e-03, 1.46376946e+00, 5.91075310e-01, 8.30694868e-01,1.41439404e+00]},\n",
    "    1000:\n",
    "        {'Ar':[1.03985469, 1.30263466, 0.63416938, 0.67890611, 1.53976874],\n",
    "         'D':[0.0049733 , 1.68528314, 0.343925  , 1.19078647, 1.61273387]},\n",
    "    2000:\n",
    "        {'Ar':[1.52733054, 1.25581977, 0.43887246, 0.65186005, 1.28167341],\n",
    "         'D':[0.00755085, 1.51087844, 0.20662619, 0.92666017, 1.68204546]}\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b81294c",
   "metadata": {},
   "outputs": [],
   "source": [
    "angles = df['Angle'].unique()\n",
    "energies = [500,1000,2000]\n",
    "projectiles = ['Ar','D']\n",
    "filling_continuous = np.linspace(0.1,.9,50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cff646e",
   "metadata": {},
   "outputs": [],
   "source": [
    "angle_ss = widgets.SelectionSlider(\n",
    "    options = angles,\n",
    "    value=45,\n",
    "    description='Angle:',\n",
    "    disabled=False,\n",
    ")\n",
    "energy_dd = widgets.Dropdown(\n",
    "    options = energies,\n",
    "    value=1000,\n",
    "    description='Energy:',\n",
    "    disabled=False,\n",
    ")\n",
    "proj_dd = widgets.Dropdown(\n",
    "    options = projectiles,\n",
    "    value='Ar',\n",
    "    description='Projectile:',\n",
    "    disabled=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5239ab5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(angle=angle_ss, energy=energy_dd, projectile=proj_dd)\n",
    "\n",
    "def plot(angle, energy, projectile):\n",
    "    \n",
    "    # Data wrangling\n",
    "    filtered_df = df[(df.Angle == angle) & \n",
    "                     (df.Energy_eV == energy) & \n",
    "                     (df.Atom == projectile) ]\n",
    "    flat = eckstein(np.deg2rad(angle), *eckstein_params[energy][projectile])\n",
    "    df_r = rough[(rough.Angle == angle) & \n",
    "                 (rough.Energy_eV == energy) & \n",
    "                 (rough.Atom == projectile) ]\n",
    "\n",
    "    # Plot sputter yields of porous structures \n",
    "    fig, ax = plt.subplots(1,1,layout='constrained')\n",
    "    ax.errorbar(filtered_df['Filling'], filtered_df['Sputter_Coeff_Mean'], \n",
    "                yerr=filtered_df['Sputter_Coeff_Std'],\n",
    "                ls='none', marker='o', color='tab:blue', markerfacecolor='w')\n",
    "    \n",
    "    # Calculate and plot linear fit through data\n",
    "    model = LinearModel()\n",
    "    params = model.make_params()\n",
    "    fit_results = model.fit(filtered_df['Sputter_Coeff_Mean'], \n",
    "                            params=params, x=filtered_df['Filling'])\n",
    "    ax.plot(filling_continuous, fit_results.eval(x=filling_continuous), \n",
    "            color='tab:blue', ls='-')\n",
    "    \n",
    "    # Plot extrapolation towards solid target\n",
    "    ax.plot(np.linspace(.9,1,50), fit_results.eval(x=np.linspace(.9,1,50)), \n",
    "            color='tab:blue', ls=':')\n",
    "\n",
    "    # Annotate/write fit parameters\n",
    "    ax.text(.1, 0.1, f'fit: $Y = m \\\\times \\\\rho_\\\\mathrm{{fil}} + b$\\n$m={fit_results.best_values[\"slope\"]:.4f}$\\n$b={fit_results.best_values[\"intercept\"]:.4f}$', \n",
    "            transform=ax.transAxes, color='tab:blue')\n",
    "    \n",
    "    # Plot flat surface, rough surface yield\n",
    "    ax.scatter([1], flat, color='tab:orange', zorder=5.)\n",
    "    ax.scatter(df_r['Filling'], df_r['Sputter_Coeff_Mean'], color='tab:red', zorder=6)\n",
    "    \n",
    "    # Annotate solid target, flat & rough yield\n",
    "    ax.axvline(1, color='grey', ls='--')\n",
    "    ax.text(1.02, df_r['Sputter_Coeff_Mean'].iloc[0], 'rough\\n$Y_\\\\mathrm{r}$', \n",
    "            color='tab:red')\n",
    "    ax.text(.9, flat, 'flat\\n$Y_\\\\mathrm{f}$', color='tab:orange')\n",
    "    ax.text(.95, df_r['Sputter_Coeff_Mean'].iloc[0]/5, 'solid target', \n",
    "            rotation='vertical', color='gray')\n",
    "\n",
    "    # Axis labels and limits\n",
    "    ax.set_xlabel('Volume filling factor $\\\\rho_\\\\mathrm{fil}$')\n",
    "    ax.set_ylabel(f'Sputter yield $Y$ (W/{projectile})')\n",
    "    ax.set_xlim(0,1.2)\n",
    "    ax.set_xticks([.2,.4,.6,.8,1])\n",
    "    ax.set_ylim(0,None)\n",
    "\n",
    "    # Title\n",
    "    ax.set_title(f'{energy}$\\,$eV {projectile} $\\\\rightarrow$ W @ {angle}°')\n",
    "    \n",
    "    # Second axis, normalised to flat surface yields\n",
    "    ax2 = ax.twinx()\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax2.spines['top'].set_visible(False)\n",
    "    ax2.spines['right'].set_color('grey')\n",
    "    (y_low, y_high) = ax.get_ylim()\n",
    "    ax2.set_ylim(y_low/flat, y_high/flat)\n",
    "    ax2.set_ylabel('Normalised sputter yield $Y/Y_\\\\mathrm{f}$', color='grey')\n",
    "    ax2.tick_params(colors='grey')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe83b135",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}