{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "90058a59-2941-4598-b7cf-9e6fe298e6ba",
   "metadata": {},
   "source": [
    "# Fit NASA9 Polynomials to Gibbs Free Energy (G vs T)\n",
    "\n",
    "This notebook reads a CSV of Gibbs free energies (J/mol) vs temperature (K)  \n",
    "and fits **NASA9-style polynomial coefficients** separately on a low-T and a high-T  \n",
    "range split at `T_split` (default: 500 K). It then:\n",
    "\n",
    "1. Saves the two sets of 9 coefficients per species to `nasa9_coefficients.txt`  \n",
    "2. Plots observed G vs T alongside the fitted curves  \n",
    "3. Plots residuals (fit − observed)\n",
    "\n",
    "---\n",
    "\n",
    "### Input data\n",
    "The file `Gibbs_Energies_All.csv` contains:\n",
    "- A one-line comment  \n",
    "- A comma-only “blank” line  \n",
    "- A header row  \n",
    "- A temperature column labeled `\"T (K)\"` (or similar; the loader is robust)  \n",
    "- One column per species with **Gibbs free energy values** in J/mol  \n",
    "\n",
    "---\n",
    "\n",
    "### Model being fitted\n",
    "We use the NASA9 functional form for **enthalpy (H)** and **entropy (S)** as functions of T  \n",
    "(with R = 8.3144621 J/(mol·K)), and compute:\n",
    "\n",
    "\\[\n",
    "G(T) = H(T) - TS(T)\n",
    "\\]\n",
    "\n",
    "We fit the 9 coefficients `a0..a8` that define H(T) and S(T).  \n",
    "(Cp equations are omitted here.)\n",
    "\n",
    "---\n",
    "\n",
    "### Output\n",
    "- `nasa9_coefficients.txt`: For each species, two lines of 9 coefficients:  \n",
    "  - Low-temperature fit (`T ≤ T_split`)  \n",
    "  - High-temperature fit (`T > T_split`)  \n",
    "- One figure per species with G-vs-T and residual plots.\n",
    "\n",
    "---\n",
    "\n",
    "### Usage\n",
    "- Place `Gibbs_Energies_All.csv` in the same folder (or update the `csv_file` path).  \n",
    "- Run the final cell calling `run_all(...)`.\n",
    "\n",
    "---\n",
    "\n",
    "### Notes & Caveats\n",
    "- This fits **G(T)** via NASA9 H(T) and S(T); ensure your G values are consistent (units: J/mol).  \n",
    "- The split temperature `T_split` can be adjusted if needed (default 500 K).  \n",
    "- The loader is robust to slight variations in the temperature column name (e.g., `T (K)` vs `T(K)`).  \n",
    "- This is a numerical fit; inspect residuals to judge fit quality and adjust ranges/constraints if necessary.  \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5c8b2232-a322-483a-9465-891b8e11774f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting Water\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 85528 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting NH3\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 66528 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting Acetlyne/ethyne\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 73596 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting HCN\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 71672 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting Formaldehyde\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 77024 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting Adenine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 76652 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting Guanine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 69104 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting Cytosine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 72460 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting H2\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 67876 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting glycine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 71760 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting alanine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 77072 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting arginine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 77064 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting asparagine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 74924 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting aspartic acid\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 74204 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting cysteine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 70548 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting glutamic acid\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 74340 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting glutamine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 74740 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting histidine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 75036 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting isoleucine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 73384 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting leucine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 71068 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting lysine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 74240 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting methionine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 74876 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting phenylalanine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 74332 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting proline\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 71552 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting serine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 70632 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting threonine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 71840 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting tryptophan\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 76456 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting tyrosine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 76496 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting valine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 69452 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting beta-alanine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 82480 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting AAH\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 73620 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting 2AH\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 66648 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "→ Fitting Ethanamine\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 67760 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "All coefficients saved to nasa9_coefficients.txt\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import least_squares\n",
    "\n",
    "# Universal gas constant in J/(mol·K)\n",
    "R_GAS = 8.3144621\n",
    "\n",
    "\n",
    "def nasa9_properties(a, T):\n",
    "    \"\"\"\n",
    "    Compute H(T) and S(T) (per mole) from NASA9-style coefficients.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    a : array-like of shape (9,)\n",
    "        NASA9 coefficients [a0, ..., a8].\n",
    "    T : array-like\n",
    "        Temperatures in K.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    H : ndarray\n",
    "        Enthalpy in J/mol.\n",
    "    S : ndarray\n",
    "        Entropy in J/(mol·K).\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    This follows the NASA9 functional forms for H/(R·T) and S/R:\n",
    "        H/(R·T) = (-a0/T^2) + a1*ln(T)/T + a2 + (a3/2)T + (a4/3)T^2 + (a5/4)T^3 + (a6/5)T^4 + a7/T\n",
    "        S/R     = (-0.5*a0/T^2) - a1/T + a2*ln(T) + a3*T + (a4/2)T^2 + (a5/3)T^3 + (a6/4)T^4 + a8\n",
    "    Cp expressions are not used in this fitting workflow.\n",
    "    \"\"\"\n",
    "    a0, a1, a2, a3, a4, a5, a6, a7, a8 = a\n",
    "    T = np.asarray(T, dtype=float)\n",
    "\n",
    "    lnT    = np.log(T)\n",
    "    T_inv  = 1.0 / T\n",
    "    T_inv2 = T_inv * T_inv\n",
    "\n",
    "    # H/(R·T)\n",
    "    H_RT = (-a0 * T_inv2\n",
    "            + a1 * lnT * T_inv\n",
    "            + a2\n",
    "            + (a3 / 2.0) * T\n",
    "            + (a4 / 3.0) * T**2\n",
    "            + (a5 / 4.0) * T**3\n",
    "            + (a6 / 5.0) * T**4\n",
    "            + a7 * T_inv)\n",
    "\n",
    "    # S/R\n",
    "    S_R = (-0.5 * a0 * T_inv2\n",
    "           - a1 * T_inv\n",
    "           + a2 * lnT\n",
    "           + a3 * T\n",
    "           + (a4 / 2.0) * T**2\n",
    "           + (a5 / 3.0) * T**3\n",
    "           + (a6 / 4.0) * T**4\n",
    "           + a8)\n",
    "\n",
    "    # Convert to dimensional H and S\n",
    "    H = H_RT * R_GAS * T      # J/mol\n",
    "    S = S_R  * R_GAS          # J/(mol·K)\n",
    "    return H, S\n",
    "\n",
    "\n",
    "def compute_gibbs(a, T):\n",
    "    \"\"\"\n",
    "    Compute Gibbs free energy: G(T) = H(T) − T·S(T).\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    a : array-like of shape (9,)\n",
    "        NASA9 coefficients.\n",
    "    T : array-like\n",
    "        Temperatures in K.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    G : ndarray\n",
    "        Gibbs free energy in J/mol.\n",
    "    \"\"\"\n",
    "    H, S = nasa9_properties(a, T)\n",
    "    return H - T * S\n",
    "\n",
    "\n",
    "def residuals(a, T, G_obs):\n",
    "    \"\"\"\n",
    "    Residuals between model-predicted G(T) and observed G(T).\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    a : array-like of shape (9,)\n",
    "        Current guess for NASA9 coefficients.\n",
    "    T : array-like\n",
    "        Temperatures in K.\n",
    "    G_obs : array-like\n",
    "        Observed Gibbs free energy in J/mol.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    r : ndarray\n",
    "        Residuals = G_model(T) − G_obs(T).\n",
    "    \"\"\"\n",
    "    return compute_gibbs(a, T) - G_obs\n",
    "\n",
    "\n",
    "def fit_nasa9(T, G, T_split=500.0):\n",
    "    \"\"\"\n",
    "    Fit two NASA9 coefficient sets: one below/at T_split, one above T_split.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    T : ndarray\n",
    "        Temperatures in K.\n",
    "    G : ndarray\n",
    "        Observed Gibbs free energies (J/mol) at the same T points.\n",
    "    T_split : float, optional\n",
    "        Temperature (K) at which to split the fit (default 500.0).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    a_low : ndarray, shape (9,)\n",
    "        NASA9 coefficients for T ≤ T_split.\n",
    "    a_high : ndarray, shape (9,)\n",
    "        NASA9 coefficients for T > T_split.\n",
    "    \"\"\"\n",
    "    mask_low  = T <= T_split\n",
    "    mask_high = T >  T_split\n",
    "\n",
    "    # Initial guess: zeros (can be customized if needed)\n",
    "    a0 = np.zeros(9)\n",
    "\n",
    "    # Independent least-squares fits on each side of the split\n",
    "    res_low  = least_squares(residuals, a0, args=(T[mask_low],  G[mask_low]))\n",
    "    res_high = least_squares(residuals, a0, args=(T[mask_high], G[mask_high]))\n",
    "    return res_low.x, res_high.x\n",
    "\n",
    "\n",
    "def load_csv(csv_file):\n",
    "    \"\"\"\n",
    "    Load the Gibbs energy CSV, returning temperatures and a dict of species data.\n",
    "\n",
    "    Expected file layout:\n",
    "      - 1st line: comment\n",
    "      - 2nd line: comma-only placeholders (blank line)\n",
    "      - 3rd line: header with 'T (K)' (or similar) and species columns\n",
    "      - Data rows follow\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    csv_file : str or Path\n",
    "        Path to the CSV file (e.g., 'Gibbs_Energies_All.csv').\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    T : ndarray\n",
    "        Temperatures in K.\n",
    "    species : dict[str, ndarray]\n",
    "        Mapping of species name -> Gibbs free energy values (J/mol) at each T.\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    - The loader robustly identifies the temperature column by looking for\n",
    "      headers that normalize to 't(k)' or 't'.\n",
    "    - Columns that are empty or 'Unnamed: ...' placeholders are dropped.\n",
    "    \"\"\"\n",
    "    # Read the actual header on the 3rd line (skip comment + blank commas line)\n",
    "    df = pd.read_csv(csv_file, skiprows=2, header=0)\n",
    "\n",
    "    # Drop empty/placeholder columns\n",
    "    drop_cols = [c for c in df.columns if str(c).startswith(\"Unnamed\") or df[c].isna().all()]\n",
    "    df = df.drop(columns=drop_cols)\n",
    "\n",
    "    # Normalize column names (trim spaces)\n",
    "    df.columns = [c.strip() for c in df.columns]\n",
    "\n",
    "    # Identify the temperature column robustly\n",
    "    temp_candidates = [c for c in df.columns if c.replace(\" \", \"\").lower() in (\"t(k)\", \"t\")]\n",
    "    temp_col = temp_candidates[0] if temp_candidates else df.columns[0]\n",
    "\n",
    "    # Ensure numeric data\n",
    "    for c in df.columns:\n",
    "        df[c] = pd.to_numeric(df[c], errors=\"coerce\")\n",
    "    df = df.dropna(how=\"all\")\n",
    "\n",
    "    # Extract arrays\n",
    "    T = df[temp_col].to_numpy()\n",
    "    species = {c: df[c].to_numpy() for c in df.columns if c != temp_col}\n",
    "\n",
    "    # Safety: in case a temp-like column slipped into species (belt & suspenders)\n",
    "    for k in list(species.keys()):\n",
    "        if k.replace(\" \", \"\").lower() in (\"t(k)\", \"t\"):\n",
    "            species.pop(k)\n",
    "\n",
    "    return T, species\n",
    "\n",
    "\n",
    "def run_all(csv_file, T_split=500.0, output_file='nasa9_coefficients.txt'):\n",
    "    \"\"\"\n",
    "    Orchestrate the full workflow:\n",
    "      1) Load data\n",
    "      2) Fit NASA9 coefficients for each species on low/high ranges\n",
    "      3) Save coefficients\n",
    "      4) Plot G-vs-T and residuals\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    csv_file : str or Path\n",
    "        Path to the CSV file (e.g., 'Gibbs_Energies_All.csv').\n",
    "    T_split : float, optional\n",
    "        Temperature (K) at which to split the fit (default 500.0).\n",
    "    output_file : str, optional\n",
    "        Text file to save coefficients (default 'nasa9_coefficients.txt').\n",
    "    \"\"\"\n",
    "    T, species = load_csv(csv_file)\n",
    "\n",
    "    with open(output_file, 'w') as fout:\n",
    "        for name, G in species.items():\n",
    "            print(f\"→ Fitting {name}\")\n",
    "\n",
    "            # Fit low/high ranges\n",
    "            a_low, a_high = fit_nasa9(T, G, T_split=T_split)\n",
    "\n",
    "            # Write coefficients for this species\n",
    "            fout.write(f\"{name}\\n\")\n",
    "            fout.write(\"[\" + \", \".join(f\"{x:.6e}\" for x in a_low)  + \"]\\n\")\n",
    "            fout.write(\"[\" + \", \".join(f\"{x:.6e}\" for x in a_high) + \"]\\n\\n\")\n",
    "\n",
    "            # Compute fitted curves for plotting\n",
    "            mask_low  = T <= T_split\n",
    "            mask_high = T >  T_split\n",
    "            G_low  = compute_gibbs(a_low,  T[mask_low])\n",
    "            G_high = compute_gibbs(a_high, T[mask_high])\n",
    "\n",
    "            # --- Plot: Observed vs Fitted and residuals ---\n",
    "            plt.figure(figsize=(10, 4))\n",
    "\n",
    "            # (1) G vs T\n",
    "            plt.subplot(1, 2, 1)\n",
    "            plt.plot(T,            G,        'ko', label='Observed G(T)')\n",
    "            plt.plot(T[mask_low],  G_low,    'b-', label='Low-T fit')\n",
    "            plt.plot(T[mask_high], G_high,   'r-', label='High-T fit')\n",
    "            plt.axvline(T_split, color='k', ls='--', label=f'T_split={T_split:g} K')\n",
    "            plt.xlabel('T (K)')\n",
    "            plt.ylabel('G (J/mol)')\n",
    "            plt.title(name)\n",
    "            plt.legend()\n",
    "\n",
    "            # (2) Residuals\n",
    "            res_low  = G_low  - G[mask_low]\n",
    "            res_high = G_high - G[mask_high]\n",
    "\n",
    "            plt.subplot(1, 2, 2)\n",
    "            plt.plot(T[mask_low],  res_low,  'bo-', label='Residuals (Low)')\n",
    "            plt.plot(T[mask_high], res_high, 'ro-', label='Residuals (High)')\n",
    "            plt.axhline(0, color='k', ls='--')\n",
    "            plt.xlabel('T (K)')\n",
    "            plt.ylabel('Fit − Observed (J/mol)')\n",
    "            plt.legend()\n",
    "\n",
    "            plt.tight_layout()\n",
    "            plt.show()\n",
    "\n",
    "    print(f\"\\nAll coefficients saved to {output_file}\")\n",
    "\n",
    "\n",
    "# === Parameters ===\n",
    "# The file named 'Gibbs_Energies_All.csv' contains all Gibbs energy data (J/mol)\n",
    "# that this code will ingest to perform the polynomial fitting.\n",
    "csv_file    = 'Gibbs_Energies_All.csv'      # path to your data file\n",
    "T_split     = 500.0                         # split temperature (K) for low/high fits\n",
    "output_file = 'nasa9_coefficients.txt'      # where to write coefficients\n",
    "\n",
    "# Run everything\n",
    "run_all(csv_file, T_split=T_split, output_file=output_file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15cd2d54-c20f-49d9-be81-afbc11d429f0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (P1_ThermochemicalModels)",
   "language": "python",
   "name": "p1_thermochemicalmodels"
  },
  "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.13.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
