Figure_5_Figure_6_post-processing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "This Source Code Form is subject to the terms of the Mozilla Public\n",
    "License, v. 2.0. If a copy of the MPL was not distributed with this\n",
    "file, You can obtain one at http://mozilla.org/MPL/2.0/.\n",
    "\n",
    "\n",
    "Authors:\n",
    "Roland Baatz roland.baatz @ zalf.de\n",
    "\n",
    "Maintainers and contact:\n",
    "Currently maintained by the authors.\n",
    "\n",
    "This file requires output of the CRNS_*.ipynb\n",
    "\n",
    "Copyright (C) Leibniz Centre for Agricultural Landscape Research (ZALF)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import time\n",
    "import scipy.io as sio\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#dfi_moderated = pd.read_csv(csv_filename, sep=';', index_col=0)  # Set index_col to 0 for row names\n",
    "### Iterate over EU Study sites:\n",
    "plog_flag_id=0\n",
    "src_path = './data'\n",
    "suffix_csv = \"rmse_u\"\n",
    "csv_filename=src_path+\"my_dataframe_moderated_\"+suffix_csv+'.csv'\n",
    "dfi_moderated = pd.read_csv(csv_filename, sep=';', index_col=0)  # Set index_col to 0 for row names\n",
    "\n",
    "#bare counter result file\n",
    "src_path = './data'\n",
    "suffix_csv = \"rmse_u\"\n",
    "csv_filename=src_path+\"my_dataframe_slow_\"+suffix_csv+'-.csv'\n",
    "dfi_slow = pd.read_csv(csv_filename, sep=';', index_col=0)  # Set index_col to 0 for row names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_rmse(ts1,ts2):\n",
    "    return np.sqrt(np.sum((ts1 - ts2)**2))\n",
    "\n",
    "def calc_mae(ts1,ts2):\n",
    "    return np.mean(abs(ts1 - ts2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Assuming dfi_moderated and dfi_slow are already defined\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 4))\n",
    "\n",
    "# First subplot\n",
    "ax1 = plt.subplot(1, 3, 1)\n",
    "ax1.scatter(dfi_moderated['beta'], dfi_slow['beta'], color='k', marker='*', alpha=0.7)\n",
    "print(\"MAE beta=\",calc_mae(dfi_moderated['beta'], dfi_slow['beta']))\n",
    "# Get the limits of the x and y axes\n",
    "x_limits = ax1.get_xlim()\n",
    "y_limits = ax1.get_ylim()\n",
    "\n",
    "# Determine the range for the 1:1 line\n",
    "min_limit = min(x_limits[0], y_limits[0])\n",
    "max_limit = max(x_limits[1], y_limits[1])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax1.plot([min_limit, max_limit], [min_limit, max_limit], 'r--', label='1-to-1 line')\n",
    "\n",
    "# Set the limits again to ensure the plot is correct\n",
    "#ax1.set_xlim(x_limits)\n",
    "#ax1.set_ylim(y_limits)\n",
    "\n",
    "# Rotate the x-axis labels\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "\n",
    "ax1.set_xlabel('beta (epithermal neutrons)')\n",
    "ax1.set_ylabel('beta (thermal neutrons)')\n",
    "\n",
    "ax1.legend()\n",
    "ax1.set_title('a) beta - atmosph. press. scaling', loc='left')\n",
    "\n",
    "# Second subplot\n",
    "ax2 = plt.subplot(1, 3, 2)\n",
    "ax2.scatter(dfi_moderated['omega'], dfi_slow['omega'], color='k', marker='*', alpha=0.7)\n",
    "print(\"MAE omega=\",calc_mae(dfi_moderated['omega'], dfi_slow['omega']))\n",
    "\n",
    "# Get the limits of the x and y axes\n",
    "x_limits = ax2.get_xlim()\n",
    "y_limits = ax2.get_ylim()\n",
    "\n",
    "# Determine the range for the 1:1 line\n",
    "min_limit = min(x_limits[0], y_limits[0])\n",
    "max_limit = max(x_limits[1], y_limits[1])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax2.plot([min_limit, max_limit], [min_limit, max_limit], 'r--')\n",
    "\n",
    "# Set the limits again to ensure the plot is correct\n",
    "#ax2.set_xlim(x_limits)\n",
    "#ax2.set_ylim(y_limits)\n",
    "\n",
    "# Rotate the x-axis labels\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "\n",
    "ax2.set_xlabel('omega (epithermal neutrons)')\n",
    "ax2.set_ylabel('omega (thermal neutrons)')\n",
    "ax2.set_title('b) omega - air humidity scaling', loc='left')\n",
    "\n",
    "# Third subplot\n",
    "ax3 = plt.subplot(1, 3, 3)\n",
    "ax3.scatter(dfi_moderated['psi'], dfi_slow['psi'], color='k', marker='*', alpha=0.7)\n",
    "print(\"MAE psi=\",calc_mae(dfi_moderated['psi'], dfi_slow['psi']))\n",
    "\n",
    "# Get the limits of the x and y axes\n",
    "x_limits = ax3.get_xlim()\n",
    "y_limits = ax3.get_ylim()\n",
    "\n",
    "# Determine the range for the 1:1 line\n",
    "min_limit = min(x_limits[0], y_limits[0])\n",
    "max_limit = max(x_limits[1], y_limits[1])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax3.plot([min_limit, max_limit], [min_limit, max_limit], 'r--')\n",
    "\n",
    "# Set the limits again to ensure the plot is correct\n",
    "#ax3.set_xlim(x_limits)\n",
    "#ax3.set_ylim(y_limits)\n",
    "\n",
    "# Rotate the x-axis labels\n",
    "ax3.tick_params(axis='x', rotation=45)\n",
    "\n",
    "ax3.set_xlabel('psi (epithermal neutrons)')\n",
    "ax3.set_ylabel('psi (thermal neutrons)')\n",
    "ax3.set_title('c) psi - incoming neutron scaling', loc='left')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "pngname=\"./plots/Figure_6_Scatter_Moderated_vs_Slow.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Assuming dfi_moderated and dfi_slow are already defined\n",
    "fig, axes = plt.subplots(1, 2, figsize=(8, 4))\n",
    "\n",
    "# First subplot\n",
    "ax1 = plt.subplot(1, 2, 1)\n",
    "xdata = dfi_moderated['number_of_days_for_regression']\n",
    "ydata = dfi_moderated['beta_sd']\n",
    "ax1.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "ax1.plot(xaxis, curve_y, color='grey', label='Fitted Line')\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "\n",
    "\n",
    "ax1.set_xlabel('Number of days', loc='left')\n",
    "ax1.set_ylabel('Uncertainty')\n",
    "\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "# First subplot\n",
    "ax2 = plt.subplot(1, 2, 2)\n",
    "xdata = dfi_moderated['number_of_days_for_regression']*dfi_moderated['mean_N']*24\n",
    "ydata = dfi_moderated['beta_sd']\n",
    "ax2.scatter(xdata, ydata, color='k', marker='*', alpha=0.7)\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "function_text = f'$f(x) = {coefficients[0]:.8f}x^3 + {coefficients[1]:.8f}x^2 + {coefficients[2]:.8f}x + {coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^{{3}}+{coefficients[1]:.8f}x^{{2}}+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "#ax2.plot(xaxis, curve_y, color='grey', label=function_text)\n",
    "ax2.plot(xaxis, curve_y, color='grey')\n",
    "\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "\n",
    "\n",
    "ax2.set_xlabel('Total neutron counts', loc='left')\n",
    "ax2.set_ylabel('Uncertainty')\n",
    "#ax2.legend()\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "pngname=\"./plots/Moderated_uncertainty_beta.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Assuming dfi_moderated and dfi_slow are already defined\n",
    "fig, axes = plt.subplots(1, 2, figsize=(8, 4))\n",
    "\n",
    "# First subplot\n",
    "ax1 = plt.subplot(1, 2, 1)\n",
    "xdata = dfi_moderated['number_of_days_for_regression']\n",
    "ydata = dfi_moderated['omega_sd']\n",
    "ax1.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "ax1.plot(xaxis, curve_y, color='grey', label='Fitted Line')\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "\n",
    "\n",
    "ax1.set_xlabel('Number of days', loc='left')\n",
    "ax1.set_ylabel('Uncertainty')\n",
    "\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "# First subplot\n",
    "ax2 = plt.subplot(1, 2, 2)\n",
    "xdata = dfi_moderated['number_of_days_for_regression']*dfi_moderated['mean_N']*24\n",
    "ydata = dfi_moderated['omega_sd']\n",
    "ax2.scatter(xdata, ydata, color='k', marker='*', alpha=0.7)\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "function_text = f'$f(x) = {coefficients[0]:.8f}x^3 + {coefficients[1]:.8f}x^2 + {coefficients[2]:.8f}x + {coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^{{3}}+{coefficients[1]:.8f}x^{{2}}+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "#ax2.plot(xaxis, curve_y, color='grey', label=function_text)\n",
    "ax2.plot(xaxis, curve_y, color='grey')\n",
    "\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "\n",
    "\n",
    "ax2.set_xlabel('Total neutron counts', loc='left')\n",
    "ax2.set_ylabel('Uncertainty')\n",
    "#ax2.legend()\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "pngname=\"./plots/Moderated_uncertainty_omega.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Assuming dfi_moderated and dfi_slow are already defined\n",
    "fig, axes = plt.subplots(1, 2, figsize=(8, 4))\n",
    "\n",
    "# First subplot\n",
    "ax1 = plt.subplot(1, 2, 1)\n",
    "xdata = dfi_moderated['number_of_days_for_regression']\n",
    "ydata = dfi_moderated['psi_sd']\n",
    "ax1.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "ax1.plot(xaxis, curve_y, color='grey', label='Fitted Line')\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "\n",
    "\n",
    "ax1.set_xlabel('Number of days', loc='left')\n",
    "ax1.set_ylabel('Uncertainty')\n",
    "\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "# First subplot\n",
    "ax2 = plt.subplot(1, 2, 2)\n",
    "xdata = dfi_moderated['number_of_days_for_regression']*dfi_moderated['mean_N']*24\n",
    "ydata = dfi_moderated['psi_sd']\n",
    "ax2.scatter(xdata, ydata, color='k', marker='*', alpha=0.7)\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "function_text = f'$f(x) = {coefficients[0]:.8f}x^3 + {coefficients[1]:.8f}x^2 + {coefficients[2]:.8f}x + {coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^{{3}}+{coefficients[1]:.8f}x^{{2}}+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "#ax2.plot(xaxis, curve_y, color='grey', label=function_text)\n",
    "ax2.plot(xaxis, curve_y, color='grey')\n",
    "\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "\n",
    "\n",
    "ax2.set_xlabel('Total neutron counts', loc='left')\n",
    "ax2.set_ylabel('Uncertainty')\n",
    "#ax2.legend()\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "pngname=\"./plots/Moderated_uncertainty_psi.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Assuming dfi_moderated and dfi_slow are already defined\n",
    "fig, axes = plt.subplots(3, 2, figsize=(8, 4*3))\n",
    "\n",
    "# First subplot\n",
    "ax1 = axes[0, 0]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']\n",
    "ydata = dfi_moderated['beta_sd']\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "ax1.plot(xaxis, curve_y, color='grey', label='Fitted Line')\n",
    "ax1.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "ax1.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1)\n",
    "ax1.scatter(1315,  0.00033, color='b', marker='+',s=93, alpha=1) #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "ax1.set_ylabel('Uncertainty')\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "\n",
    "\n",
    "# First subplot\n",
    "ax2 = axes[0, 1]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']*dfi_moderated['mean_N']*24\n",
    "ydata = dfi_moderated['beta_sd']\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "function_text = f'$f(x) = {coefficients[0]:.8f}x^3 + {coefficients[1]:.8f}x^2 + {coefficients[2]:.8f}x + {coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^{{3}}+{coefficients[1]:.8f}x^{{2}}+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "#ax2.plot(xaxis, curve_y, color='grey', label=function_text)\n",
    "\n",
    "ax2.plot(xaxis, curve_y, color='grey')\n",
    "ax2.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "ax2.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1)\n",
    "ax2.scatter(1315*24*445.903,0.00033, color='b', marker='+', s=93, alpha=1) #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "ax1.set_title('a) Beta uncertainty (days obs.)', loc='left')\n",
    "ax2.set_title('b) Beta uncertainty (neutrons obs.)', loc='left')\n",
    "\n",
    "#ax2.legend()\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "ax1 = axes[1, 0]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']\n",
    "ydata = dfi_moderated['omega_sd']\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "\n",
    "ax1.plot(xaxis, curve_y, color='grey', label='Fitted line')\n",
    "ax1.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7, label='COSMOS EU')\n",
    "ax1.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1, label='Alento(2)')\n",
    "ax1.scatter(1315,0.0013, color='b', marker='+', alpha=1, s=93, label='Synthetic case') #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "ax1.legend()\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "ax1.set_ylabel('Uncertainty')\n",
    "\n",
    "# First subplot\n",
    "ax2 = axes[1, 1]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']*dfi_moderated['mean_N']*24\n",
    "ydata = dfi_moderated['omega_sd']\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "function_text = f'$f(x) = {coefficients[0]:.8f}x^3 + {coefficients[1]:.8f}x^2 + {coefficients[2]:.8f}x + {coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^{{3}}+{coefficients[1]:.8f}x^{{2}}+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "#ax2.plot(xaxis, curve_y, color='grey', label=function_text)\n",
    "ax2.plot(xaxis, curve_y, color='grey')\n",
    "ax2.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "ax2.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1)\n",
    "ax2.scatter(1315*445.903*24,0.0013, color='b', marker='+',s=93, alpha=1) #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "\n",
    "ax1.set_title('c) Omega uncertainty (days obs.)', loc='left')\n",
    "ax2.set_title('d) Omega uncertainty (neutrons obs.)', loc='left')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# First subplot\n",
    "ax1 = axes[2, 0]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']\n",
    "ydata = dfi_moderated['psi_sd']\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "ax1.plot(xaxis, curve_y, color='grey', label='Fitted Line')\n",
    "ax1.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "ax1.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1)\n",
    "ax1.scatter(1315,0.18, color='b', marker='+',s=93, alpha=1) #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "ax1.set_ylabel('Uncertainty')\n",
    "ax1.set_xlabel('Number of days', loc='left')\n",
    "\n",
    "# First subplot\n",
    "ax2 = axes[2, 1]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']*dfi_moderated['mean_N']*24\n",
    "ydata = dfi_moderated['psi_sd']\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "function_text = f'$f(x) = {coefficients[0]:.8f}x^3 + {coefficients[1]:.8f}x^2 + {coefficients[2]:.8f}x + {coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^{{3}}+{coefficients[1]:.8f}x^{{2}}+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "#ax2.plot(xaxis, curve_y, color='grey', label=function_text)\n",
    "\n",
    "ax2.plot(xaxis, curve_y, color='grey')\n",
    "ax2.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "ax2.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1)\n",
    "ax2.scatter(1315*445.903*24,0.18, color='b', marker='+',s=93, alpha=1) #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "\n",
    "ax1.set_title('e) Psi uncertainty (days obs.)', loc='left')\n",
    "ax2.set_title('f) Psi uncertainty (neutrons obs.)', loc='left')\n",
    "\n",
    "ax2.set_xlabel('Total neutron counts', loc='left')\n",
    "#ax2.legend()\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "pngname=\"./plots/Moderated_uncertainty_all.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "from scipy.optimize import least_squares\n",
    "# Define the objective function to minimize (sum of absolute errors)\n",
    "def objective_function(params, x, y):\n",
    "    a, b = params\n",
    "    return np.sum(np.abs(y - logarithmic_function(x, a, b)))\n",
    "\n",
    "# Define a logarithmic function for fitting\n",
    "def logarithmic_function(x, a, b):\n",
    "    return a * np.log(x) + b\n",
    "\n",
    "# Define a residual function for least squares fitting\n",
    "def residuals(params, x, y):\n",
    "    a, b = params\n",
    "    return y - logarithmic_function(x, a, b)\n",
    "\n",
    "\n",
    "# Define an exponential function for fitting\n",
    "def exponential_function(x, a, b):\n",
    "    return a * np.exp(b * x)\n",
    "    \n",
    "# Assuming dfi_moderated and dfi_slow are already defined\n",
    "fig, axes = plt.subplots(3, 2, figsize=(8, 4*3))\n",
    "\n",
    "# First subplot\n",
    "ax1 = axes[0, 0]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']\n",
    "ydata = dfi_moderated['beta_sd']\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "# Initial guess for parameters\n",
    "initial_guess = [1, 1]\n",
    "\n",
    "# Perform robust fitting using least squares with Huber loss\n",
    "result = least_squares(residuals, initial_guess, args=(xdata, ydata), loss='huber')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "\n",
    "#ax1.plot(xaxis, curve_y, color='grey', label='Fitted Line')\n",
    "\n",
    "# Perform robust fitting using least squares with Soft L1 loss (LAD approximation)\n",
    "result = least_squares(residuals, initial_guess, args=(xdata, ydata), loss='soft_l1')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "\n",
    "# Perform the minimization using absolute error\n",
    "result = minimize(objective_function, initial_guess, args=(xdata, ydata), method='Nelder-Mead')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "\n",
    "ax1.plot(xaxis, curve_y, color='grey', label='Fitted Line')\n",
    "\n",
    "ax1.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7, label='COSMOS EU')\n",
    "ax1.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1, label='Alento(2)')\n",
    "ax1.scatter(1315,  0.00023, color='b', marker='+',s=93, alpha=1, label='Synthetic case') #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "ax1.legend()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ax1.set_ylabel('Uncertainty')\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "\n",
    "\n",
    "# First subplot\n",
    "ax2 = axes[0, 1]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']*dfi_moderated['mean_N']*24\n",
    "ydata = dfi_moderated['beta_sd']\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "function_text = f'$f(x) = {coefficients[0]:.8f}x^3 + {coefficients[1]:.8f}x^2 + {coefficients[2]:.8f}x + {coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^{{3}}+{coefficients[1]:.8f}x^{{2}}+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "#ax2.plot(xaxis, curve_y, color='grey', label=function_text)\n",
    "\n",
    "# Perform the minimization using absolute error\n",
    "result = minimize(objective_function, initial_guess, args=(xdata, ydata), method='Nelder-Mead')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "\n",
    "ax2.plot(xaxis, curve_y, color='grey')\n",
    "ax2.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "ax2.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1)\n",
    "ax2.scatter(1315*24*445.903,0.00023, color='b', marker='+', s=93, alpha=1) #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "ax1.set_title('a) Beta uncertainty (days obs.)', loc='left')\n",
    "ax2.set_title('b) Beta uncertainty (neutrons obs.)', loc='left')\n",
    "\n",
    "#ax2.legend()\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "ax1 = axes[1, 0]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']\n",
    "ydata = dfi_moderated['omega_sd']\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Finding the index of the dot with the highest uncertainty (highest omega_sd)\n",
    "max_uncertainty_index = dfi_moderated['omega_sd'].idxmax()\n",
    "# Removing the row with the highest uncertainty\n",
    "ydata = dfi_moderated.drop(max_uncertainty_index)\n",
    "#xdata = xdata.drop(max_uncertainty_index)\n",
    "\n",
    "\n",
    "# Ensure xdata and ydata are numeric\n",
    "xdata = dfi_moderated['number_of_days_for_regression'].astype(float)\n",
    "ydata = dfi_moderated['omega_sd'].astype(float)\n",
    "\n",
    "# Drop the row with the highest uncertainty for both xdata and ydata\n",
    "xdata = xdata.drop(max_uncertainty_index)\n",
    "ydata = ydata.drop(max_uncertainty_index)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "\n",
    "# Initial guess for parameters\n",
    "initial_guess = [1, 1]\n",
    "\n",
    "# Perform robust fitting using least squares with Huber loss\n",
    "result = least_squares(residuals, initial_guess, args=(xdata, ydata), loss='huber')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "\n",
    "#ax1.plot(xaxis, curve_y, color='grey', label='Fitted line')\n",
    "\n",
    "\n",
    "# Perform the minimization using absolute error\n",
    "result = minimize(objective_function, initial_guess, args=(xdata, ydata), method='Nelder-Mead')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "\n",
    "ax1.plot(xaxis, curve_y, color='grey', label='Fitted Line')\n",
    "\n",
    "ax1.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7, label='COSMOS EU')\n",
    "ax1.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1, label='Alento(2)')\n",
    "ax1.scatter(1315,0.0013, color='b', marker='+', alpha=1, s=93, label='Synthetic case') #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "#ax1.legend()\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "ax1.set_ylabel('Uncertainty')\n",
    "\n",
    "# First subplot\n",
    "ax2 = axes[1, 1]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']*dfi_moderated['mean_N']*24\n",
    "ydata = dfi_moderated['omega_sd']\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Finding the index of the dot with the highest uncertainty (highest omega_sd)\n",
    "max_uncertainty_index = dfi_moderated['omega_sd'].idxmax()\n",
    "# Removing the row with the highest uncertainty\n",
    "ydata = dfi_moderated.drop(max_uncertainty_index)\n",
    "# Calculate xdata as a pandas Series and then drop the same index\n",
    "xdata = (dfi_moderated['number_of_days_for_regression'] * dfi_moderated['mean_N'] * 24)\n",
    "xdata = xdata.drop(max_uncertainty_index)\n",
    "\n",
    "\n",
    "\n",
    "# Ensure xdata and ydata are numeric\n",
    "xdata = (dfi_moderated['number_of_days_for_regression'] * dfi_moderated['mean_N'] * 24).astype(float)\n",
    "ydata = dfi_moderated['omega_sd'].astype(float)\n",
    "\n",
    "# Drop the row with the highest uncertainty for both xdata and ydata\n",
    "xdata = xdata.drop(max_uncertainty_index)\n",
    "ydata = ydata.drop(max_uncertainty_index)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "function_text = f'$f(x) = {coefficients[0]:.8f}x^3 + {coefficients[1]:.8f}x^2 + {coefficients[2]:.8f}x + {coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^{{3}}+{coefficients[1]:.8f}x^{{2}}+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "\n",
    "\n",
    "\n",
    "# Perform the minimization using absolute error\n",
    "result = minimize(objective_function, initial_guess, args=(xdata, ydata), method='Nelder-Mead')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "\n",
    "\n",
    "#ax2.plot(xaxis, curve_y, color='grey', label=function_text)\n",
    "ax2.plot(xaxis, curve_y, color='grey')\n",
    "ax2.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "ax2.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1)\n",
    "ax2.scatter(1315*445.903*24,0.0013, color='b', marker='+',s=93, alpha=1) #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "\n",
    "ax1.set_title('c) Omega uncertainty (days obs.)', loc='left')\n",
    "ax2.set_title('d) Omega uncertainty (neutrons obs.)', loc='left')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# First subplot\n",
    "ax1 = axes[2, 0]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']\n",
    "ydata = dfi_moderated['psi_sd']\n",
    "\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "# Initial guess for parameters\n",
    "initial_guess = [1, 1]\n",
    "\n",
    "# Perform robust fitting using least squares with Huber loss\n",
    "result = least_squares(residuals, initial_guess, args=(xdata, ydata), loss='huber')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "\n",
    "\n",
    "\n",
    "ax1.plot(xaxis, curve_y, color='grey', label='Fitted Line')\n",
    "ax1.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "ax1.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1)\n",
    "ax1.scatter(1315,0.15, color='b', marker='+',s=93, alpha=1) #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "ax1.set_ylabel('Uncertainty')\n",
    "ax1.set_xlabel('Number of days', loc='left')\n",
    "\n",
    "# First subplot\n",
    "ax2 = axes[2, 1]\n",
    "xdata = dfi_moderated['number_of_days_for_regression']*dfi_moderated['mean_N']*24\n",
    "ydata = dfi_moderated['psi_sd']\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = p(xaxis)\n",
    "\n",
    "function_text = f'$y={coefficients[0]:.8f}x^3+{coefficients[1]:.8f}x^2+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "function_text = f'$f(x) = {coefficients[0]:.8f}x^3 + {coefficients[1]:.8f}x^2 + {coefficients[2]:.8f}x + {coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "function_text = f'$y={coefficients[0]:.8f}x^{{3}}+{coefficients[1]:.8f}x^{{2}}+{coefficients[2]:.8f}x+{coefficients[3]:.8f}$'\n",
    "print(function_text)\n",
    "#ax2.plot(xaxis, curve_y, color='grey', label=function_text)\n",
    "\n",
    "# Initial guess for parameters\n",
    "initial_guess = [1, 1]\n",
    "\n",
    "# Perform robust fitting using least squares with Huber loss\n",
    "result = least_squares(residuals, initial_guess, args=(xdata, ydata), loss='huber')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "\n",
    "\n",
    "# Perform the minimization using absolute error\n",
    "result = minimize(objective_function, initial_guess, args=(xdata, ydata), method='Nelder-Mead')\n",
    "\n",
    "# Extract fitted parameters\n",
    "a, b = result.x\n",
    "\n",
    "# Create an x-axis for plotting the fitted curve\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "\n",
    "# Evaluate the fitted logarithmic function\n",
    "curve_y = logarithmic_function(xaxis, a, b)\n",
    "ax2.plot(xaxis, curve_y, color='grey')\n",
    "ax2.scatter(xdata,  ydata, color='k', marker='*', alpha=0.7)\n",
    "ax2.scatter(xdata[4],  ydata[4], color='r', marker='*', alpha=1)\n",
    "ax2.scatter(1315*445.903*24,0.18, color='b', marker='+',s=93, alpha=1) #SYN: 1315, 445.903*24*1315\n",
    "\n",
    "ax2.tick_params(axis='x', rotation=45)\n",
    "\n",
    "ax1.set_title('e) Psi uncertainty (days obs.)', loc='left')\n",
    "ax2.set_title('f) Psi uncertainty (neutrons obs.)', loc='left')\n",
    "\n",
    "ax2.set_xlabel('Total neutron counts', loc='left')\n",
    "#ax2.legend()\n",
    "#ax1.set_title('Uncertainty', loc='left')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "pngname=\"./plots/Figure_5_Moderated_uncertainty_all.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9 (py_env_lda)",
   "language": "python",
   "name": "py_env_lda"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}