{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0b595594",
   "metadata": {},
   "source": [
    "# Imports, Definitions, Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e6cb0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.gridspec import GridSpec\n",
    "import numpy as np\n",
    "from lmfit.models import LinearModel\n",
    "from ipywidgets import interact, widgets\n",
    "\n",
    "blue, orange, green, red, purple, pink =  '#49C1DF', '#F08000', '#00dd6a', '#E21A45', '#A982FF', '#EB259B'\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))))\n",
    "\n",
    "plt.style.use('./single_column.mplstyle')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "301da4a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ar_recoils = pd.read_csv('Data/Polar_angles_of_backsputtered_recoils_Ar.txt', skiprows=1, sep='\\s+', \n",
    "                        names=['angle', 'freq'])\n",
    "d_recoils = pd.read_csv('Data/Polar_angles_of_backsputtered_recoils_D.txt', skiprows=1, sep='\\s+', \n",
    "                        names=['angle', 'freq'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e530d89",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in all the evaluated results from SDTrimSP-3D simulations\n",
    "df_porous = pd.read_csv('Data/Porous_Sputter_Yields.csv')\n",
    "df_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_porous['Angle'].unique()\n",
    "energies = [500,1000,2000]\n",
    "projectiles = ['Ar','D']\n",
    "filling_continuous = np.linspace(0.1,.9,50)\n",
    "angles_sim = np.linspace(0,85,50)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1bbfb89",
   "metadata": {},
   "source": [
    "# Figure 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b706dd43",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up Figure\n",
    "fig, (ax, ax2) = plt.subplots(2, 1, figsize=(3.39, 3.5), sharex=True, layout='constrained')\n",
    "\n",
    "# Plot and annotate flat data line\n",
    "ax.plot(angles_sim, eckstein(np.deg2rad(angles_sim), *eckstein_params[1000]['Ar']), color=orange)\n",
    "ax.text(20, 1.25, 'flat', color=orange, rotation=20)\n",
    "\n",
    "# Plot exemplary fuzz data\n",
    "df_porous[(df_porous['Energy_eV'] == 1000) & \n",
    "          (df_porous['Atom']=='Ar') & \n",
    "          (df_porous['Filling'] == 0.2)].plot(x='Angle', y='Sputter_Coeff_Mean', \n",
    "                                              ax=ax, label=0.2, color=green)\n",
    "df_porous[(df_porous['Energy_eV'] == 1000) & \n",
    "          (df_porous['Atom']=='Ar') & \n",
    "          (df_porous['Filling'] == 0.5)].plot(x='Angle', y='Sputter_Coeff_Mean', \n",
    "                                              ax=ax, label=0.5, color=red)\n",
    "df_porous[(df_porous['Energy_eV'] == 1000) & \n",
    "          (df_porous['Atom']=='Ar') & \n",
    "          (df_porous['Filling'] == 0.8)].plot(x='Angle', y='Sputter_Coeff_Mean', \n",
    "                                              ax=ax, label=0.8, color=blue)\n",
    "\n",
    "# Draw the arrows and text annotating the fuzz lines\n",
    "con_style = \"angle, angleA=180, angleB=-90\"\n",
    "ax.annotate(\"\",\n",
    "            xy=(40,.88), xycoords='data',\n",
    "            xytext=(50,1.08), textcoords='data',\n",
    "            arrowprops=dict(arrowstyle=\"-|>\", color=blue,\n",
    "                            connectionstyle=con_style))\n",
    "ax.annotate(\"\",\n",
    "            xy=(42,.72), xycoords='data',\n",
    "            xytext=(50,1.05), textcoords='data',\n",
    "            arrowprops=dict(arrowstyle=\"-|>\", color=red,\n",
    "                            connectionstyle=con_style))\n",
    "ax.annotate(\"\",\n",
    "            xy=(44,.59), xycoords='data',\n",
    "            xytext=(50,1.02), textcoords='data',\n",
    "            arrowprops=dict(arrowstyle=\"-|>\", color=green,\n",
    "                            connectionstyle=con_style))\n",
    "ax.text(50, 1.02, 'Fuzz')\n",
    "\n",
    "# Labels and limits\n",
    "ax.set_ylim(0)\n",
    "#ax.set_xlabel('Incidence angle $\\\\alpha$ (°)')\n",
    "\n",
    "# Legend and panel title\n",
    "ax.legend(title='Volume filling $\\\\rho_\\mathrm{fil}$', ncol=3)\n",
    "ax.text(.05, .875, f'a) $1\\,$keV Ar $\\\\rightarrow$ W', transform=ax.transAxes)\n",
    "ax.get_legend().remove()\n",
    "\n",
    "# Repeat for D in second axis\n",
    "energy_D = 1000\n",
    "ax2.plot(angles_sim, eckstein(np.deg2rad(angles_sim), *eckstein_params[energy_D]['D']), color=orange)\n",
    "\n",
    "df_porous[(df_porous['Energy_eV'] == energy_D) & \n",
    "          (df_porous['Atom']=='D') & \n",
    "          (df_porous['Filling'] == 0.2)].plot(x='Angle', y='Sputter_Coeff_Mean', \n",
    "                                              ax=ax2, label=0.2, color=green)\n",
    "df_porous[(df_porous['Energy_eV'] == energy_D) & \n",
    "          (df_porous['Atom']=='D') & \n",
    "          (df_porous['Filling'] == 0.5)].plot(x='Angle', y='Sputter_Coeff_Mean', \n",
    "                                              ax=ax2, label=0.5, color=red)\n",
    "df_porous[(df_porous['Energy_eV'] == energy_D) & \n",
    "          (df_porous['Atom']=='D') & \n",
    "          (df_porous['Filling'] == 0.8)].plot(x='Angle', y='Sputter_Coeff_Mean', \n",
    "                                              ax=ax2, label=0.8, color=blue)\n",
    "\n",
    "ax2.legend(title='Volume filling $\\\\rho_\\mathrm{fil}$', ncol=3, loc='lower right', frameon=False)\n",
    "ax2.set_ylim(0, None)\n",
    "ax2.set_xlabel('Incidence angle $\\\\alpha$ (°)')\n",
    "ax2.ticklabel_format(axis='y', style='sci', scilimits=(0,0), useMathText=True)\n",
    "ax2.text(.05,.875, f'b) $1\\,$keV D $\\\\rightarrow$ W', transform=ax2.transAxes)\n",
    "\n",
    "# Create and center supylabel\n",
    "bbox1 = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())\n",
    "bbox2 = ax2.get_window_extent().transformed(fig.dpi_scale_trans.inverted())\n",
    "fig.supylabel('Sputter yield $Y$ (W/ion)', y=(bbox2.y0+bbox1.y1)/2/bbox1.y1, va='center')\n",
    "\n",
    "plt.savefig('Figures/exemplary_results_1_keV_Ar.pdf')\n",
    "fig.savefig('Figures/Figure_3.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea41d9ce",
   "metadata": {},
   "source": [
    "# Figure 4\n",
    "\n",
    "Plot the example cases for bombardment under 45° for 1 keV Ar in the left panel and 0.5 keV D in the right panel.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c5cea6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up plot\n",
    "fig, (ax1, ax3) = plt.subplots(1,2,layout='constrained')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "794e13fe",
   "metadata": {},
   "source": [
    "## Figure 4a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1040bf9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define point in parameter space\n",
    "angle = 45\n",
    "projectile = 'Ar'\n",
    "energy = 1000\n",
    "\n",
    "# Data wrangling\n",
    "filtered_df = df_porous[(df_porous.Angle == angle) & \n",
    "                        (df_porous.Energy_eV == energy) & \n",
    "                        (df_porous.Atom == projectile) ]\n",
    "flat = eckstein(np.deg2rad(angle), *eckstein_params[energy][projectile])\n",
    "df_r = df_rough[(df_rough.Angle == angle) & \n",
    "                (df_rough.Energy_eV == energy) & \n",
    "                (df_rough.Atom == projectile) ]\n",
    "\n",
    "# Prepare linear fit\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",
    "r_squ = fit_results.summary()['rsquared']\n",
    "\n",
    "# Plot porous data\n",
    "ax1.errorbar(filtered_df['Filling'], filtered_df['Sputter_Coeff_Mean'], \n",
    "             yerr=filtered_df['Sputter_Coeff_Std'],\n",
    "            ls='none', marker='o', markerfacecolor='w')\n",
    "\n",
    "# Plot linear fit\n",
    "ax1.plot(filling_continuous, fit_results.eval(x=filling_continuous), color=blue, ls='-')\n",
    "# Plot extrapolation of linear fit\n",
    "ax1.plot(np.linspace(.9,1,50), fit_results.eval(x=np.linspace(.9,1,50)), color=blue, ls=':')\n",
    "\n",
    "# Plot data for solid target (flat and rough)\n",
    "ax1.scatter([1], flat, color=orange, zorder=5)\n",
    "ax1.scatter(df_r['Filling'], df_r['Sputter_Coeff_Mean'], color=red, zorder=6)\n",
    "\n",
    "# Annotate region for solid target\n",
    "ax1.axvline(1, color='grey', ls='--')\n",
    "ax1.text(1.01, .47, 'solid', rotation='vertical', color='gray')\n",
    "\n",
    "# Draw annotation and arrow for flat and rough data points\n",
    "ax1.text(1.03, df_r['Sputter_Coeff_Mean'].iloc[0]*1.16, 'rough $Y_\\\\mathrm{r}$', \n",
    "         color=red, rotation=90)\n",
    "ax1.annotate('', xy=(1,flat), xytext=(.85, flat),\n",
    "             arrowprops=dict(arrowstyle='-|>', \n",
    "                             color=orange))\n",
    "ax1.text(.7, flat*.95, 'flat $Y_\\\\mathrm{f}$', color=orange)\n",
    "ax1.annotate('', xy=(1,df_r['Sputter_Coeff_Mean'].iloc[0]), xytext=(1.07, 1.2),\n",
    "             arrowprops=dict(arrowstyle='-|>', \n",
    "                             connectionstyle='angle, angleA=90',\n",
    "                             color=red))\n",
    "\n",
    "# Labels, ticks and limits\n",
    "ax1.set_ylabel(f'Sputter yield $Y$ (W/ion)')\n",
    "ax1.set_xlim(0,1.13)\n",
    "ax1.set_xticks([.2,.4,.6,.8,1])\n",
    "ax1.set_ylim(None,1.7)\n",
    "\n",
    "# Text annotations for the fit line\n",
    "ax1.text(.1, 0.62, f'$Y = m \\\\times \\\\rho_\\\\mathrm{{fil}} + b$', color=blue, rotation=17)\n",
    "ax1.text(.05, .64, \n",
    "         f'$m={fit_results.best_values[\"slope\"]:.4f}$\\n$b={fit_results.best_values[\"intercept\"]:.4f}$\\n$R^2={r_squ:.4f}$', \n",
    "        transform=ax1.transAxes, color=blue)\n",
    "\n",
    "# Create twin axis with coordinates relative to flat yield\n",
    "ax2 = ax1.twinx()\n",
    "ax2.spines['top'].set_visible(True)\n",
    "ax2.spines['right'].set_color('grey')\n",
    "(y_low, y_high) = ax1.get_ylim()\n",
    "ax2.set_ylim(y_low/flat, y_high/flat)\n",
    "ax2.tick_params(colors='grey')\n",
    "\n",
    "# Panel caption\n",
    "ax1.text(.05, .9, \n",
    "         f'a) ${int(energy/1000)}\\,$keV {projectile} $\\\\rightarrow$ W @ {angle}°', \n",
    "         transform=ax1.transAxes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "947c02f9",
   "metadata": {},
   "source": [
    "## Figure 4b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb7a5738",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Define point in parameter space\n",
    "angle = 45\n",
    "projectile = 'D'\n",
    "energy = 500\n",
    "\n",
    "# Data wrangling\n",
    "filtered_df = df_porous[(df_porous.Angle == angle) & \n",
    "                        (df_porous.Energy_eV == energy) & \n",
    "                        (df_porous.Atom == projectile) ]\n",
    "flat = eckstein(np.deg2rad(angle), *eckstein_params[energy][projectile])\n",
    "df_r = df_rough[(df_rough.Angle == angle) & \n",
    "                (df_rough.Energy_eV == energy) & \n",
    "                (df_rough.Atom == projectile) ]\n",
    "\n",
    "# Prepare linear fit\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",
    "r_squ = fit_results.summary()['rsquared']\n",
    "\n",
    "# Plot porous data\n",
    "ax3.errorbar(filtered_df['Filling'], filtered_df['Sputter_Coeff_Mean'], \n",
    "             yerr=filtered_df['Sputter_Coeff_Std'],\n",
    "            ls='none', marker='o', markerfacecolor='w')\n",
    "\n",
    "# Plot linear fit\n",
    "ax3.plot(filling_continuous, fit_results.eval(x=filling_continuous), color=blue, ls='-')\n",
    "# Plot extrapolation of linear fit\n",
    "ax3.plot(np.linspace(.9,1,50), fit_results.eval(x=np.linspace(.9,1,50)), color=blue, ls=':')\n",
    "\n",
    "# Plot data for solid target (flat and rough)\n",
    "ax3.scatter([1], flat, color=orange, zorder=5)\n",
    "ax3.scatter(df_r['Filling'], df_r['Sputter_Coeff_Mean'], color=red, zorder=6)\n",
    "\n",
    "# Annotate region for solid target\n",
    "ax3.axvline(1, color='grey', ls='--')\n",
    "\n",
    "# Annotate flat and rough data\n",
    "ax3.text(1.03, df_r['Sputter_Coeff_Mean'].iloc[0]*.95, '$Y_\\\\mathrm{r}$', color=red)\n",
    "ax3.text(.92, flat*.92, '$Y_\\\\mathrm{f}$', color=orange)\n",
    "\n",
    "# Labels, ticks and limits\n",
    "ax3.set_xlim(0,1.13)\n",
    "ax3.set_xticks([.2,.4,.6,.8,1])\n",
    "ax3.set_ylim(None,None)\n",
    "ax3.ticklabel_format(axis='y', style='sci', scilimits=(0,0), useMathText=True)\n",
    "\n",
    "# Text annotations for the fit line\n",
    "ax3.text(0.05, .64,  \n",
    "         f'$m={fit_results.best_values[\"slope\"]:.4f}$\\n$b={fit_results.best_values[\"intercept\"]:.4f}$\\n$R^2={r_squ:.4f}$', \n",
    "         transform=ax3.transAxes, color=blue)\n",
    "\n",
    "# Create twin axis with coordinates relative to flat yield\n",
    "ax4 = ax3.twinx()\n",
    "ax4.spines['top'].set_visible(True)\n",
    "ax4.spines['right'].set_color('grey')\n",
    "(y_low, y_high) = ax3.get_ylim()\n",
    "ax4.set_ylim(y_low/flat, y_high/flat)\n",
    "ax4.set_ylabel('Normalised sputter yield $Y/Y_\\\\mathrm{f}$', color='grey')\n",
    "ax4.tick_params(colors='grey')\n",
    "\n",
    "# Panel caption\n",
    "ax3.text(.05, .9, f'b) ${energy/1000}\\,$keV {projectile} $\\\\rightarrow$ W @ {angle}°', \n",
    "         transform=ax3.transAxes)\n",
    "\n",
    "# Global x label for the entire figure\n",
    "fig.supxlabel('Volume filling factor $\\\\rho_\\\\mathrm{fil}$')\n",
    "\n",
    "# Set width to two-column layout of most journals\n",
    "fig.set_figwidth(182/25.4)\n",
    "fig.set_figheight(2.4)\n",
    "\n",
    "# Save\n",
    "fig.savefig('Figures/exemplary_results_over_porosity.pdf')\n",
    "fig.savefig('Figures/Figure_4.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc5a9cbb",
   "metadata": {},
   "source": [
    "# Figure 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca101520",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_filled = pd.read_csv('Data/Filled_Surface_Sputter_Yields.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6123a10e",
   "metadata": {},
   "outputs": [],
   "source": [
    "angle = 45\n",
    "projectile = 'D'\n",
    "energy = 1000\n",
    "\n",
    "fig, ax1 = plt.subplots(1,1,constrained_layout=True)\n",
    "\n",
    "# Data wrangling\n",
    "filtered_df = df_porous[(df_porous.Angle == angle) & \n",
    "                        (df_porous.Energy_eV == energy) & \n",
    "                        (df_porous.Atom == projectile) ]\n",
    "flat = eckstein(np.deg2rad(angle), *eckstein_params[energy][projectile])\n",
    "df_r = df_rough[(df_rough.Angle == angle) & \n",
    "                (df_rough.Energy_eV == energy) & \n",
    "                (df_rough.Atom == projectile) ]\n",
    "filtered_filled = df_filled[(df_filled.Angle == angle) & \n",
    "                        (df_filled.Energy_eV == energy) & \n",
    "                        (df_filled.Atom == projectile) ]\n",
    "\n",
    "# Plot porous data\n",
    "ax1.errorbar(filtered_df['Filling'], filtered_df['Sputter_Coeff_Mean'], \n",
    "             yerr=filtered_df['Sputter_Coeff_Std'],\n",
    "            ls='none', marker='o', markerfacecolor='w',\n",
    "            label='Porous', capsize=2)\n",
    "ax1.errorbar(filtered_filled['Filling'], filtered_filled['Sputter_Coeff_Mean'], \n",
    "             yerr=filtered_filled['Sputter_Coeff_Std'],\n",
    "            ls='none', marker='o', markerfacecolor='w',\n",
    "            label='Filled surface', zorder=9, capsize=2)\n",
    "\n",
    "ax1.legend()\n",
    "\n",
    "# Labels, ticks and limits\n",
    "ax1.set_ylabel(f'Sputter yield $Y$ (W/ion)')\n",
    "ax1.set_xlabel('Volume filling factor $\\\\rho_\\\\mathrm{fil}$')\n",
    "ax1.set_xlim(None,.85)\n",
    "ax1.text(.05, .9, f'$1\\,$keV {projectile} $\\\\rightarrow$ W @ {angle}°', transform=ax1.transAxes)\n",
    "ax1.legend(loc='lower right', frameon=False, labelspacing=.8)\n",
    "ax1.ticklabel_format(axis='y', style='sci', scilimits=(0,0), useMathText=True)\n",
    "\n",
    "fig.savefig('Figures/rough_vs_porous')\n",
    "fig.savefig('Figures/Figure_5.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f90be58",
   "metadata": {},
   "source": [
    "# Figure 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2075617",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract slopes from ALL linear fits, not just example cases above\n",
    "\n",
    "slopes = []\n",
    "\n",
    "for energy in energies:\n",
    "    for angle in angles:\n",
    "        for proj in projectiles:\n",
    "            # Wrangle data\n",
    "            filtered_df = df_porous[(df_porous.Angle == angle) & \n",
    "                                    (df_porous.Energy_eV == energy) & \n",
    "                                    (df_porous.Atom == proj) ]\n",
    "            \n",
    "            # Calculate fit and extract the parameters\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",
    "            slopes.append(dict(energy=energy, \n",
    "                               angle=angle, \n",
    "                               projectile=proj, \n",
    "                               slope=fit_results.best_values['slope'], \n",
    "                               intercept=fit_results.best_values['intercept']))\n",
    "\n",
    "slopes = pd.DataFrame(slopes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3d5426e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set styles for later\n",
    "colors={500:blue, 1000:orange, 2000:green}\n",
    "ls = {'Ar':'--', 'D':'-'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1684afbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up figure \n",
    "fig, ax = plt.subplots(1,1,layout='constrained', figsize=(3.39, 2.4))\n",
    "\n",
    "norm_slopes = {500:{}, 1000:{}, 2000:{}}\n",
    "\n",
    "# Go through parameter space, calculate and plot slopes normalised to rough surface yield\n",
    "for energy in energies:\n",
    "    for proj in projectiles:\n",
    "        norm_slopes_series = slopes[(slopes['projectile'] == proj) & \n",
    "                                    (slopes['energy'] == energy)]['slope'] \\\n",
    "                / np.asarray(df_rough[(df_rough.Energy_eV == energy) & \n",
    "                           (df_rough.Atom == proj) ]['Sputter_Coeff_Mean'])\n",
    "        \n",
    "        norm_slopes[energy][proj] = pd.DataFrame({'angle':angles, \n",
    "                                                  'normalised_slope':norm_slopes_series})\n",
    "        \n",
    "        ax.plot(angles, norm_slopes_series, \n",
    "                label=f'{energy/1000:.1g}', color=colors[energy], ls=ls[proj])\n",
    "\n",
    "# Labels, limits\n",
    "ax.set_xlabel('Incidence angle $\\\\alpha$ (°)')\n",
    "ax.set_ylabel('Normalised slope $m / Y_\\\\mathrm{r}$')\n",
    "ax.set_ylim(0, None)\n",
    "\n",
    "# Do legends manually:\n",
    "# First, use line colours to denote energy\n",
    "lines = ax.get_lines()\n",
    "legend1 = plt.legend([lines[i] for i in [1,3,5]], [lines[i].get_label() for i in [1,3,5]],\n",
    "                        title='Energy (keV)',\n",
    "                        loc='upper left', bbox_to_anchor=(.05, .37),\n",
    "                        frameon=False, handletextpad=.8)\n",
    "for legobj in legend1.legend_handles:\n",
    "    legobj.set_linewidth(7.0)\n",
    "\n",
    "# Second legend: create new axis object that does not plot anything.\n",
    "# Just use it to display another legend, using line styles to denote projectile species\n",
    "ax2 = ax.twinx()\n",
    "ax2.plot(np.NaN, np.NaN, ls='-', color='k', label='Ar')\n",
    "ax2.plot(np.NaN, np.NaN, ls='--', color='k', label='D')\n",
    "ax2.get_yaxis().set_visible(False)\n",
    "ax2.legend(title='Projectile', loc='upper left', bbox_to_anchor=(.33, .37), \n",
    "           frameon=False)\n",
    "\n",
    "# Calculate the average and mean for data <45°\n",
    "for_average = pd.DataFrame()\n",
    "for energy in energies:\n",
    "    for proj in projectiles:\n",
    "        for_average = pd.concat([for_average, \n",
    "                                 norm_slopes[energy][proj]\n",
    "                                 [norm_slopes[energy][proj]['angle'] <= 45]['normalised_slope']])\n",
    "mean, std = for_average.mean(), 2*for_average.std()\n",
    "slope_low = mean-std\n",
    "slope_high = mean+std\n",
    "\n",
    "# Plot shaded band \n",
    "ax.fill_between([0,15,30,45,85], mean-std, mean+std, \n",
    "                color='grey', zorder=0, alpha=.15, linewidth=0)\n",
    "\n",
    "# Preparations to annotate shaded band\n",
    "low = (mean-std)*1.03\n",
    "high =( mean+std)*.97\n",
    "length = .06\n",
    "x_pos = 20\n",
    "\n",
    "# Annotate shaded band\n",
    "ax.annotate(f'', xy=(x_pos,low), xytext=(x_pos, low-length-.03),\n",
    "            arrowprops=dict(arrowstyle='-|>', color='k',\n",
    "                            linewidth=.75))\n",
    "ax.annotate(f'', xy=(x_pos,high), xytext=(x_pos, high+length),\n",
    "            arrowprops=dict(arrowstyle='-|>', color='k',\n",
    "                            linewidth=.75))\n",
    "ax.text(x_pos,.43,f'${mean.iloc[0]:.2f}\\pm{std.iloc[0]:.2f}$', \n",
    "        ha='center', va='top')\n",
    "ax.hlines(mean.iloc[0], 0,85,color='k', zorder=0, alpha=.6)\n",
    "\n",
    "ax.set_ylim(0, .65)\n",
    "\n",
    "# Save\n",
    "fig.savefig('Figures/norm_slope_over_incidence_angle.pdf')\n",
    "fig.savefig('Figures/Figure_7.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4239f3dd",
   "metadata": {},
   "source": [
    "# Figure 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "240a68a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the sputter yield data\n",
    "# These data were either provided by the authors or digitsed from the publication\n",
    "lop_exp_1keV = pd.read_csv('Data/Lopez_Cupak_Exp_1keV_Ar.txt', sep='\\s+')\n",
    "lop_afm_1keV = pd.read_csv('Data/Lopez_Cupak_AFM_1keV_Ar.txt',\n",
    "                           names=['angle', 'yield'])\n",
    "lop_exp_2keV = pd.read_csv('Data/Lopez_Cupak_Exp_2keV_Ar.txt', sep='\\s+')\n",
    "lop_afm_2keV = pd.read_csv('Data/Lopez_Cupak_AFM_2keV_Ar.txt', \n",
    "                           names=['angle', 'yield'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbde693a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This work's linear model prediction for 1 keV\n",
    "m = mean.iloc[0] * lop_afm_1keV['yield']\n",
    "b = lop_afm_1keV['yield'] - m\n",
    "prediction_1keV = pd.DataFrame()\n",
    "prediction_1keV['yield'] = m*.58+b\n",
    "prediction_1keV['angle'] = lop_afm_1keV['angle']\n",
    "\n",
    "# This work's linear model prediction for 2 keV\n",
    "m = mean.iloc[0] * lop_afm_2keV['yield']\n",
    "b = lop_afm_2keV['yield'] - m\n",
    "prediction_2keV = pd.DataFrame()\n",
    "prediction_2keV['yield'] = m * .58 + b\n",
    "prediction_2keV['angle'] = lop_afm_2keV['angle']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7228a40",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the point in parameter space \n",
    "angle = 45\n",
    "energy = 1000\n",
    "projectile = 'Ar'\n",
    "\n",
    "# Wrangle data accordingly\n",
    "filtered_df = df_porous[(df_porous.Angle == angle) & \n",
    "                        (df_porous.Energy_eV == energy) & \n",
    "                        (df_porous.Atom == projectile) ]\n",
    "flat = eckstein(np.deg2rad(angle), *eckstein_params[energy][projectile])\n",
    "df_r = df_rough[(df_rough.Angle == angle) & \n",
    "                (df_rough.Energy_eV == energy) & \n",
    "                (df_rough.Atom == projectile) ]\n",
    "\n",
    "# Formulate prediction according to linear model\n",
    "lop_afm_example = lop_afm_1keV[lop_afm_1keV['angle']==45]['yield'].iloc[0]\n",
    "lop_exp_example = lop_exp_1keV[lop_exp_1keV['angle']==45]\n",
    "exp_slope = mean.iloc[0]*lop_afm_example\n",
    "exp_intercept = lop_afm_example - exp_slope\n",
    "rho = np.linspace(0.1,1,50)\n",
    "\n",
    "# Create two subfigures: \n",
    "# One (top) for the first panel spanning the whole column,\n",
    "# the second one (bottom) will include two subplots aligned using GridSpec\n",
    "fig = plt.figure(layout='constrained', figsize=(3.39, 3.6))\n",
    "(subfig_t, subfig_b) = fig.subfigures(2, 1, hspace=0.05, height_ratios=[1.2, 1])\n",
    "gs = subfig_b.add_gridspec(nrows=1, ncols=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "246258ae",
   "metadata": {},
   "source": [
    "## Figure 8a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea12cdf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create axis\n",
    "ax=subfig_t.subplots()\n",
    "\n",
    "# Plot fuzz data from this study\n",
    "# Individual points\n",
    "ax.errorbar(filtered_df['Filling'], filtered_df['Sputter_Coeff_Mean'], \n",
    "            yerr=filtered_df['Sputter_Coeff_Std'],\n",
    "            ls='none', marker='o', markerfacecolor='w', capsize=2)\n",
    "# Linear model\n",
    "model = LinearModel()\n",
    "params = model.make_params()\n",
    "fit_results = model.fit(filtered_df['Sputter_Coeff_Mean'], params=params, x=filtered_df['Filling'])\n",
    "ax.plot(filling_continuous, fit_results.eval(x=filling_continuous), color=blue, ls='-')\n",
    "# Extrapolate linear model\n",
    "ax.plot(np.linspace(.9,1,50), fit_results.eval(x=np.linspace(.9,1,50)), color=blue, ls=':')\n",
    "\n",
    "# Annotate and plot solid target data \n",
    "ax.axvline(1, color='grey', ls='--')\n",
    "ax.text(.95, df_r['Sputter_Coeff_Mean'].iloc[0]/6, 'solid', rotation='vertical', color='gray')\n",
    "ax.scatter([1], flat, color=orange, zorder=5.)\n",
    "ax.scatter(df_r['Filling'], df_r['Sputter_Coeff_Mean'], color=red, zorder=6)\n",
    "ax.text(1.02, df_r['Sputter_Coeff_Mean'].iloc[0]*.98, '$Y_\\\\mathrm{r}$', color=red)\n",
    "ax.text(.94, flat*.95, '$Y_\\\\mathrm{f}$', color=orange)\n",
    "\n",
    "# Plot our model prediction\n",
    "ax.plot(rho, rho*exp_slope+exp_intercept, ls='--', color=pink)\n",
    "\n",
    "# Plot data from Lopez-Cazalilla et al.\n",
    "ax.errorbar([.58], lop_exp_example['yield'].iloc[0], yerr=lop_exp_example['error'].iloc[0], \n",
    "            fmt='o', color=green, capsize=2, markerfacecolor='w', lw=1, markeredgewidth=1)\n",
    "ax.scatter([1], lop_afm_example, color=green, zorder=10, facecolor='w')\n",
    "\n",
    "# Labels, Limits, Ticks\n",
    "ax_lab = 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.13)\n",
    "ax.set_ylim(0,None)\n",
    "ax.set_xticks([.2,.4,.6,.8,1])\n",
    "ax.set_yticks([0,.5,1,1.5])\n",
    "\n",
    "# Further annotations \n",
    "ax.text(.4,.81, 'Fuzz (this study)', color=blue, rotation=13)\n",
    "ax.text(.65,.3,'NCW exp.', color=green)\n",
    "ax.annotate('', (.58, .58),\n",
    "            xytext=(.65, .33), color=green,\n",
    "            arrowprops=dict(arrowstyle='-|>', color=green,\n",
    "                            connectionstyle='angle, angleA=-180, angleB=90'))\n",
    "ax.text(1.02, .7, 'NCW\\nAFM', color=green)\n",
    "\n",
    "# Set panel title\n",
    "ax.text(.05, .9, f'a) $1\\,$keV Ar $\\\\rightarrow$ W @ 45°', transform=ax.transAxes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0e89abf",
   "metadata": {},
   "source": [
    "## Figure 8b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ec01625",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create axis\n",
    "ax2 = subfig_b.add_subplot(gs[0,0])\n",
    "\n",
    "# Plot and annotate flat surface data \n",
    "ax2.plot(angles_sim, eckstein(np.deg2rad(angles_sim), *eckstein_params[1000]['Ar']), \n",
    "         color=orange)\n",
    "ax2.text(20,1.25, 'flat', color=orange, rotation=20)\n",
    "\n",
    "# Plot our model prediction\n",
    "ax2.plot(prediction_1keV['angle'], prediction_1keV['yield'], '--', \n",
    "         color=pink, lw=1.3, label='Prediction')\n",
    "\n",
    "# Plot Lopez-Cazalilla et al. data\n",
    "ax2.errorbar(lop_exp_1keV[lop_exp_1keV['angle'] <= 45]['angle'], \n",
    "             lop_exp_1keV[lop_exp_1keV['angle'] <= 45]['yield'], \n",
    "             yerr=lop_exp_1keV[lop_exp_1keV['angle'] <= 45]['error'],\n",
    "             fmt='o', color=green, capsize=2, markerfacecolor='None', \n",
    "             lw=1, markeredgewidth=1, elinewidth=0, label='NCW')\n",
    "\n",
    "# Labels, Limits\n",
    "ax2.set_ylabel('Sputter yield $Y$ (W/Ar)')\n",
    "ax2.set_xlim(None, 47)\n",
    "ax2.set_ylim(0)\n",
    "\n",
    "# Set panel title\n",
    "ax2.text(.05, .85, f'b) $1\\,$keV Ar', transform=ax2.transAxes)\n",
    "\n",
    "# Legend\n",
    "ax2.legend(frameon=False, ncols=2, fontsize='small', loc=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5175b5c",
   "metadata": {},
   "source": [
    "## Figure 8c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f03e582d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax3 = subfig_b.add_subplot(gs[0,1])\n",
    "\n",
    "# Plot flat surface data \n",
    "ax3.plot(angles_sim, eckstein(np.deg2rad(angles_sim), *eckstein_params[2000]['Ar']), \n",
    "         color=orange)\n",
    "\n",
    "# Plot model prediction\n",
    "ax3.plot(prediction_2keV['angle'], prediction_2keV['yield'], '--', color=pink, lw=1.3)\n",
    "\n",
    "# Plot Lopez-Cazalilla et al. data\n",
    "ax3.errorbar(lop_exp_2keV['angle'], lop_exp_2keV['yield'], yerr=lop_exp_2keV['error'],\n",
    "             fmt='o', color=green, capsize=2, markerfacecolor='None', \n",
    "             lw=1, markeredgewidth=1, elinewidth=0)\n",
    "\n",
    "# Limits\n",
    "ax3.set_xlim(None, 47)\n",
    "\n",
    "# Set panel title\n",
    "ax3.text(.05, .85, f'c) $2\\,$keV {projectile}', transform=ax3.transAxes)\n",
    "\n",
    "# set centered supylabel\n",
    "fig.canvas.draw_idle()\n",
    "ax2_left = ax2.get_position().x0\n",
    "ax3_right = ax3.get_position().x1\n",
    "label_pos = ax2_left + (ax3_right - ax2_left)/2\n",
    "\n",
    "subfig_b.supxlabel('Incidence angle $\\\\alpha$ (°)', x=label_pos, ha='center')\n",
    "\n",
    "# Save\n",
    "fig.savefig('Figures/comparison_lopez_cazalilla_et_al.pdf')\n",
    "fig.savefig('Figures/Figure_8.pdf')"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
