{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3b501644",
   "metadata": {},
   "source": [
    "### What this script does \n",
    "\n",
    "- Loads all 10-min merged sonic/mast/radiation datasets and 10-min rain files.\n",
    "\n",
    "- Merges rain onto the sonic time series.\n",
    "\n",
    "- Flags/filters out rainy or inconsistent records (temp/wind sanity checks).\n",
    "\n",
    "- Computes u*, θ*, Monin–Obukhov length L, ζ = z/L, vertical gradients, and stability functions ϕm and ϕh.\n",
    "\n",
    "- Compares measured ϕm, ϕh to Businger–Dyer curves (stable/unstable), reports deviations, and plots nice side-by-side panels.\n",
    "\n",
    "- Saves the final figure to disk.\n",
    "\n",
    "#### Lines you’ll likely need to edit:\n",
    "1) Your data locations:\n",
    "\n",
    "    - DATA_ROOT = r\"C:\\path\\to\\your\\Sonic\"   <-- change if Sonic path differs\n",
    "    - RAIN_ROOT = r\"C:\\path\\to\\your\\Cloud_radar\"  <-- change if Rain path differs\n",
    "\n",
    "2) Which months to include (folders like 2024-03, 2024-04, …)\n",
    "   \n",
    "   DATE_GLOB = '2024-0[3-5]'  <-- e.g., '2024-0[3-6]' to include June\n",
    "\n",
    "3) Site / constants (use your sonic height)\n",
    "- Z_SONIC = 2.99     <-- set to your sonic height (m)\n",
    "- KAPPA, G, theta0 = 0.4, 9.81, 288  <--leave unless you need different constants\n",
    "\n",
    "4) Expected file names inside rain folders: rain_pattern = os.path.join(RAIN_ROOT, DATE_GLOB, '**', 'Rain_10min_Averages.csv')<-- change 'Rain_10min_Averages.csv' if your filename differs\n",
    "\n",
    "5) QC thresholds (tune if needed)\n",
    "\n",
    "6) Optional clear-sky filter \n",
    "\n",
    "7) Output plot location/name: output_path = os.path.join(DATA_ROOT, 'stability_functions.jpg')  <-- change name/path if desired"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c407b11",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, glob\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.polynomial.polynomial import Polynomial\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b298785d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Edit before running!!!\n",
    "# 1) Your data locations\n",
    "# Where your merged_data_10min.csv files live (tower/sonic)\n",
    "# Example: DATA_ROOT = r\"D:\\MyThesis\\data\\Sonic\"\n",
    "DATA_ROOT = r\"C:\\path\\to\\your\\Sonic\"      # <-- change if Sonic path differs\n",
    "# Where your rain files live\n",
    "# Example: RAIN_ROOT = r\"D:\\MyThesis\\data\\Cloud_radar\"\n",
    "RAIN_ROOT = r\"C:\\path\\to\\your\\Cloud_radar\" # <-- change if Rain path differs\n",
    "\n",
    "\n",
    "# Date range of interest (folders named 2024-03, 2024-04, 2024-05, 2024-06)\n",
    "DATE_GLOB = '2024-0[3-5]'\n",
    "\n",
    "# Sonic height, von Kármán constant, gravity\n",
    "Z_SONIC = 2.99\n",
    "KAPPA, G,theta0 = 0.4, 9.81, 288\n",
    "\n",
    "# ── 1) LOAD ALL TOWER/SONIC DATA (10-min merged) ───────────────────────────────\n",
    "\n",
    "sonic_pattern = os.path.join(DATA_ROOT, '**', 'merged_data_10min.csv')\n",
    "sonic_files = glob.glob(sonic_pattern, recursive=True)\n",
    "if not sonic_files:\n",
    "    raise RuntimeError(f\"No sonic files found under {DATA_ROOT}\")\n",
    "\n",
    "df = pd.concat(\n",
    "    (pd.read_csv(fn, parse_dates=['TIMESTAMP']) for fn in sonic_files),\n",
    "    ignore_index=True\n",
    ").sort_values('TIMESTAMP')\n",
    "print(f\"Loaded {len(sonic_files)} sonic files → {len(df)} rows\")\n",
    "\n",
    "# ── 2) LOAD ALL RAIN DATA (10-min averages) ────────────────────────────────────\n",
    "\n",
    "rain_pattern = os.path.join(RAIN_ROOT, DATE_GLOB, '**', 'Rain_10min_Averages.csv')\n",
    "rain_files = glob.glob(rain_pattern, recursive=True)\n",
    "if not rain_files:\n",
    "    raise RuntimeError(f\"No rain files found under {RAIN_ROOT}/{DATE_GLOB}\")\n",
    "\n",
    "rain = pd.concat(\n",
    "    (pd.read_csv(fn, parse_dates=['TIMESTAMP']) for fn in rain_files),\n",
    "    ignore_index=True\n",
    ").sort_values('TIMESTAMP')\n",
    "print(f\"Loaded {len(rain_files)} rain files → {len(rain)} rows\")\n",
    "\n",
    "# ── 3) MERGE RAIN INTO SONIC BY TIMESTAMP ───────────────────────────────────────\n",
    "\n",
    "df = pd.merge_asof(\n",
    "    rain, \n",
    "    df,\n",
    "    on='TIMESTAMP',\n",
    "    direction='nearest'\n",
    ")\n",
    "\n",
    "print(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f22744a",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1824c0e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# ── After merging rain into df via merge_asof ───────────────────────────────────\n",
    "\n",
    "# 3) Flag rainy periods\n",
    "df['Flag_Rain'] = (df['Rain'] > 0).astype(int)\n",
    "\n",
    "# 4) Basic QC filters\n",
    "#    Flag large T mismatches between sonic and corrected mast temp\n",
    "df['Temp_Diff']   = df['Temperature_K_2.99'] - df['Average_Temperature_Corr']\n",
    "df['Flag_Temp']   = (df['Temp_Diff'].abs() > 1).astype(int)\n",
    "# Windspeed mismatch flag: sonic vs. mast\n",
    "df['Windspeed_Diff']  = df['WS_ms_D15014_Avg'] - df['Wind_Speed']\n",
    "df['Flag_Windspeed']  = (df['Windspeed_Diff'].abs() > 1).astype(int)\n",
    "\n",
    "# 2) Drop any record flagged for rain, temp, or wind‐speed\n",
    "df = df[\n",
    "    (df['Flag_Rain']      == 0) &\n",
    "    (df['Flag_Temp']      == 0) &\n",
    "    (df['Flag_Windspeed'] == 0)\n",
    "]\n",
    "\n",
    "#    Keep only dry, non-spurious records\n",
    "#df = df[(df['Flag_Rain'] == 0) & (df['Flag_Temp'] == 0)]\n",
    "\n",
    "# 4) drop NaNs/infs in key fields\n",
    "df = df.replace([np.inf, -np.inf], np.nan)\n",
    "df = df.dropna(subset=[\n",
    "    'uw_flux_corr','vw_flux_corr','wT_Flux',\n",
    "    'Virtual_Dry_Static_Energy_2','Virtual_Dry_Static_Energy_10',\n",
    "    'WS_ms_D15008_Avg','WS_ms_D15463_Avg','Dry_Static_Energy_2.99','Average_Temperature_Corr'\n",
    "])\n",
    "\n",
    "# 5) compute u*\n",
    "#df['u_star'] = ((df['uw_flux']**2 + df['vw_flux']**2)**0.5)**0.5\n",
    "\n",
    "\n",
    "# θ at 2 m and sonic (2.99 m)\n",
    "df['theta_2']      = df['Virtual_Dry_Static_Energy_2']    \n",
    "df['theta_sonic']  = df['Virtual_Dry_Static_Energy_10'] \n",
    "\n",
    "# 5) Compute flux scales and stability parameters\n",
    "\n",
    "# friction velocity\n",
    "df['u_star']     = ((df['uw_flux_corr']**2 + df['vw_flux_corr']**2)**0.5)**0.5\n",
    "\n",
    "# temperature scale\n",
    "df['theta_star'] = - df['wT_Flux'] / df['u_star']\n",
    "\n",
    "# Monin–Obukhov length (using 10 m temp in K)\n",
    "df['L']          = - df['u_star']**3 / (KAPPA * G / df['Average_Temperature_Corr']  * df['wT_Flux'])\n",
    "\n",
    "# stability parameter\n",
    "df['zeta']       = Z_SONIC / df['L']\n",
    "\n",
    "def classify_stability(zeta):\n",
    "    if zeta < 0:\n",
    "        return 'Unstable'\n",
    "    elif zeta > 0:\n",
    "        return 'Stable'\n",
    "    else:\n",
    "        return 'Near-neutral'\n",
    "\n",
    "df['stability_class'] = df['zeta'].apply(classify_stability)\n",
    "print(df['stability_class'].value_counts(normalize=True) * 100)\n",
    "\n",
    "\n",
    "# 6) Vertical gradients by finite differences\n",
    "\n",
    "# ∂U/∂z between 2 m and 10 m\n",
    "df['dU_dz'] = (df['WS_ms_D15463_Avg'] - df['WS_ms_D15008_Avg']) / (10.0 - 2.0)\n",
    "\n",
    "# ∂θ/∂z between 2 m and 10 m (Temperatures already in K)\n",
    "df['dT_dz'] = (df['theta_sonic'] - df['theta_2'])     / (10.0  - 2.0)\n",
    "\n",
    "\n",
    "# 7) Stability functions and final QC\n",
    "\n",
    "df['phi_m'] = KAPPA * Z_SONIC / df['u_star']     * df['dU_dz']\n",
    "df['phi_h'] = KAPPA * Z_SONIC / df['theta_star'] * df['dT_dz']\n",
    "\n",
    "# remove low‐turbulence and extremes\n",
    "df = df[df['u_star'] > 0.1]\n",
    "df = df[(df['zeta'] > -1) & (df['zeta'] < 1)]\n",
    "df = df[(df['phi_m'] > 0) & (df['phi_m'] < 10)]\n",
    "df = df[(df['phi_h'] > 0) & (df['phi_h'] < 10)]\n",
    "\n",
    "\n",
    "# 8) Campaign‐wide stability‐function plot\n",
    "\n",
    "ζ       = np.linspace(-1, 1, 400)\n",
    "φm_unst = (1 - 15*ζ)**(-1/4)\n",
    "φh_unst = 0.74*((1 - 9*ζ)**(-1/2))\n",
    "φm_st   = 1 + 4.7*ζ\n",
    "φh_st   = 0.74 + 4.7*ζ\n",
    "\n",
    "#df = df[(df['phi_m'] > 0) & (df['phi_h'] > 0) ]#& (df['u_star'] > 0.1)]\n",
    "\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5), sharey=True)\n",
    "\n",
    "# φ_m panel\n",
    "ax1.scatter(df['zeta'], df['phi_m'], s=8, alpha=0.3, label='measured')\n",
    "ax1.plot(ζ[ζ<0], φm_unst[ζ<0], 'k--', label='BD unstable')\n",
    "ax1.plot(ζ[ζ>0], φm_st[ζ>0],   'k-',  label='BD stable')\n",
    "ax1.set(xlabel=r'$\\zeta=z/L$', ylabel=r'$\\phi_m$', title='Momentum Stability')\n",
    "ax1.legend(); ax1.grid()\n",
    "\n",
    "# φ_h panel\n",
    "ax2.scatter(df['zeta'], df['phi_h'], s=8, alpha=0.3, label='measured')\n",
    "ax2.plot(ζ[ζ<0], φh_unst[ζ<0], 'k--', label='BD unstable')\n",
    "ax2.plot(ζ[ζ>0], φh_st[ζ>0],   'k-',  label='BD stable')\n",
    "ax2.set(xlabel=r'$\\zeta=z/L$', ylabel=r'$\\phi_h$', title='Heat Stability')\n",
    "ax2.legend(); ax2.grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23033ca2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- Evaluate deviation from Businger–Dyer curve for phi_m in unstable regime (ζ < -0.1) ---\n",
    "df_unstable = df[df['zeta'] < 0].copy()\n",
    "df_unstable['phi_m_BD'] = (1 - 15 * df_unstable['zeta'])**(-1/4)\n",
    "df_unstable['rel_diff'] = (df_unstable['phi_m'] - df_unstable['phi_m_BD']) / df_unstable['phi_m_BD']\n",
    "\n",
    "# Median and mode (via binning)\n",
    "median_rel_diff = df_unstable['rel_diff'].median()\n",
    "mode_bin = df_unstable['rel_diff'].round(2).mode().iloc[0]\n",
    "\n",
    "print(f\"Median deviation: {median_rel_diff:.2%}\")\n",
    "print(f\"Mode deviation: {mode_bin:.2%}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3987ab01",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate deviation from Businger–Dyer curve for phi_m in stable regime (ζ > 0)\n",
    "df_stable = df[df['zeta'] > 0].copy()\n",
    "df_stable['phi_m_BD'] = 1 + 4.7 * df_stable['zeta']\n",
    "df_stable['rel_diff'] = (df_stable['phi_m'] - df_stable['phi_m_BD']) / df_stable['phi_m_BD']\n",
    "\n",
    "# Median and mode (via binning)\n",
    "median_rel_diff_stable = df_stable['rel_diff'].median()\n",
    "mode_bin_stable = df_stable['rel_diff'].round(2).mode().iloc[0]\n",
    "num_points_stable = len(df_stable)\n",
    "\n",
    "median_rel_diff_stable, mode_bin_stable, num_points_stable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b396389f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_stable['phi_h_BD']= 0.74 + 4.7*df_stable['zeta']\n",
    "\n",
    "df_stable['diff'] = (df_stable['phi_h'] - df_stable['phi_h_BD']) / df_stable['phi_h_BD']\n",
    "print(# Median and mode (via binning)\n",
    "df_stable['diff'].median(),\n",
    "df_stable['diff'].round(2).mode().iloc[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6ea540e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import mode\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd82f534",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ensure clear-sky filter is available and remove missing or out-of-bounds values\n",
    "df_cs = df[\n",
    "  #  (df['CSI'] > 0.7) &\n",
    "    (df['zeta'] > -1) & (df['zeta'] < 1) &\n",
    "    (df['phi_h'] > 0) & (df['phi_h'] < 10)\n",
    "].copy()\n",
    "\n",
    "# Define theoretical BD curve for φh\n",
    "df_cs['phi_h_BD'] = np.where(\n",
    "    df_cs['zeta'] < 0,\n",
    "    0.74 * (1 - 9 * df_cs['zeta'])**(-0.5),\n",
    "    0.74 + 4.7 * df_cs['zeta']\n",
    ")\n",
    "\n",
    "\n",
    "# Compute percent deviation\n",
    "df_cs['phi_h_dev'] = 100 * (df_cs['phi_h'] - df_cs['phi_h_BD']) / df_cs['phi_h_BD']\n",
    "\n",
    "# Median and mode\n",
    "median_dev = df_cs['phi_h_dev'].median()\n",
    "mode_dev = mode(df_cs['phi_h_dev'].round()).mode  # Round to nearest integer for mode\n",
    "\n",
    "print(f\"Median deviation: {median_dev:.2f}%\")\n",
    "print(f\"Mode deviation: {mode_dev:.2f}%\")\n",
    "print(f\"N points: {len(df_cs)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ac65302",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define clear-sky daytime: CSI > 0.7 and SW_down > 50 W/m²\n",
    "df_clear_sky = df[(df['CSI'] > 0.7)]# & (df['SR15D1Dn_Irr'] > 50)]  # Replace 'SW↓' with actual column name\n",
    "df_cs_unstable = df_clear_sky[df_clear_sky['zeta'] < -0.1].copy()\n",
    "df_cs_unstable['phi_m_BD'] = (1 - 15 * df_cs_unstable['zeta'])**(-1/4)\n",
    "df_cs_unstable['rel_diff'] = (df_cs_unstable['phi_m'] - df_cs_unstable['phi_m_BD']) / df_cs_unstable['phi_m_BD']\n",
    "median_cs = df_cs_unstable['rel_diff'].median()\n",
    "mode_cs = df_cs_unstable['rel_diff'].round(2).mode().iloc[0]\n",
    "within_15_cs = (df_cs_unstable['rel_diff'].abs() <= 0.15).mean() * 100\n",
    "print(f\"Median deviation: {median_cs:.2%}\")\n",
    "print(f\"Mode deviation: {mode_cs:.2%}\")\n",
    "print(within_15_cs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a063368",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Create figure and subplots\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 7), sharey=True)\n",
    "\n",
    "# Common plot parameters\n",
    "scatter_kwargs = dict(s=20, alpha=0.5, edgecolors='none')\n",
    "line_kwargs_stable = dict(color='black', linewidth=2.5, label='BD stable')\n",
    "line_kwargs_unstable = dict(color='black', linestyle='--', linewidth=2.5, label='BD unstable')\n",
    "\n",
    "# --- φ_m subplot ---\n",
    "ax1.scatter(df['zeta'], df['phi_m'], **scatter_kwargs, label='Measured')\n",
    "ax1.plot(ζ[ζ < 0], φm_unst[ζ < 0], **line_kwargs_unstable)\n",
    "ax1.plot(ζ[ζ > 0], φm_st[ζ > 0], **line_kwargs_stable)\n",
    "ax1.set_xlabel(r'$\\zeta = z/L$', fontsize=20)\n",
    "ax1.set_ylabel(r'$\\phi_m$', fontsize=20)\n",
    "ax1.set_title('Momentum Stability Function', fontsize=20)\n",
    "ax1.tick_params(labelsize=18)\n",
    "ax1.legend(loc='upper left', fontsize=18)\n",
    "ax1.grid(True)\n",
    "\n",
    "# --- φ_h subplot ---\n",
    "ax2.scatter(df['zeta'], df['phi_h'], **scatter_kwargs, label='Measured')\n",
    "ax2.plot(ζ[ζ < 0], φh_unst[ζ < 0], **line_kwargs_unstable)\n",
    "ax2.plot(ζ[ζ > 0], φh_st[ζ > 0], **line_kwargs_stable)\n",
    "ax2.set_xlabel(r'$\\zeta = z/L$', fontsize=20)\n",
    "ax2.set_ylabel(r'$\\phi_h$', fontsize=20)\n",
    "ax2.set_title('Heat Stability Function', fontsize=20)\n",
    "ax2.tick_params(labelsize=18)\n",
    "ax2.legend(loc='upper left', fontsize=18)\n",
    "ax2.grid(True)\n",
    "\n",
    "# Final layout\n",
    "#fig.suptitle('Stability Functions Compared to Businger–Dyer Curves', fontsize=22)\n",
    "plt.tight_layout(rect=[0, 0, 1, 0.95])\n",
    "output_path = os.path.join(DATA_ROOT, 'stability_functions.jpg')\n",
    "plt.savefig(output_path, dpi=300)  # save at high resolution for report\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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}