{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0fb7f66",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import matplotlib.patheffects as path_effects\n",
    "plt.style.use('../single_column.mplstyle')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f895c3c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = [1,2,3]\n",
    "pors = [1,2,3,4,5,6,7,8,9]\n",
    "rms_dict = {sample:{por/10:0 for por in pors} for sample in samples}\n",
    "siad_dict = {sample:{por/10:0 for por in pors} for sample in samples}\n",
    "deltam_dict = {sample:{por/10:0 for por in pors} for sample in samples}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5a61a92",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for ii, sample in enumerate(samples):\n",
    "    fig_surf = plt.figure(constrained_layout=True, figsize=(6.69, 13/2.54))\n",
    "    for jj, por in enumerate(pors):\n",
    "        df = pd.read_csv(f'Rough{sample}/0.{por}/surface_xyz.txt', \n",
    "                        names=['x', 'y', 'z'], sep='\\s+')\n",
    "        z_mean = df['z'].mean()\n",
    "        devs = [(z - z_mean)**2 for z in df['z']]\n",
    "        rms = np.sqrt(np.asarray(devs).mean())\n",
    "        rms_dict[sample][por/10] = rms\n",
    "        \n",
    "        siad = pd.read_csv(f'Rough{sample}/0.{por}/siad.txt',names=['ang', 'freq'],skiprows=3, sep='\\s+')\n",
    "        siad_dict[sample][por/10] = siad\n",
    "        delta_m = siad['ang']@siad['freq']/siad['freq'].sum()\n",
    "        deltam_dict[sample][por/10] = delta_m\n",
    "        \n",
    "        surface = df.pivot(index='y', columns='x', values='z')\n",
    "        x_coords = df['x'].unique()\n",
    "        y_coords = df['y'].unique()\n",
    "        X, Y = np.meshgrid(x_coords/10, y_coords/10)\n",
    "        \n",
    "        ax = fig_surf.add_subplot(3,3,jj+1)\n",
    "        #ax.ticklabel_format(axis='both', style='sci', scilimits=(0,0), useMathText=True)\n",
    "        pc = ax.pcolormesh(X, Y, surface/10,shading='gouraud')\n",
    "        ax.set_aspect('equal')\n",
    "        cb = fig_surf.colorbar(pc, ax=ax, shrink=1.0)#, label='$z$ (nm)')\n",
    "        if jj == 0: cb.ax.set_title('$z$ (nm)', fontsize=8)\n",
    "\n",
    "        txt = ax.text(0.05, .9, f'$\\\\rho_\\\\mathrm{{fil}}={por/10}$', transform=ax.transAxes)\n",
    "        txt.set_path_effects([\n",
    "            path_effects.Stroke(linewidth=4, foreground='white'),\n",
    "            path_effects.Normal()\n",
    "        ])\n",
    "        \n",
    "        txt = ax.text(.05, .07, f'$\\\\delta_\\\\mathrm{{m}}={deltam_dict[sample][por/10]:.0f}$°', transform=ax.transAxes)\n",
    "        txt.set_path_effects([\n",
    "            path_effects.Stroke(linewidth=3.5, foreground='white'),\n",
    "            path_effects.Normal()\n",
    "        ])\n",
    "        \n",
    "        txt = ax.text(.05, .16, f'RMS$={rms_dict[sample][por/10]/10:.0f}\\,$nm', transform=ax.transAxes)\n",
    "        txt.set_path_effects([\n",
    "            path_effects.Stroke(linewidth=3.5, foreground='white'),\n",
    "            path_effects.Normal()\n",
    "        ])\n",
    "        \n",
    "        if jj == 3 : ax.set_ylabel('$y$ (nm)')\n",
    "        if jj == 7 : ax.set_xlabel('$x$ (nm)')\n",
    "\n",
    "\n",
    "    fig_surf.savefig(f'Rough{sample}/surfaces{sample}')\n",
    "    if ii == 0: fig_surf.savefig(f'../Figures/Figure_2.pdf')\n",
    "    plt.close(fig_surf)\n",
    "    \n",
    "pors = np.asarray(pors)/10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ea083cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(constrained_layout=True, figsize=(6.69, 10/2.54))\n",
    "ax_list = []\n",
    "for ii, por in enumerate(pors):\n",
    "    ax = fig.add_subplot(3,3,ii+1)\n",
    "    \n",
    "    ax.plot(siad_dict[1][por]['ang'], siad_dict[1][por]['freq'], label='sample 1')\n",
    "    ax.plot(siad_dict[2][por]['ang'], siad_dict[2][por]['freq'], label='sample 2')\n",
    "    ax.plot(siad_dict[3][por]['ang'], siad_dict[3][por]['freq'], label='sample 3')\n",
    "    \n",
    "    ax.set_xlim(0,90)\n",
    "    ax.set_ylim(0, None)\n",
    "    ax.text(0.05, 0.85, f'$\\\\rho_\\\\mathrm{{fil}} = {por}$', transform=ax.transAxes)\n",
    "    if ii == 0: ax.legend()\n",
    "    if ii == 3: ax.set_ylabel('frequency')\n",
    "    if ii == 7: ax.set_xlabel('local inclination angle $\\\\theta$ (°)')\n",
    "    if ii not in [6,7,8]: ax.get_xaxis().set_ticks([20, 40, 60, 80], ['', '', '', ''])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c88eb968",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1,constrained_layout=True)\n",
    "cmap = mpl.colormaps['viridis']\n",
    "colors = cmap(np.linspace(0, 1, 10))\n",
    "\n",
    "for por, color in zip(pors, colors):\n",
    "    ax.plot(siad_dict[3][por]['ang'], siad_dict[3][por]['freq'], color=color, label=f'{por}')\n",
    "ax.set_xlim(1, 90)\n",
    "ax.set_ylim(0, .12)\n",
    "\n",
    "\n",
    "ins2 = ax.inset_axes([.25,0.35,.5,.6])\n",
    "ins2.set_xlim(37,53)\n",
    "ins2.set_ylim(None, .032)\n",
    "\n",
    "for por, color in zip(pors, colors):\n",
    "    ins2.plot(siad_dict[3][por]['ang'], siad_dict[3][por]['freq'], color=color, label=f'{por}')\n",
    "\n",
    "ax.indicate_inset_zoom(ins2, edgecolor=\"black\")\n",
    "ax.set_xlabel('local inclination angle $\\\\theta$ (°)')\n",
    "ax.set_ylabel('frequency')\n",
    "ax.legend( loc='upper left', ncols=1, title='$\\\\rho_\\\\mathrm{fil}$', \n",
    "          bbox_to_anchor = (1.0,1.027), frameon=False)\n",
    "fig.savefig('../Figures/sample3_siads')\n",
    "fig.savefig('../Figures/Figure_6.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
}