Figure_4_sensitivity.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b064147",
   "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",
    "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,
   "id": "e45a361d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "current_dir = os.getcwd()\n",
    "lib_dir = os.path.join(current_dir, 'lib')\n",
    "sys.path.append(lib_dir)\n",
    "\n",
    "import numpy as np\n",
    "import CRNS_library as CRNS_lib\n",
    "import time\n",
    "import scipy.io as sio\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import CRNS_optimizer as CRNS_opt\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "\n",
    "def func(x, p1, p2):\n",
    "    return p1 * np.cos(p2 * x) + p2 * np.sin(p1 * x)\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d75c1888",
   "metadata": {},
   "outputs": [],
   "source": [
    "src_path = './data'\n",
    "##Read in data for one site as example:\n",
    "#read incoming cosmic ray intensity\n",
    "file_path=src_path+'/Scaling/JUNG1_Data_UTC.TXT'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Incoming'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Jung = new_data.set_index('DateTime')\n",
    "\n",
    "src_path = 'D:/Projekte_ongoing/CRNS_ML/'\n",
    "##Read in data for one site as example:\n",
    "#read incoming cosmic ray intensity\n",
    "file_path=src_path+'/Scaling/OULU_1h_UTC.TXT'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Oulu'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Oulu = new_data.set_index('DateTime')\n",
    "\n",
    "\n",
    "src_path = 'D:/Projekte_ongoing/CRNS_ML/'\n",
    "##Read in data for one site as example:\n",
    "#read incoming cosmic ray intensity\n",
    "file_path=src_path+'/Scaling/APTY_1h_UTC.TXT'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Apty'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Apty = new_data.set_index('DateTime')\n",
    "\n",
    "src_path = 'D:/Projekte_ongoing/CRNS_ML/'\n",
    "##Read in data for one site as example:\n",
    "#read incoming cosmic ray intensity\n",
    "file_path=src_path+'/Scaling/AthensR8.53Alt260 m.TXT'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Athens'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Athens = new_data.set_index('DateTime')\n",
    "\n",
    "\n",
    "src_path = 'D:/Projekte_ongoing/CRNS_ML/'\n",
    "##Read in data for one site as example:\n",
    "#read incoming cosmic ray intensity\n",
    "file_path=src_path+'/Scaling/MXCO_1h_UTC8.28Alt2274m.TXT'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Mxco'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Mxco = new_data.set_index('DateTime')\n",
    "\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd468eb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#ID: 33 Reading: MEBCK001\n",
    "MEBCK001_site_ID = 33\n",
    "ii=33\n",
    "nn=\"MEBCK001\"\n",
    "bd = 1.31\n",
    "#for ii, nn in enumerate(sites_names[59:-2], start=59):\n",
    "print(\"ID:\",ii, \"Reading:\",nn)\n",
    "No = 1205\n",
    "\n",
    "psi = 0.7\n",
    "omega = 0.0054\n",
    "beta_ref = -0.0074\n",
    "beta = beta_ref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fefe8fd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta=beta_ref\n",
    "\n",
    "## Schäfertal ID 3 removed; and columns renamed; asterisk remarks removed; BUC001 data gaps were removed\n",
    "fname='EU_Study_sites.csv'\n",
    "fname=src_path+\"/Scaling/COSMOS_Europe_Data/Additional_information.csv\"\n",
    "site_info = pd.read_csv(fname)\n",
    "fname=src_path+\"/Scaling/COSMOS_Europe_Data/General_information.csv\"\n",
    "site_general_info = pd.read_csv(fname)\n",
    "sites_names = site_info.StationID.to_list()\n",
    "CutoffRigidity = site_info.CutoffRigidity.to_list()\n",
    "sites_Bulkdensity = site_info.Bulkdensity.to_list()\n",
    "sites_Altitude = site_general_info.Altitude.to_list()\n",
    "sites_Precip = site_general_info.Meanannualprecipitation.to_list()\n",
    "sites_Temperature = site_general_info.Meanairtemperature.to_list()\n",
    "sites_lat = site_general_info.Latitude.to_list()\n",
    "dfi = pd.DataFrame({\n",
    "    'StationID': sites_names,\n",
    "    'CutoffRigidity': CutoffRigidity,\n",
    "    'sites_Precip': sites_Precip,\n",
    "    'sites_Temperature': sites_Temperature,\n",
    "    'Bulkdensity': sites_Bulkdensity,\n",
    "    'Altitude': sites_Altitude\n",
    "})\n",
    "\n",
    "# Add new columns (beta, omega, psi) with initial values (you can update these later)\n",
    "dfi['beta'] = np.nan\n",
    "dfi['beta_sd'] = np.nan\n",
    "dfi['beta_mae'] = np.nan\n",
    "dfi['beta_rmse'] = np.nan\n",
    "dfi['beta_oulu'] = np.nan\n",
    "dfi['beta_apty'] = np.nan\n",
    "dfi['beta_mxco'] = np.nan\n",
    "dfi['beta_athens'] = np.nan\n",
    "dfi['omega'] = np.nan\n",
    "dfi['omega_sd'] = np.nan\n",
    "dfi['omega_oulu'] = np.nan\n",
    "dfi['omega_apty'] = np.nan\n",
    "dfi['omega_mxco'] = np.nan\n",
    "dfi['omega_athens'] = np.nan\n",
    "dfi['omega_mae'] = np.nan\n",
    "dfi['omega_rmse'] = np.nan\n",
    "dfi['psi'] = np.nan\n",
    "dfi['psi_sd'] = np.nan\n",
    "dfi['psi_oulu'] = np.nan\n",
    "dfi['psi_apty'] = np.nan\n",
    "dfi['psi_athens'] = np.nan\n",
    "dfi['psi_mxco'] = np.nan\n",
    "dfi['psi_mae'] = np.nan\n",
    "dfi['psi_rmse'] = np.nan\n",
    "dfi['BetaDesilets'] = np.nan\n",
    "dfi['mean_P'] = np.nan\n",
    "\n",
    "DF = CRNS_lib.load_EU_data(src_path+'/Scaling/COSMOS_Europe_Data/',nn)\n",
    "DF['Abs_h'] = CRNS_opt.Absolute_conv(DF.RH, DF.TEMP)\n",
    "#resample to 1 hour exactly, for merging with 1h inc. radiation\n",
    "DF = DF.resample('1H').mean()\n",
    "\n",
    "DF = DF.join(Jung)\n",
    "DF = DF.join(Apty)\n",
    "DF = DF.join(Oulu)\n",
    "DF = DF.join(Athens)\n",
    "DF = DF.join(Mxco)\n",
    "\n",
    "#calculate relative scaling for inc. radiation\n",
    "DF[\"Incoming\"] = DF[\"Incoming\"]/np.nanmean(DF[\"Incoming\"])#Jung\n",
    "DF[\"Apty\"] = DF[\"Apty\"]/np.nanmean(DF[\"Apty\"])\n",
    "DF[\"Oulu\"] = DF[\"Oulu\"]/np.nanmean(DF[\"Oulu\"])\n",
    "DF[\"Mxco\"] = DF[\"Mxco\"]/np.nanmean(DF[\"Mxco\"])\n",
    "DF[\"Athens\"] = DF[\"Athens\"]/np.nanmean(DF[\"Athens\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7178cf6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f72ee56",
   "metadata": {},
   "outputs": [],
   "source": [
    "#counts when SWC is observed after correction\n",
    "beta=beta_ref\n",
    "#Npih = (0.0808/(SWC / bd + 0.115) + 0.372 )* No"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5210b17",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Constants\n",
    "\n",
    "pref = 1011  # Reference Atmospheric pressure in hPa\n",
    "ense_size = 500\n",
    "# Define ranges for P and beta\n",
    "P_range = np.linspace(963, 1043, ense_size)\n",
    "beta_range = np.linspace(-0.0063, -0.0078, ense_size)\n",
    "\n",
    "# Define ranges for Inc and psi\n",
    "# Calculate reference Npih and SWC\n",
    "P_ref = pref  # Reference Atmospheric pressure\n",
    "SWC_ref = 0.3  # Reference SWC value\n",
    "cor_ref = CRNS_opt.corr_exp_pressure(P_ref, pref, beta_ref)\n",
    "Npih_ref= (0.0808/(SWC_ref / bd + 0.115) + 0.372 )* No\n",
    "print(Npih_ref)\n",
    "print(CRNS_opt.Calculate_swc_from_Npih(Npih_ref, No, bd))\n",
    "\n",
    "\n",
    "#  Min = 963.6\n",
    "#  Max = 1043.6\n",
    "#  Standard Deviation = 8.980325633854576\n",
    "#  95th Percentile = 1026.1\n",
    "\n",
    "# Initialize a matrix to store SWC differences\n",
    "SWC_diff = np.zeros((len(P_range), len(beta_range)))\n",
    "\n",
    "# Calculate SWC for each combination of P and beta\n",
    "for i, P in enumerate(P_range):\n",
    "    for j, beta in enumerate(beta_range):\n",
    "        ##Calculate observed neutron flux using beta:\n",
    "        cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "        N_ij = Npih_ref / cor_coeff\n",
    "        ## Two SWC will be possibly estimated.\n",
    "        # a) true SWC (0.1, 0.2 or 0.3)\n",
    "        # b) SWC_estimated with wrong beta using wrong cor_coeff(P)\n",
    "        #Only store plot the difference\n",
    "        SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij*cor_coeff, No, bd)\n",
    "        #Calculate SWC with real beta value:\n",
    "        cor_coeff_case = CRNS_opt.corr_exp_pressure(P, pref, beta_ref)\n",
    "        SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij*cor_coeff_case, No, bd)\n",
    "        SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "# Plot the heatmap\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.contourf(beta_range, P_range, SWC_diff, levels=50, cmap='viridis')\n",
    "plt.colorbar(label='SWC Difference')\n",
    "plt.xlabel('Beta')\n",
    "plt.ylabel('Atmospheric Pressure (hPa)')\n",
    "plt.title('Heatmap of SWC Difference to Reference Value')\n",
    "plt.show()\n",
    "\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65772c97",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d87b180",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize subplots\n",
    "fig, axes = plt.subplots(1, 3, figsize=(24, 8))\n",
    "\n",
    "# Reference SWC values\n",
    "SWC_refs = [0.1, 0.2, 0.3]\n",
    "titles = ['SWC_ref = 0.1', 'SWC_ref = 0.2', 'SWC_ref = 0.3']\n",
    "\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((len(P_range), len(beta_range)))\n",
    "\n",
    "    # Calculate SWC for each combination of P and beta\n",
    "    for i, P in enumerate(P_range):\n",
    "        for j, beta in enumerate(beta_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            ## Two SWC will be possibly estimated.\n",
    "            # a) true SWC (0.1, 0.2 or 0.3)\n",
    "            # b) SWC_estimated with wrong beta using wrong cor_coeff(P)\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_exp_pressure(P, pref, beta_ref)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "            #cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[idx]\n",
    "    c = ax.contourf(beta_range, P_range, SWC_diff, levels=50, cmap='RdYlGn_r')\n",
    "    fig.colorbar(c, ax=ax, label='SWC Difference')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(beta_range, P_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['blue','k', 'red'], linewidths=2)\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=8)\n",
    "\n",
    "    ax.set_xlabel('Beta')\n",
    "    ax.set_ylabel('Atmospheric Pressure (hPa)')\n",
    "    ax.set_title(titles[idx])\n",
    "\n",
    "plt.suptitle('Heatmap of SWC Difference to Reference Values')\n",
    "plt.tight_layout(rect=[0, 0, 1, 0.96])\n",
    "plt.show()\n",
    "\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "269a2b62",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "Inc_ref = 1.0  # Reference incoming cosmic ray intensity\n",
    "psi_ref = 0.7\n",
    "\n",
    "# Define ranges for Inc and psi\n",
    "Inc_range = np.linspace(0.92, 1.08, ense_size)\n",
    "\n",
    "psi_range = np.linspace(0.2, 1.0, ense_size)\n",
    "\n",
    "# Initialize subplots\n",
    "fig, axes = plt.subplots(2, 4, figsize=(24, 12))\n",
    "\n",
    "# Reference SWC values\n",
    "SWC_refs = [0.1, 0.2, 0.3, 0.4]\n",
    "titles = ['SWC_ref = 0.1', 'SWC_ref = 0.2', 'SWC_ref = 0.3','SWC_ref = 0.4']\n",
    "\n",
    "# Plot for beta and P\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((len(P_range), len(beta_range)))\n",
    "\n",
    "    # Calculate SWC for each combination of P and beta\n",
    "    for i, P in enumerate(P_range):\n",
    "        for j, beta in enumerate(beta_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            ## Two SWC will be possibly estimated.\n",
    "            # a) true SWC (0.1, 0.2 or 0.3)\n",
    "            # b) SWC_estimated with wrong beta using wrong cor_coeff(P)\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_exp_pressure(P, pref, beta_ref)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "            #cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "            #cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[0, idx]\n",
    "    c = ax.contourf(beta_range, P_range, SWC_diff, levels=50, cmap='RdYlGn_r')\n",
    "    fig.colorbar(c, ax=ax, label='SWC Difference')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(beta_range, P_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['blue','k', 'red'], linewidths=2)\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "\n",
    "    ax.set_xlabel('Beta')\n",
    "    ax.set_ylabel('Atmospheric Pressure (hPa)')\n",
    "    ax.set_title(titles[idx])\n",
    "\n",
    "# Plot for Inc and psi\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((len(Inc_range), len(psi_range)))\n",
    "\n",
    "    # Calculate SWC for each combination of Inc and psi\n",
    "    for i, Inc in enumerate(Inc_range):\n",
    "        for j, psi in enumerate(psi_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi_ref)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            ## Two SWC will be possibly estimated.\n",
    "            # a) true SWC (0.1, 0.2 or 0.3)\n",
    "            # b) SWC_estimated with wrong psi using wrong cor_coeff(P)\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "            \n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[1, idx]\n",
    "    c = ax.contourf(psi_range, Inc_range, SWC_diff, levels=50, cmap='RdYlGn_r')\n",
    "    fig.colorbar(c, ax=ax, label='SWC Difference')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(psi_range, Inc_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['blue','k', 'red'], linewidths=2)\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "\n",
    "    ax.set_xlabel('Psi')\n",
    "    ax.set_ylabel('Incoming Cosmic Ray Intensity')\n",
    "    ax.set_title(titles[idx])\n",
    "\n",
    "plt.suptitle('Heatmap of SWC Difference to Reference Values')\n",
    "plt.tight_layout(rect=[0, 0, 1, 0.96])\n",
    "plt.show()\n",
    "\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n",
    "\n",
    "\n",
    "#pngname=\"./sensitivity.png\"\n",
    "#print(\"Save file as \"+pngname)\n",
    "#fig.savefig(pngname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b4bbf0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "Inc_ref = 1.0  # Reference incoming cosmic ray intensity\n",
    "\n",
    "psi_ref = 0.68\n",
    "omega_ref = -0.0054\n",
    "abs_ref = 7.0\n",
    "ense_size = 500\n",
    "\n",
    "# Define ranges for Inc and psi\n",
    "Inc_range = np.linspace(0.92, 1.08, ense_size)\n",
    "abs_range = np.linspace(1, 20, ense_size)\n",
    "\n",
    "psi_range = np.linspace(0.2, 1.0, ense_size)\n",
    "omega_range = np.linspace(-0.003, -0.01, ense_size)\n",
    "\n",
    "# Initialize subplots\n",
    "fig, axes = plt.subplots(3, 4, figsize=(24, 16))\n",
    "\n",
    "# Reference SWC values\n",
    "SWC_refs = [0.1, 0.2, 0.3, 0.4]\n",
    "titles = ['SWC_ref = 0.1', 'SWC_ref = 0.2', 'SWC_ref = 0.3','SWC_ref = 0.4']\n",
    "\n",
    "# Plot for beta and P\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((len(P_range), len(beta_range)))\n",
    "\n",
    "    # Calculate SWC for each combination of P and beta\n",
    "    for i, P in enumerate(P_range):\n",
    "        for j, beta in enumerate(beta_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            ## Two SWC will be possibly estimated.\n",
    "            # a) true SWC (0.1, 0.2 or 0.3)\n",
    "            # b) SWC_estimated with wrong beta using wrong cor_coeff(P)\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_exp_pressure(P, pref, beta_ref)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "            #cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "            #cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[0, idx]\n",
    "    c = ax.contourf(beta_range, P_range, SWC_diff, levels=50, cmap='RdYlGn_r')\n",
    "    fig.colorbar(c, ax=ax, label='SWC Difference')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(beta_range, P_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['blue','k', 'red'], linewidths=2)\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "\n",
    "    ax.set_xlabel('Beta')\n",
    "    ax.set_ylabel('Atmospheric Pressure (hPa)')\n",
    "    ax.set_title(titles[idx])\n",
    "\n",
    "# Plot for Inc and psi\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((ense_size,ense_size))\n",
    "\n",
    "    # Calculate SWC for each combination of Inc and psi\n",
    "    for i, Inc in enumerate(Inc_range):\n",
    "        for j, psi in enumerate(psi_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi_ref)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            ## Two SWC will be possibly estimated.\n",
    "            # a) true SWC (0.1, 0.2 or 0.3)\n",
    "            # b) SWC_estimated with wrong psi using wrong cor_coeff(P)\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[1, idx]\n",
    "    c = ax.contourf(psi_range, Inc_range, SWC_diff, levels=50, cmap='RdYlGn_r')\n",
    "    fig.colorbar(c, ax=ax, label='SWC Difference')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(psi_range, Inc_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['blue','k', 'red'], linewidths=2)\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "\n",
    "    ax.set_xlabel('Psi')\n",
    "    ax.set_ylabel('Incoming Cosmic Ray Intensity')\n",
    "    ax.set_title(titles[idx])\n",
    "\n",
    "# Plot for Inc and psi\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((ense_size,ense_size))\n",
    "\n",
    "    # Calculate SWC for each combination of Inc and psi\n",
    "    for i, abs_val in enumerate(abs_range):\n",
    "        for j, omega in enumerate(omega_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_lin_humidity(abs_val, abs_ref, omega_ref)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            ## Two SWC will be possibly estimated.\n",
    "            # a) true SWC (0.1, 0.2 or 0.3)\n",
    "            # b) SWC_estimated with wrong psi using wrong cor_coeff(P)\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_lin_humidity(abs_val, abs_ref, omega)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "            \n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[2, idx]\n",
    "    c = ax.contourf(omega_range, abs_range, SWC_diff, levels=50, cmap='RdYlGn_r')\n",
    "    fig.colorbar(c, ax=ax, label='SWC Difference')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(omega_range, abs_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['blue','k', 'red'], linewidths=2)\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "\n",
    "    ax.set_xlabel('Omega')\n",
    "    ax.set_ylabel('Air Humidity')\n",
    "    ax.set_title(titles[idx])\n",
    "\n",
    "plt.suptitle('Heatmap of SWC Difference to Reference Values')\n",
    "plt.tight_layout(rect=[0, 0, 1, 0.96])\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n",
    "\n",
    "#pngname=\"./sensitivityA.png\"\n",
    "#print(\"Save file as \"+pngname)\n",
    "#fig.savefig(pngname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "991cfaf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "Inc_ref = 1.0  # Reference incoming cosmic ray intensity\n",
    "\n",
    "psi_ref = 0.68\n",
    "omega_ref = -0.0054\n",
    "abs_ref = 7.0\n",
    "ense_size = 500\n",
    "\n",
    "# Define ranges for Inc and psi\n",
    "Inc_range = np.linspace(0.92, 1.08, ense_size)\n",
    "abs_range = np.linspace(1, 20, ense_size)\n",
    "\n",
    "psi_range = np.linspace(0.2, 1.0, ense_size)\n",
    "omega_range = np.linspace(-0.003, -0.01, ense_size)\n",
    "\n",
    "# Initialize subplots\n",
    "#fig, axes = plt.subplots(3, 4, figsize=(24, 16))\n",
    "#SWC_refs = [0.1, 0.2, 0.3, 0.4]\n",
    "#titles = ['SWC_ref = 0.1', 'SWC_ref = 0.2', 'SWC_ref = 0.3','SWC_ref = 0.4']\n",
    "\n",
    "# Initialize subplots\n",
    "fig, axes = plt.subplots(3, 3, figsize=(15, 12))\n",
    "SWC_refs = [0.1, 0.25,  0.4]\n",
    "titles = ['Ref. Soil Moisture = 0.1 m³/m³', 'Ref. Soil Moisture = 0.25 m³/m³', 'Ref. Soil Moisture  = 0.4 m³/m³']\n",
    "\n",
    "# Plot for beta and P\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((len(P_range), len(beta_range)))\n",
    "\n",
    "    # Calculate SWC for each combination of P and beta\n",
    "    for i, P in enumerate(P_range):\n",
    "        for j, beta in enumerate(beta_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            ## Two SWC will be possibly estimated.\n",
    "            # a) true SWC (0.1, 0.2 or 0.3)\n",
    "            # b) SWC_estimated with wrong beta using wrong cor_coeff(P)\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_exp_pressure(P, pref, beta_ref)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "            #cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "            #cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[0, idx]\n",
    "    c = ax.contourf(beta_range, P_range, SWC_diff, levels=50, cmap='RdYlBu_r')\n",
    "    fig.colorbar(c, ax=ax, label='Difference (soil moisture)')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(beta_range, P_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['k','k', 'k'], linewidths=[1,1,1], linestyle=['--','-','--'])\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "    ax.tick_params(axis='x', rotation=45)\n",
    "\n",
    "    ax.set_xlabel('Beta (atmospheric pressure)')\n",
    "    ax.set_ylabel('Atmospheric pressure (hPa)')\n",
    "    ax.set_title(titles[idx])\n",
    "\n",
    "# Plot for Inc and psi\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((ense_size,ense_size))\n",
    "\n",
    "    # Calculate SWC for each combination of Inc and psi\n",
    "    for i, Inc in enumerate(Inc_range):\n",
    "        for j, psi in enumerate(psi_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi_ref)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            ## Two SWC will be possibly estimated.\n",
    "            # a) true SWC (0.1, 0.2 or 0.3)\n",
    "            # b) SWC_estimated with wrong psi using wrong cor_coeff(P)\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[1, idx]\n",
    "    c = ax.contourf(psi_range, Inc_range, SWC_diff, levels=50, cmap='RdYlBu_r')\n",
    "    fig.colorbar(c, ax=ax, label='Difference (soil moisture)')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(psi_range, Inc_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['k','k', 'k'], linewidths=[1,1,1], linestyle=['--','-','--'])\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "    ax.tick_params(axis='x', rotation=45)\n",
    "\n",
    "    ax.set_xlabel('Psi (incoming)')\n",
    "    ax.set_ylabel('Incoming cosmic ray intensity')\n",
    "    ax.set_title(titles[idx])\n",
    "\n",
    "# Plot for Inc and psi\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((ense_size,ense_size))\n",
    "\n",
    "    # Calculate SWC for each combination of Inc and psi\n",
    "    for i, abs_val in enumerate(abs_range):\n",
    "        for j, omega in enumerate(omega_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_lin_humidity(abs_val, abs_ref, omega_ref)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            ## Two SWC will be possibly estimated.\n",
    "            # a) true SWC (0.1, 0.2 or 0.3)\n",
    "            # b) SWC_estimated with wrong psi using wrong cor_coeff(P)\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_lin_humidity(abs_val, abs_ref, omega)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "            \n",
    "            #SWC = CRNS_opt.Calculate_swc_from_Npih(Npih_ref * cor_coeff, No, bd)\n",
    "            #SWC_diff[i, j] = SWC - SWC_ref\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[2, idx]\n",
    "    c = ax.contourf(omega_range, abs_range, SWC_diff, levels=50, cmap='RdYlBu_r')\n",
    "    fig.colorbar(c, ax=ax, label='Difference (soil moisture)')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(omega_range, abs_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['k','k', 'k'], linewidths=[1,1,1], linestyle=['--','-','--'])\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "    ax.tick_params(axis='x', rotation=45)\n",
    "    ax.set_xlabel('Omega (air humidity)')\n",
    "    ax.set_ylabel('Air humidity')\n",
    "    ax.set_title(titles[idx])\n",
    "\n",
    "plt.tight_layout(rect=[0, 0, 1, 0.96])\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "pngname=\"./sensitivityb.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b5b2d6d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f687b39",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to print ranges, standard deviation, and 95th percentiles of individual columns\n",
    "def print_column_statistics(df):\n",
    "    for column in df.columns:\n",
    "        min_val = df[column].min()\n",
    "        max_val = df[column].max()\n",
    "        std_dev = df[column].std()\n",
    "        percentile_05 = df[column].quantile(0.05)\n",
    "        percentile_50 = df[column].quantile(0.50)\n",
    "        percentile_95 = df[column].quantile(0.95)\n",
    "        print(f\"Column {column}:\")\n",
    "        print(f\"  Min = {min_val}\")\n",
    "        print(f\"  Max = {max_val}\")\n",
    "        print(f\"   5th Percentile = {percentile_05}\")\n",
    "        print(f\"  50th Percentile = {percentile_50}\")\n",
    "        print(f\"  95th Percentile = {percentile_95}\")\n",
    "        print()\n",
    "\n",
    "# Call the function\n",
    "print_column_statistics(DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d1f9171",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(beta)\n",
    "\n",
    "print(DF[\"Mxco\"])\n",
    "print(max(DF[\"Mxco\"]))\n",
    "plt.plot(DF[\"Mxco\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78053a79",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "baf08897",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cb524e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_prefix = ['a) ', 'b) ', 'c) ', 'd) ', 'e) ', 'f) ', 'g) ', 'e) ', 'h) ']\n",
    "pref_ct = 0\n",
    "Inc_ref = 1.0  # Reference incoming cosmic ray intensity\n",
    "\n",
    "psi_ref = 0.68\n",
    "omega_ref = -0.0054\n",
    "abs_ref = 7.0\n",
    "ense_size = 500\n",
    "\n",
    "# Define ranges for Inc and psi\n",
    "Inc_range = np.linspace(0.92, 1.08, ense_size)\n",
    "abs_range = np.linspace(1, 20, ense_size)\n",
    "\n",
    "psi_range = np.linspace(0.2, 1.0, ense_size)\n",
    "omega_range = np.linspace(-0.003, -0.01, ense_size)\n",
    "\n",
    "# Initialize subplots\n",
    "#fig, axes = plt.subplots(3, 4, figsize=(24, 16))\n",
    "#SWC_refs = [0.1, 0.2, 0.3, 0.4]\n",
    "#titles = ['SWC_ref = 0.1', 'SWC_ref = 0.2', 'SWC_ref = 0.3','SWC_ref = 0.4']\n",
    "\n",
    "# Initialize subplots\n",
    "fig, axes = plt.subplots(3, 3, figsize=(15, 12))\n",
    "SWC_refs = [0.1, 0.25,  0.4]\n",
    "titles = ['True soil moisture = 0.1 m³/m³', 'True soil moisture = 0.25 m³/m³', 'True soil moisture  = 0.4 m³/m³']\n",
    "\n",
    "global_min = float('inf')\n",
    "global_max = float('-inf')\n",
    "\n",
    "# Plot for beta and P\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((len(P_range), len(beta_range)))\n",
    "\n",
    "    # Calculate SWC for each combination of P and beta\n",
    "    for i, P in enumerate(P_range):\n",
    "        for j, beta in enumerate(beta_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_exp_pressure(P, pref, beta_ref)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "    # Update global min and max values\n",
    "    global_min = min(global_min, SWC_diff.min())\n",
    "    global_max = max(global_max, SWC_diff.max())\n",
    "\n",
    "# Plot for Inc and psi\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((ense_size,ense_size))\n",
    "\n",
    "    # Calculate SWC for each combination of Inc and psi\n",
    "    for i, Inc in enumerate(Inc_range):\n",
    "        for j, psi in enumerate(psi_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi_ref)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "    # Update global min and max values\n",
    "    global_min = min(global_min, SWC_diff.min())\n",
    "    global_max = max(global_max, SWC_diff.max())\n",
    "\n",
    "\n",
    "# Plot for Inc and psi\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((ense_size,ense_size))\n",
    "\n",
    "    # Calculate SWC for each combination of Inc and psi\n",
    "    for i, abs_val in enumerate(abs_range):\n",
    "        for j, omega in enumerate(omega_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_lin_humidity(abs_val, abs_ref, omega_ref)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_lin_humidity(abs_val, abs_ref, omega)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "    # Update global min and max values\n",
    "    global_min = min(global_min, SWC_diff.min())\n",
    "    global_max = max(global_max, SWC_diff.max())\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "# Plot for beta and P\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((len(P_range), len(beta_range)))\n",
    "\n",
    "    # Calculate SWC for each combination of P and beta\n",
    "    for i, P in enumerate(P_range):\n",
    "        for j, beta in enumerate(beta_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_exp_pressure(P, pref, beta)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_exp_pressure(P, pref, beta_ref)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[0, idx]\n",
    "    c = ax.contourf(beta_range, P_range, SWC_diff, levels=np.linspace(global_min, global_max, 50), cmap='RdYlBu_r')\n",
    "    fig.colorbar(c, ax=ax, label='Difference (soil moisture)')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(beta_range, P_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['k','k', 'k'], linewidths=[1,1,1], linestyle=['--','-','--'])\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "    ax.tick_params(axis='x', rotation=45)\n",
    "\n",
    "    ax.set_xlabel('Beta (atmospheric pressure)')\n",
    "    ax.set_ylabel('Atmospheric pressure (hPa)')\n",
    "    ax.set_title(my_prefix[pref_ct]+titles[idx],loc=\"left\")\n",
    "    pref_ct = pref_ct + 1\n",
    "    \n",
    "# Plot for Inc and psi\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((ense_size,ense_size))\n",
    "\n",
    "    # Calculate SWC for each combination of Inc and psi\n",
    "    for i, Inc in enumerate(Inc_range):\n",
    "        for j, psi in enumerate(psi_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi_ref)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_lin_incoming_correction(Inc, Inc_ref, psi)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[1, idx]\n",
    "    c = ax.contourf(psi_range, Inc_range, SWC_diff, levels=np.linspace(global_min, global_max, 50), cmap='RdYlBu_r')\n",
    "    fig.colorbar(c, ax=ax, label='Difference (soil moisture)')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(psi_range, Inc_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['k','k', 'k'], linewidths=[1,1,1], linestyle=['--','-','--'])\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "    ax.tick_params(axis='x', rotation=45)\n",
    "\n",
    "    ax.set_xlabel('Psi (incoming)')\n",
    "    ax.set_ylabel('Incoming cosmic ray intensity')\n",
    "    ax.set_title(my_prefix[pref_ct]+titles[idx],loc=\"left\")\n",
    "    pref_ct = pref_ct + 1\n",
    "    \n",
    "# Plot for Inc and psi\n",
    "for idx, SWC_ref in enumerate(SWC_refs):\n",
    "    # Calculate reference Npih\n",
    "    Npih_ref = (0.0808 / (SWC_ref / bd + 0.115) + 0.372) * No\n",
    "    \n",
    "    # Initialize a matrix to store SWC differences\n",
    "    SWC_diff = np.zeros((ense_size,ense_size))\n",
    "\n",
    "    # Calculate SWC for each combination of Inc and psi\n",
    "    for i, abs_val in enumerate(abs_range):\n",
    "        for j, omega in enumerate(omega_range):\n",
    "            ##Calculate observed neutron flux using beta:\n",
    "            cor_coeff = CRNS_opt.corr_lin_humidity(abs_val, abs_ref, omega_ref)\n",
    "            N_ij = Npih_ref / cor_coeff\n",
    "            #Only store plot the difference\n",
    "            SWC_true = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff, No, bd)\n",
    "            #Calculate SWC with real beta value:\n",
    "            cor_coeff_case = CRNS_opt.corr_lin_humidity(abs_val, abs_ref, omega)\n",
    "            SWC_est = CRNS_opt.Calculate_swc_from_Npih(N_ij * cor_coeff_case, No, bd)\n",
    "            SWC_diff[i, j] = SWC_est - SWC_true\n",
    "\n",
    "    # Plot the heatmap\n",
    "    ax = axes[2, idx]\n",
    "    c = ax.contourf(omega_range, abs_range, SWC_diff, levels=np.linspace(global_min, global_max, 50), cmap='RdYlBu_r')\n",
    "    fig.colorbar(c, ax=ax, label='Difference (soil moisture)')\n",
    "    \n",
    "    # Add contour lines at 0.02, 0.0 and -0.02\n",
    "    contour_lines = ax.contour(omega_range, abs_range, SWC_diff, levels=[-0.02, 0.0, 0.02], colors=['k','k', 'k'], linewidths=[1,1,1], linestyle=['--','-','--'])\n",
    "    ax.clabel(contour_lines, fmt={-0.02: '-0.02',0.00: '0.0', 0.02: '0.02'}, inline=True, fontsize=12)\n",
    "    ax.tick_params(axis='x', rotation=45)\n",
    "    ax.set_xlabel('Omega (air humidity)')\n",
    "    ax.set_ylabel('Air humidity')\n",
    "    ax.set_title(my_prefix[pref_ct]+titles[idx],loc=\"left\")\n",
    "    pref_ct = pref_ct + 1\n",
    "    \n",
    "plt.tight_layout(rect=[0, 0, 1, 0.96])\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "pngname=\"./Figure_4_sensitivity_b1.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb007100",
   "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": 5
}