{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bca6f2bd",
   "metadata": {},
   "source": [
    "\n",
    "# Cantera Equilibrium Skeleton (TP, VCS Minimization)\n",
    "\n",
    "This notebook is a minimal, self-contained scaffold showing how we set up and run\n",
    "thermodynamic equilibrium calculations in Cantera for aqueous systems (e.g., Titan\n",
    "impact-pond analogs). It is intended for readers and reviewers who want to understand\n",
    "the exact simulation flow and reproduce results with their own species YAML files.\n",
    "\n",
    "**What this does**\n",
    "- Loads an aqueous phase from a provided YAML (with NASA-9 polynomials).\n",
    "- Defines pressure, a temperature list (in °C), and an initial composition.\n",
    "- Normalizes the composition, sets the `TPX` state, and equilibrates using the VCS solver.\n",
    "- Records final mole fractions for all species vs. temperature to a CSV file.\n",
    "\n",
    "**What you must provide**\n",
    "- A Cantera YAML file containing the relevant species and an `aqueous` phase definition.\n",
    "- The folder path containing your YAML(s). If your YAMLs depend on relative includes,\n",
    "  place them in the same directory and call `ct.add_directory(...)` accordingly.\n",
    "\n",
    "**Assumptions**\n",
    "- Thermodynamic data are represented as NASA-9 polynomials fitted from a database\n",
    "  (e.g., CHNOSZ-derived fits). The equilibrium is obtained by **Gibbs energy minimization**\n",
    "  using Cantera's `VCS` algorithm under constant temperature and pressure (`\"TP\"`).\n",
    "\n",
    "> **Tip:** This is a skeleton. Start here, then expand (e.g., NH3 sweeps, multiple YAMLs,\n",
    "> alternate pressures, additional outputs).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cd9f87f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# =============================================================\n",
    "# Minimal Cantera Equilibrium Skeleton\n",
    "# =============================================================\n",
    "# Purpose: demonstrate the essential steps to (1) load a phase,\n",
    "# (2) define state & composition, (3) equilibrate with VCS at TP,\n",
    "# (4) export species mole fractions vs. temperature.\n",
    "#\n",
    "# Expand as needed for your study (parameter sweeps, plots, etc.).\n",
    "# -------------------------------------------------------------\n",
    "\n",
    "from pathlib import Path\n",
    "import cantera as ct\n",
    "import pandas as pd\n",
    "\n",
    "# ---------------- User Configuration (EDIT THESE) ----------------\n",
    "# Path that contains your YAML files (and any includes they need).\n",
    "ROOT_DIR   = Path.cwd()  # e.g., Path('Input files/with_NH3')\n",
    "PHASE_FILE = 'filename.yaml'   # <-- replace with your YAML filename\n",
    "PHASE_NAME = 'aqueous'              # <-- replace if your phase name differs\n",
    "\n",
    "PRESSURE_PA = 1.5e5                 # Example: 1.5 bar in Pascals\n",
    "TEMPS_C     = [0, 25, 50, 75, 100]  # Temperatures to evaluate (°C)\n",
    "\n",
    "# Initial amounts. Use any convenient totals;\n",
    "# they will be normalized to mole fractions automatically.\n",
    "BASE_COMP = {\n",
    "    'H2O(l)'  : 1.00,\n",
    "    'HCN(aq)' : 0.02,\n",
    "    'C2H2(aq)': 0.04,\n",
    "    # 'NH3(aq)': 0.05,  # Uncomment/edit to include ammonia, etc.\n",
    "    # Add/remove species keys to match your YAML species names.\n",
    "}\n",
    "\n",
    "OUTPUT_CSV = 'equilibrium_results_skeleton.csv'  # saved in notebook folder\n",
    "# ----------------------------------------------------------------\n",
    "\n",
    "# Make the YAML directory discoverable (important for relative includes)\n",
    "ct.add_directory(str(ROOT_DIR))\n",
    "\n",
    "# Load the aqueous phase\n",
    "sol = ct.Solution(PHASE_FILE, PHASE_NAME)\n",
    "\n",
    "# Helper: convert unnormalized amounts to mole fractions\n",
    "def normalize_amounts(amounts: dict) -> dict:\n",
    "    total = sum(amounts.values())\n",
    "    if total <= 0:\n",
    "        raise ValueError(\"Total initial amount must be positive.\")\n",
    "    return {k: v/total for k, v in amounts.items()}\n",
    "\n",
    "# Prepare output table columns: Temperature + all species in phase\n",
    "columns = ['Temperature_C'] + list(sol.species_names)\n",
    "rows = []\n",
    "\n",
    "# Loop over temperatures and run TP equilibrium (VCS)\n",
    "for T_C in TEMPS_C:\n",
    "    T_K = T_C + 273.15\n",
    "    # Fresh solution object each iteration (keeps each state independent)\n",
    "    aq = ct.Solution(PHASE_FILE, PHASE_NAME)\n",
    "\n",
    "    # Normalize composition and set state\n",
    "    comp_X = normalize_amounts(BASE_COMP)\n",
    "    aq.TPX = T_K, PRESSURE_PA, comp_X\n",
    "\n",
    "    # Equilibrate using VCS (robust for multi-component systems)\n",
    "    aq.equilibrate('TP', solver='vcs', max_steps=100000)\n",
    "\n",
    "    # Collect mole fractions for every species at this temperature\n",
    "    row = [T_C] + [aq.X[aq.species_index(sp)] for sp in aq.species_names]\n",
    "    rows.append(row)\n",
    "\n",
    "# Save to CSV for provenance/reuse\n",
    "df = pd.DataFrame(rows, columns=columns)\n",
    "df.to_csv(OUTPUT_CSV, index=False)\n",
    "\n",
    "print(f\"  Done. Wrote {len(df)} rows × {len(df.columns)} columns to {OUTPUT_CSV}\")\n",
    "print(\"   Phase file:\", PHASE_FILE, \"| Phase name:\", PHASE_NAME)\n",
    "print(\"   Temperatures (°C):\", TEMPS_C)\n",
    "\n",
    "# ---------------- Optional: quick look (uncomment to preview) ----------------\n",
    "# import matplotlib.pyplot as plt\n",
    "# plt.figure(figsize=(9,5))\n",
    "# for sp in sol.species_names:\n",
    "#     plt.plot(df['Temperature_C'], df[sp], label=sp)\n",
    "# plt.yscale('log')\n",
    "# plt.xlabel('Temperature (°C)')\n",
    "# plt.ylabel('Equilibrium mole fraction (log scale)')\n",
    "# plt.title('Equilibrium composition vs. Temperature (skeleton run)')\n",
    "# plt.legend(bbox_to_anchor=(1.02, 1), loc='upper left')\n",
    "# plt.tight_layout()\n",
    "# plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
