{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from parameters_U233 import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from jitcdde import t\n",
    "from msrDynamics import Node, System\n",
    "import pandas as pd\n",
    "import math\n",
    "import matplotlib.patches as mpatches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = 8\n",
    "\n",
    "# unpack ORNL data\n",
    "df_ORNL_mag = pd.read_csv(f\"./data/ORNL_msre_{int(P)}MW_U233_magnitude.csv\",names=['f','mag'])\n",
    "df_ORNL_mag = df_ORNL_mag.sort_values(df_ORNL_mag.columns[0])\n",
    "df_ORNL_phase = pd.read_csv(f\"./data/ORNL_msre_{int(P)}MW_U233_phase.csv\",names=['f','phase'])\n",
    "df_ORNL_phase = df_ORNL_phase.sort_values(df_ORNL_phase.columns[0])\n",
    "df_ORNL_mag_t = pd.read_csv(f\"./data/ORNL_msre_{int(P)}MW_U233_magnitude_theoretical.csv\",names=['f','mag'])\n",
    "df_ORNL_mag_t = df_ORNL_mag_t.sort_values(df_ORNL_mag_t.columns[0])\n",
    "df_ORNL_phase_t = pd.read_csv(f\"./data/ORNL_msre_{int(P)}MW_U233_phase_theoretical.csv\",names=['f','phase'])\n",
    "df_ORNL_phase_t = df_ORNL_phase_t.sort_values(df_ORNL_phase_t.columns[0])\n",
    "df_jit = pd.read_csv(f\"./data/frequency_response_results_{P}_MW.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 179304 characters]",
      "text/plain": [
       "<Figure size 600x1200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,axs = plt.subplots(2,1,figsize=(6,12))\n",
    "\n",
    "colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown']\n",
    "\n",
    "# function to update the style of each axis\n",
    "def update_axis_style(ax, title, x_label='', y_label='', x_ticks=True):\n",
    "    ax.set_title(title,fontsize=20)\n",
    "    ax.set_xlabel(x_label,fontsize=14)\n",
    "    ax.set_ylabel(y_label,fontsize=14)\n",
    "    ax.grid(True, which='both', linestyle='--', linewidth=0.5)\n",
    "    ax.tick_params(axis='x', which='both', bottom=x_ticks, top=False, labelbottom=x_ticks)\n",
    "    ax.tick_params(axis='y', which='both', left=True, right=False, labelleft=True)\n",
    "\n",
    "update_axis_style(axs[0],\"Magnitude\")\n",
    "axs[0].plot(df_jit['Frequency'],df_jit['Gain'],label=\"msrDynamics\")\n",
    "axs[0].plot(df_ORNL_mag_t['f'],df_ORNL_mag_t['mag'],label=\"ORNL: theoretical\",color=colors[1],linestyle=\"--\")\n",
    "axs[0].scatter(df_ORNL_mag['f'],df_ORNL_mag['mag'],marker=\"2\",color=colors[2],label=\"ORNL: measured\")\n",
    "axs[0].set_yscale(\"log\")\n",
    "axs[0].set_xscale(\"log\")\n",
    "axs[0].set_xlim([9e-3,1e0])\n",
    "axs[0].set_ylim([1e2,3.5e3])\n",
    "axs[0].set_ylabel(r\"$\\frac{\\delta n}{\\delta \\rho n_0}$\")\n",
    "axs[0].legend()\n",
    "\n",
    "update_axis_style(axs[1],\"Phase\",x_label='Frequency (rad/s)',y_label=\"Phase(deg)\")\n",
    "axs[1].plot(df_jit['Frequency'],df_jit['Phase Shift'],label=\"msrDynamics\")\n",
    "axs[1].plot(df_ORNL_phase_t['f'],df_ORNL_phase_t['phase'],label=\"ORNL: theoretical\",color=colors[1],linestyle=\"--\")\n",
    "axs[1].scatter(df_ORNL_phase['f'],df_ORNL_phase['phase'],marker=\"2\",color=colors[2],label=\"ORNL-TM-2997\")\n",
    "axs[1].set_xscale(\"log\")\n",
    "axs[1].set_xlim([9e-3,1e0])\n",
    "axs[1].legend()\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_adj = pd.read_csv('frequency_response_results_8_MW_double_tau.csv')\n",
    "\n",
    "\n",
    "fig,axs = plt.subplots(2,1,figsize=(10,10))\n",
    "\n",
    "colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown']\n",
    "\n",
    "# function to update the style of each axis\n",
    "def update_axis_style(ax, title, x_label='', y_label='', x_ticks=True):\n",
    "    ax.set_title(title)\n",
    "    ax.set_xlabel(x_label)\n",
    "    ax.set_ylabel(y_label)\n",
    "    ax.grid(True, which='both', linestyle='--', linewidth=0.5)\n",
    "    ax.tick_params(axis='x', which='both', bottom=x_ticks, top=False, labelbottom=x_ticks)\n",
    "    ax.tick_params(axis='y', which='both', left=True, right=False, labelleft=True)\n",
    "\n",
    "update_axis_style(axs[0],r\"Magnitude, $\\tau = \\tau_l + \\tau_c$\")\n",
    "axs[0].plot(df_jit['Frequency'],df_jit['Gain'],label=r\"$\\tau = \\tau_l + \\tau_c$\") \n",
    "axs[0].set_yscale(\"log\")\n",
    "axs[0].set_xscale(\"log\")\n",
    "axs[0].set_ylabel(r\"$\\frac{\\delta n}{\\delta \\rho n_0}$\")\n",
    "axs[0].set_xlim([1e-2,2e0])\n",
    "p = tau_c+tau_l\n",
    "f_dip = (2*math.pi)/(p)\n",
    "f_dip2 = (2*math.pi)/((p)/2)\n",
    "f_peak = (2*math.pi)/(2*(p))\n",
    "f_peak2 = (2*math.pi)/((2/3)*(p))\n",
    "axs[0].axvline(x=f_peak,linestyle=\"--\",label=r\"$f = 2\\pi / (2\\tau)$\",color=colors[1])\n",
    "axs[0].axvline(x=f_dip,linestyle=\"--\",label=r\"$f = 2\\pi / \\tau$\",color=colors[2])\n",
    "axs[0].axvline(x=f_peak2,linestyle=\"--\",label=r\"$f = 2\\pi / (2\\tau/3)$\",color=colors[4])\n",
    "axs[0].axvline(x=f_dip2,linestyle=\"--\",label=r\"$f = 2\\pi / (\\tau/2)$\",color=colors[3])\n",
    "axs[0].legend()\n",
    "\n",
    "update_axis_style(axs[1],r\"Magnitude, $\\tau = 2\\tau_l + 2\\tau_c$\")\n",
    "axs[1].plot(df_adj['Frequency'],df_adj['Gain'],label=r\"$\\tau = 2(\\tau_l + \\tau_c)$\") \n",
    "axs[1].set_yscale(\"log\")\n",
    "axs[1].set_xscale(\"log\")\n",
    "axs[1].set_ylabel(r\"$\\frac{\\delta n}{\\delta \\rho n_0}$\")\n",
    "axs[1].set_xlim([1e-2,2e0])\n",
    "p2 = 2*(tau_c+tau_l)\n",
    "f_dip = (2*math.pi)/(p2)\n",
    "f_dip2 = (2*math.pi)/((p2)/2)\n",
    "f_peak = (2*math.pi)/(2*(p2))\n",
    "f_peak2 = (2*math.pi)/((2/3)*(p2))\n",
    "axs[1].axvline(x=f_peak,linestyle=\"--\",label=r\"$f = 2\\pi / (2\\tau)$\",color=colors[1])\n",
    "axs[1].axvline(x=f_dip,linestyle=\"--\",label=r\"$f = 2\\pi / \\tau$\",color=colors[2])\n",
    "axs[1].axvline(x=f_peak2,linestyle=\"--\",label=r\"$f = 2\\pi / (2\\tau/3)$\",color=colors[4])\n",
    "axs[1].axvline(x=f_dip2,linestyle=\"--\",label=r\"$f = 2\\pi / (\\tau/2)$\",color=colors[3])\n",
    "axs[1].legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_adj = pd.read_csv('frequency_response_results_8_MW_double_tau.csv')\n",
    "\n",
    "\n",
    "fig,axs = plt.subplots(2,1,figsize=(10,10))\n",
    "\n",
    "colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown']\n",
    "\n",
    "# function to update the style of each axis\n",
    "def update_axis_style(ax, title, x_label='', y_label='', x_ticks=True):\n",
    "    ax.set_title(title)\n",
    "    ax.set_xlabel(x_label)\n",
    "    ax.set_ylabel(y_label)\n",
    "    ax.grid(True, which='both', linestyle='--', linewidth=0.5)\n",
    "    ax.tick_params(axis='x', which='both', bottom=x_ticks, top=False, labelbottom=x_ticks)\n",
    "    ax.tick_params(axis='y', which='both', left=True, right=False, labelleft=True)\n",
    "\n",
    "# Filling between ranges where inc_dec is True\n",
    "inc_dec = np.ones(len(df_jit))\n",
    "for f in enumerate(df_jit['Frequency']):\n",
    "    delay_input = np.cos(f[1]*((-p))) # t = 2pi-tau\n",
    "    if (delay_input > 0):\n",
    "        inc_dec[f[0]] = 0\n",
    "\n",
    "# Find start and end points of True and False segments in inc_dec and fill between them\n",
    "start_true = None\n",
    "start_false = None\n",
    "for i in range(len(df_jit['Frequency'])):\n",
    "    if inc_dec[i]:\n",
    "        if start_false is not None:\n",
    "            # Fill the False segment before starting the True segment\n",
    "            axs[0].fill_between(df_jit['Frequency'][start_false:i], 0, df_jit['Gain'][start_false:i], alpha=0.23, color='red',interpolate=True)\n",
    "            start_false = None\n",
    "        if start_true is None:\n",
    "            start_true = i  # Mark the start of a True segment\n",
    "    else:\n",
    "        if start_true is not None:\n",
    "            # Fill the True segment before starting the False segment\n",
    "            axs[0].fill_between(df_jit['Frequency'][start_true:i], 0, df_jit['Gain'][start_true:i], alpha=0.23, color='green',interpolate=True)\n",
    "            start_true = None\n",
    "        if start_false is None:\n",
    "            start_false = i  # Mark the start of a False segment\n",
    "\n",
    "# Check if there is a segment that goes until the end\n",
    "if start_true is not None:\n",
    "    axs[0].fill_between(df_jit['Frequency'][start_true:], 0, df_jit['Gain'][start_true:], alpha=0.3, color='green',interpolate=True)\n",
    "if start_false is not None:\n",
    "    axs[0].fill_between(df_jit['Frequency'][start_false:], 0, df_jit['Gain'][start_false:], alpha=0.3, color='red',interpolate=True)\n",
    "\n",
    "update_axis_style(axs[0],r\"Magnitude, $\\tau = \\tau_l + \\tau_c$\")\n",
    "axs[0].plot(df_jit['Frequency'],df_jit['Gain'],label=r\"$\\tau = \\tau_l + \\tau_c$\") \n",
    "axs[0].set_yscale(\"log\")\n",
    "axs[0].set_xscale(\"log\")\n",
    "axs[0].set_ylabel(r\"$\\frac{\\delta n}{\\delta \\rho n_0}$\")\n",
    "axs[0].set_xlim([1e-2,2e0])\n",
    "\n",
    "red_patch = mpatches.Patch(color='red', alpha=0.3, label=r\"$\\frac{d\\rho_{in}}{dt}(t-2 \\pi /f) \\cdot \\frac{d\\rho_{in}}{dt}(t) > 0$\")\n",
    "green_patch = mpatches.Patch(color='green', alpha=0.3, label=r\"$\\frac{d\\rho_{in}}{dt}(t-2 \\pi /f) \\cdot \\frac{d\\rho_{in}}{dt}(t) < 0$\")\n",
    "\n",
    "axs[0].legend(handles=[axs[0].get_lines()[0], red_patch, green_patch])\n",
    "\n",
    "######################################################################################\n",
    "# Double tau\n",
    "######################################################################################\n",
    "\n",
    "update_axis_style(axs[1],r\"Magnitude, $\\tau = 2\\tau_l + 2\\tau_c$\")\n",
    "axs[1].plot(df_adj['Frequency'],df_adj['Gain'],label=r\"$\\tau = 2(\\tau_l + \\tau_c)$\") \n",
    "axs[1].set_yscale(\"log\")\n",
    "axs[1].set_xscale(\"log\")\n",
    "axs[1].set_ylabel(r\"$\\frac{\\delta n}{\\delta \\rho n_0}$\")\n",
    "axs[1].set_xlim([1e-2,2e0])\n",
    "axs[1].legend(handles=[axs[1].get_lines()[0], red_patch, green_patch])\n",
    "\n",
    "# Filling between ranges where inc_dec is True\n",
    "inc_dec = np.ones(len(df_adj))\n",
    "for f in enumerate(df_adj['Frequency']):\n",
    "    delay_input = np.cos(f[1]*(-p2)) # t = 2pi-tau\n",
    "    if (delay_input > 0):\n",
    "        inc_dec[f[0]] = 0\n",
    "\n",
    "# Find start and end points of True and False segments in inc_dec and fill between them\n",
    "start_true = None\n",
    "start_false = None\n",
    "for i in range(len(df_adj['Frequency'])):\n",
    "    if inc_dec[i]:\n",
    "        if start_false is not None:\n",
    "            # Fill the False segment before starting the True segment\n",
    "            axs[1].fill_between(df_adj['Frequency'][start_false:i], 0, df_adj['Gain'][start_false:i], alpha=0.23, color='red',interpolate=True)\n",
    "            start_false = None\n",
    "        if start_true is None:\n",
    "            start_true = i  # Mark the start of a True segment\n",
    "    else:\n",
    "        if start_true is not None:\n",
    "            # Fill the True segment before starting the False segment\n",
    "            axs[1].fill_between(df_adj['Frequency'][start_true:i], 0, df_adj['Gain'][start_true:i], alpha=0.23, color='green',interpolate=True)\n",
    "            start_true = None\n",
    "        if start_false is None:\n",
    "            start_false = i  # Mark the start of a False segment\n",
    "\n",
    "# Check if there is a segment that goes until the end\n",
    "if start_true is not None:\n",
    "    axs[1].fill_between(df_adj['Frequency'][start_true:], 0, df_adj['Gain'][start_true:], alpha=0.3, color='green',interpolate=True)\n",
    "if start_false is not None:\n",
    "    axs[1].fill_between(df_adj['Frequency'][start_false:], 0, df_adj['Gain'][start_false:], alpha=0.3, color='red',interpolate=True)\n",
    "\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "thesis_env",
   "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.11.6"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
