Na₂SO₄ Phase Diagram Tutorial
Steiger & Asmussen (2008) theory × IoT T–RH monitoring of Huang Tu-shui's Water Buffaloes

"""
steiger_na2so4.py
=================
Na2SO4-H2O 相圖 (Steiger & Asmussen 2008, Geochim. Cosmochim. Acta 72, 4291-4306)
Python 移植版 — 直接翻譯自論文附件 mmc2.xls 中的 VBA 程式 (Modul1.bas)

模型: Pitzer 離子交互作用模型
  - aphi(T)      : Debye-Hueckel 滲透係數參數
  - beta0/1/2, cphi : Na2SO4 的 Pitzer 參數 (溫度函數)
  - phi(T,m)     : 滲透係數
  - aw(T,m)      : 水活性  =>  平衡相對濕度 RH_eq = 100 * aw
  - lnK_V / lnK_III / lnK_10 / lnK_7 : 各固相的溶解平衡常數
      V   = thenardite (無水芒硝, Na2SO4(V))
      III = Na2SO4(III) (介穩)
      10  = mirabilite (芒硝, Na2SO4·10H2O)
      7   = heptahydrate (七水合物, 介穩)

作者: 依 Tai-Sheng 需求由 Claude 產生, 2026-07-07
"""

import numpy as np
from scipy.optimize import brentq

# ------------------------------------------------------------------
# 常數 (與 VBA 完全相同)
# ------------------------------------------------------------------
MW   = 18.01528      # 水分子量 g/mol
NUM  = 2             # Na+ 個數
NUX  = 1             # SO4^2- 個數
NU   = 3             # 總離子數
ABSZ = 2             # |z+ * z-|
IFAKT = 3            # 離子強度因子 I = 3m
ALPHA1 = 1.2
ALPHA2 = 0.2
TR = 298.15          # 參考溫度 K

# ------------------------------------------------------------------
# Pitzer 參數 (溫度函數), Q(1..6) 依 VBA
# ------------------------------------------------------------------
def _Qfunc(T, Q):
    return (Q[0] + Q[1]*(1/T - 1/TR) + Q[2]*np.log(T/TR)
            + Q[3]*(T - TR) + Q[4]*(T**2 - TR**2)
            + (Q[5]*np.log(T - 225) if Q[5] != 0 else 0.0))

def aphi(T):
    """Debye-Hueckel 參數 A_phi"""
    return (-0.8176532151 - 0.8685275541/(T - 222) + 19250.991343/T**2
            + 0.005251284381*T - 7.149396961e-6*T**2 + 9.338558897e-12*T**4)

def beta0(T):
    return _Qfunc(T, [0.500271751403524, 7362.93162373699,
                      46.9944253174359, -7.94134431636842e-2, 0, 0])

def beta1(T):
    return _Qfunc(T, [1.02168434322046, -4221.18387050915,
                      -13.2556426105993, 0, 0, 0])

def beta2(T):
    return _Qfunc(T, [-2.11840490018125, -10429.7675818801, -69.1082762017204,
                      0.118455594720063, 0, 0.314460235765225])

def cphi(T):
    return _Qfunc(T, [-1.43154740436143e-2, 0, -0.125170094651518,
                      5.22666755731633e-4, 0, 0])

# ------------------------------------------------------------------
# 滲透係數與水活性
# ------------------------------------------------------------------
def _I(m):
    return IFAKT * m

def phi(T, m):
    """滲透係數 (osmotic coefficient)"""
    x1 = ALPHA1 * np.sqrt(_I(m))
    x2 = ALPHA2 * np.sqrt(_I(m))
    fphi = -aphi(T) * np.sqrt(_I(m)) / (1 + 1.2*np.sqrt(_I(m)))
    Bphi = beta0(T) + beta1(T)*np.exp(-x1) + beta2(T)*np.exp(-x2)
    return (1 + ABSZ*fphi + 2*NUM*NUX/NU * m * Bphi
            + 2*(NUM*NUX)**1.5/NU * m**2 * cphi(T))

def aw(T, m):
    """水活性 a_w;  平衡 RH(%) = 100 * a_w"""
    return np.exp(-phi(T, m) * MW/1000 * NU * m)

def lngamma(T, m):
    """平均離子活性係數 ln(gamma±)"""
    sqI = np.sqrt(_I(m))
    x1 = ALPHA1 * sqI
    x2 = ALPHA2 * sqI
    fg = -aphi(T) * (sqI/(1 + 1.2*sqI) + 2/1.2*np.log(1 + 1.2*sqI))
    Bg = (2*beta0(T)
          + 2*beta1(T)/x1**2 * (1 - (1 + x1 - 0.5*x1**2)*np.exp(-x1))
          + 2*beta2(T)/x2**2 * (1 - (1 + x2 - 0.5*x2**2)*np.exp(-x2)))
    return (ABSZ*fg + 2*NUM*NUX/NU * m * Bg
            + 1.5*2*(NUM*NUX)**1.5/NU * m**2 * cphi(T))

# ------------------------------------------------------------------
# 溶解平衡常數 ln K (VBA: lnK_V, lnK_III, lnK_10, lnK_7)
# ------------------------------------------------------------------
def _lnK(T, p):
    return p[0] + p[1]*(1/T - 1/TR) + p[2]*np.log(T/TR) + p[3]*(T - TR)

def lnK_V(T):     # thenardite
    return _lnK(T, [-0.707283559813531, -9860.32992397018, -33.9727082671554, 0])

def lnK_III(T):   # Na2SO4(III), 介穩
    return _lnK(T, [-0.140054987730789, -9565.2275844851, -34.1939432311946, 0])

def lnK_10(T):    # mirabilite
    return _lnK(T, [-2.84554475448134, 0, 79.7760744619095, -0.16100201693158])

def lnK_7(T):     # heptahydrate, 介穩
    return _lnK(T, [-1.61986647528533, -20733.6510183963, -48.9522629908378, 0])

# ------------------------------------------------------------------
# 相邊界計算
# ------------------------------------------------------------------
def sat_molality(T, lnK_func, n_hydrate):
    """求飽和溶解度 m_sat:
       ln K = ln4 + 3 ln(m) + 3 ln(gamma±) + n ln(aw)"""
    f = lambda m: (np.log(4) + 3*np.log(m) + 3*lngamma(T, m)
                   + n_hydrate*np.log(aw(T, m)) - lnK_func(T))
    # 高溫時 f(m) 可能有多個根 (介穩外插區), 取最小的根;
    # 並避開 Pitzer 模型外插至非物理區 (aw > 1)
    grid = np.concatenate([np.array([1e-4]), np.arange(0.05, 12.01, 0.05)])
    fv = np.array([f(m) for m in grid])
    sign_change = np.where(np.sign(fv[:-1]) != np.sign(fv[1:]))[0]
    if len(sign_change) == 0:
        raise ValueError(f'no saturation root at T={T}')
    i = sign_change[0]
    return brentq(f, grid[i], grid[i+1], xtol=1e-10)

def rh_deliq(T, lnK_func, n_hydrate):
    """飽和溶液的潮解相對濕度 RH_eq(%)"""
    m = sat_molality(T, lnK_func, n_hydrate)
    return 100.0 * aw(T, m)

def rh_mirabilite(TC):
    """mirabilite (Na2SO4·10H2O) 潮解線, TC = 攝氏溫度"""
    return rh_deliq(TC + 273.15, lnK_10, 10)

def rh_thenardite(TC):
    """thenardite (Na2SO4, V) 潮解線"""
    return rh_deliq(TC + 273.15, lnK_V, 0)

def rh_heptahydrate(TC):
    """heptahydrate (Na2SO4·7H2O) 潮解線 (介穩)"""
    return rh_deliq(TC + 273.15, lnK_7, 7)

def rh_then_mir(TC):
    """thenardite <-> mirabilite 水合/脫水相邊界:
       Na2SO4·10H2O = Na2SO4(V) + 10 H2O(g)
       ln K10 - ln KV = 10 ln aw  =>  RH = 100*exp((lnK10-lnKV)/10)"""
    T = TC + 273.15
    return 100.0 * np.exp((lnK_10(T) - lnK_V(T)) / 10.0)

def triple_point():
    """mirabilite-thenardite-溶液 三相點 (兩條藍線交點)"""
    f = lambda TC: rh_mirabilite(TC) - rh_then_mir(TC)
    TC = brentq(f, 25, 34)
    return TC, rh_then_mir(TC)

# ------------------------------------------------------------------
# 資料點分區與跨線事件
# ------------------------------------------------------------------
def classify(TC, RH):
    """依相圖將 (T, RH) 分區:
       'solution'   : RH 在潮解線之上   -> 鹽全溶解 (再結晶風險來源)
       'mirabilite' : 兩線之間          -> 芒硝穩定
       'thenardite' : th-mi 線之下       -> 無水芒硝穩定 (脫水)"""
    TC = np.asarray(TC, dtype=float)
    RH = np.asarray(RH, dtype=float)
    Ttp, _ = triple_point()
    upper = np.where(TC <= Ttp,
                     np.interp(TC, _grid_T, _grid_mir),
                     np.interp(TC, _grid_T, _grid_the))
    lower = np.where(TC <= Ttp,
                     np.interp(TC, _grid_T, _grid_tm),
                     np.interp(TC, _grid_T, _grid_the))
    out = np.full(TC.shape, 'mirabilite', dtype=object)
    out[RH >= upper] = 'solution'
    out[RH < lower] = 'thenardite'
    return out

# 預先建表 (0.05 度解析度, 內插用)
_grid_T = np.arange(-0.0, 50.001, 0.05)
_grid_mir = np.array([rh_mirabilite(t) if t <= 35 else np.nan for t in _grid_T])
_grid_the = np.array([rh_thenardite(t) for t in _grid_T])
_grid_tm  = np.array([rh_then_mir(t) for t in _grid_T])
# 三相點以上 mirabilite 線無意義, 以 thenardite 線取代
_Ttp = brentq(lambda tc: rh_mirabilite(tc) - rh_then_mir(tc), 25, 34)
_grid_mir[_grid_T > _Ttp] = _grid_the[_grid_T > _Ttp]

def crossing_events(times, TC, RH):
    """回傳跨越相邊界的事件表 (DataFrame)"""
    import pandas as pd
    reg = classify(TC, RH)
    idx = np.where(reg[1:] != reg[:-1])[0]
    return pd.DataFrame({
        'time_before': np.asarray(times)[idx],
        'time_after':  np.asarray(times)[idx + 1],
        'from_phase':  reg[idx],
        'to_phase':    reg[idx + 1],
        'T_after':     np.asarray(TC)[idx + 1],
        'RH_after':    np.asarray(RH)[idx + 1],
    })

# ------------------------------------------------------------------
# 驗證
# ------------------------------------------------------------------
if __name__ == '__main__':
    print('=== 模型驗證 ===')
    print(f'25C mirabilite 飽和溶解度 m = {sat_molality(298.15, lnK_10, 10):.4f} mol/kg')
    print(f'25C thenardite 飽和溶解度 m = {sat_molality(298.15, lnK_V, 0):.4f} mol/kg')
    print(f'25C mirabilite 潮解 RH = {rh_mirabilite(25):.2f} %')
    print(f'25C thenardite-mirabilite 相邊界 RH = {rh_then_mir(25):.2f} %')
    print(f'20C mirabilite 潮解 RH = {rh_mirabilite(20):.2f} %')
    Ttp, RHtp = triple_point()
    print(f"三相點: T = {Ttp:.2f} C, RH = {RHtp:.2f} %")

# -*- coding: utf-8 -*-
"""
steiger_analysis.py — 黃土水《南國(水牛群像)》IoT 溫濕度資料 vs Na2SO4 相圖
找出資料點跨越兩條藍線 (mirabilite 潮解線 / thenardite-mirabilite 相邊界) 的日期時間
"""
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib import font_manager
for fp in ['/usr/share/fonts/opentype/noto/NotoSansCJK-Regular.ttc',
           '/usr/share/fonts/opentype/noto/NotoSerifCJK-Regular.ttc']:
    import os
    if os.path.exists(fp):
        font_manager.fontManager.addfont(fp)
        plt.rcParams['font.family'] = font_manager.FontProperties(fname=fp).get_name()
        break
plt.rcParams['axes.unicode_minus'] = False
from steiger_na2so4 import (rh_mirabilite, rh_thenardite, rh_then_mir,
                            rh_heptahydrate, triple_point, classify,
                            crossing_events)

# ---------- 讀取資料 (已由 xlsx 轉出) ----------
df = pd.read_csv('iot_raw.csv')

def parse_tw_time(s):
    s = str(s).strip()
    if '上午' in s:
        d, t = s.split('上午')
        return pd.to_datetime(d.strip() + ' ' + t.strip() + ' AM',
                              format='%Y/%m/%d %I:%M:%S %p', errors='coerce')
    if '下午' in s:
        d, t = s.split('下午')
        return pd.to_datetime(d.strip() + ' ' + t.strip() + ' PM',
                              format='%Y/%m/%d %I:%M:%S %p', errors='coerce')
    return pd.to_datetime(s, errors='coerce')

df['datetime'] = df['DataRecordTime'].map(parse_tw_time)
df = df.dropna(subset=['datetime']).sort_values('datetime').reset_index(drop=True)
print('rows:', len(df), 'time range:', df['datetime'].min(), '->', df['datetime'].max())

sensors = {'micro':      ('Micro_T', 'Micro_H', 'Microenvironment (文物微環境)'),
           'peripheral': ('Peri_T',  'Peri_H',  'Peripheral (周邊環境)')}

Ttp, RHtp = triple_point()
Tgrid = np.linspace(0, Ttp, 200)
Tgrid2 = np.linspace(Ttp, 50, 100)
mir_line = np.array([rh_mirabilite(t) for t in Tgrid])
tm_line  = np.array([rh_then_mir(t)  for t in Tgrid])
the_line = np.array([rh_thenardite(t) for t in Tgrid2])

# 相邊界 CSV (供 R / 其他程式使用)
Tall = np.arange(0, 50.5, 0.5)
bound = pd.DataFrame({
    'T': Tall,
    'RH_mirabilite_deliq': [rh_mirabilite(t) if t <= Ttp else np.nan for t in Tall],
    'RH_then_mir':         [rh_then_mir(t)  if t <= Ttp else np.nan for t in Tall],
    'RH_thenardite_deliq': [rh_thenardite(t) for t in Tall],
    'RH_heptahydrate':     [rh_heptahydrate(t) if t <= 23 else np.nan for t in Tall],
})
bound.round(3).to_csv('Na2SO4_RH_vs_T_steiger.csv', index=False)

summary = {}
for key, (tc, hc, label) in sensors.items():
    d = df[['datetime', tc, hc]].dropna().copy()
    d.columns = ['datetime', 'T', 'RH']
    d = d[(d['T'] > -5) & (d['T'] < 55) & (d['RH'] > 0) & (d['RH'] <= 100)]
    d['phase'] = classify(d['T'].values, d['RH'].values)
    ev = crossing_events(d['datetime'].values, d['T'].values, d['RH'].values)
    ev.to_csv(f'crossing_events_{key}.csv', index=False)
    cnt = d['phase'].value_counts()
    # 每一次「進入 solution 區」= 一次潮解事件; 「離開 solution 區」= 一次再結晶事件
    diss  = ev[ev['to_phase'] == 'solution']
    recry = ev[ev['from_phase'] == 'solution']
    dehyd = ev[(ev['from_phase'] == 'mirabilite') & (ev['to_phase'] == 'thenardite')]
    hyd   = ev[(ev['from_phase'] == 'thenardite') & (ev['to_phase'] == 'mirabilite')]
    summary[key] = dict(label=label, n=len(d), counts=cnt.to_dict(),
                        n_cross=len(ev), diss=diss, recry=recry,
                        dehyd=dehyd, hyd=hyd, data=d, events=ev)
    print(f'--- {label} ---')
    print('  區域分布:', cnt.to_dict())
    print(f'  跨線事件總數: {len(ev)}  (潮解 {len(diss)}, 再結晶 {len(recry)}, 脫水 {len(dehyd)}, 水合 {len(hyd)})')
    # 每月跨線統計
    ev2 = ev.copy()
    ev2['month'] = pd.to_datetime(ev2['time_after']).dt.to_period('M')
    print(ev2.groupby('month').size().to_string())

# ---------- 靜態相圖 ----------
for key, s in summary.items():
    d = s['data']
    fig, ax = plt.subplots(figsize=(9, 6), dpi=110)
    colors = {'solution': '#d62728', 'mirabilite': '#2ca02c', 'thenardite': '#ff7f0e'}
    for ph, g in d.groupby('phase'):
        ax.scatter(g['T'], g['RH'], s=4, alpha=.35, color=colors[ph],
                   label=f"{ph} (n={len(g)})")
    ax.plot(Tgrid, mir_line, 'b-', lw=2)
    ax.plot(Tgrid, tm_line, 'b-', lw=2)
    ax.plot(Tgrid2, the_line, 'b-', lw=2)
    ax.axhline(75, color='orange', lw=2)
    ax.text(35, 96, 'Na$_2$SO$_4$ solution'); ax.text(10, 82, 'Mirabilite')
    ax.text(8, 50, 'Thenardite'); ax.text(42, 77, 'NaCl')
    ax.set_xlim(0, 50); ax.set_ylim(15, 100)
    ax.set_xlabel('Temperature (°C)'); ax.set_ylabel('Relative Humidity (%)')
    ax.set_title(f"2021-2022 Huang Water Buffalo — {s['label']}")
    ax.legend(loc='lower right', fontsize=8)
    fig.tight_layout(); fig.savefig(f'phase_diagram_{key}.png'); plt.close(fig)

# ---------- 時間序列圖: 何時在哪個相區 ----------
for key, s in summary.items():
    d = s['data']
    fig, ax = plt.subplots(figsize=(12, 5), dpi=110)
    colors = {'solution': '#d62728', 'mirabilite': '#2ca02c', 'thenardite': '#ff7f0e'}
    for ph, g in d.groupby('phase'):
        ax.scatter(g['datetime'], g['RH'], s=2, color=colors[ph], label=ph)
    ax.set_ylabel('Relative Humidity (%)')
    ax.set_title(f"相區隨時間變化 — {s['label']} (紅=溶解, 綠=芒硝, 橘=無水芒硝)")
    ax.legend(markerscale=6, fontsize=8)
    fig.autofmt_xdate(); fig.tight_layout()
    fig.savefig(f'timeline_{key}.png'); plt.close(fig)

# ---------- 互動式 Plotly (滑鼠移到點上即顯示日期時間) ----------
import plotly.graph_objects as go
for key, s in summary.items():
    d = s['data']
    step = max(1, len(d)//20000)  # 抽樣避免檔案過大
    dd = d.iloc[::step]
    fig = go.Figure()
    colors = {'solution': 'red', 'mirabilite': 'green', 'thenardite': 'orange'}
    for ph, g in dd.groupby('phase'):
        fig.add_trace(go.Scattergl(
            x=g['T'], y=g['RH'], mode='markers', name=f'{ph} (n={len(g)})',
            marker=dict(size=3, color=colors[ph], opacity=.5),
            text=g['datetime'].dt.strftime('%Y-%m-%d %H:%M'),
            hovertemplate='%{text}<br>T=%{x:.2f}°C RH=%{y:.1f}%<extra>'+ph+'</extra>'))
    fig.add_trace(go.Scatter(x=Tgrid, y=mir_line, mode='lines', name='Mirabilite 潮解線',
                             line=dict(color='blue', width=3)))
    fig.add_trace(go.Scatter(x=Tgrid, y=tm_line, mode='lines', name='Thenardite–Mirabilite 相邊界',
                             line=dict(color='blue', width=3, dash='dash')))
    fig.add_trace(go.Scatter(x=Tgrid2, y=the_line, mode='lines', name='Thenardite 潮解線',
                             line=dict(color='blue', width=3)))
    fig.add_hline(y=75, line_color='orange', annotation_text='NaCl 潮解線 75%')
    fig.update_layout(title=f"2021–2022 黃土水《南國(水牛群像)》 — {s['label']}<br><sub>Steiger & Asmussen (2008) Na₂SO₄ 相圖 — 滑鼠移到資料點顯示日期時間</sub>",
                      xaxis_title='Temperature (°C)', yaxis_title='Relative Humidity (%)',
                      xaxis_range=[0,50], yaxis_range=[15,100], template='plotly_white',
                      legend=dict(orientation='h', y=-0.15))
    fig.write_html(f'interactive_phase_{key}.html', include_plotlyjs='cdn')

# 儲存摘要供 HTML 使用
import json
out = {}
for key, s in summary.items():
    ev = s['events'].copy()
    ev['time_after'] = pd.to_datetime(ev['time_after']).dt.strftime('%Y-%m-%d %H:%M')
    ev['time_before'] = pd.to_datetime(ev['time_before']).dt.strftime('%Y-%m-%d %H:%M')
    out[key] = dict(label=s['label'], n=int(s['n']),
                    counts={k: int(v) for k, v in s['counts'].items()},
                    n_cross=int(s['n_cross']),
                    events=ev.round(2).to_dict('records'))
json.dump(out, open('summary.json', 'w'), ensure_ascii=False, indent=1)
print('DONE')

# =====================================================================
# steiger_na2so4.R
# Na2SO4-H2O 相圖 (Steiger & Asmussen 2008, GCA 72, 4291-4306)
# R 移植版 — 直接翻譯自論文附件 mmc2.xls 的 VBA 程式 (Modul1.bas)
# 並分析黃土水《南國(水牛群像)》IoT 溫濕度資料跨越相邊界的日期時間
# =====================================================================

# ---------------- Pitzer 模型常數 (與 VBA 相同) ----------------
MW <- 18.01528; NUM <- 2; NUX <- 1; NU <- 3
ABSZ <- 2; IFAKT <- 3; ALPHA1 <- 1.2; ALPHA2 <- 0.2; TR <- 298.15

Qfunc <- function(T, Q) {
  v <- Q[1] + Q[2]*(1/T - 1/TR) + Q[3]*log(T/TR) +
       Q[4]*(T - TR) + Q[5]*(T^2 - TR^2)
  if (Q[6] != 0) v <- v + Q[6]*log(T - 225)
  v
}

aphi <- function(T)
  -0.8176532151 - 0.8685275541/(T-222) + 19250.991343/T^2 +
   0.005251284381*T - 7.149396961e-6*T^2 + 9.338558897e-12*T^4

beta0 <- function(T) Qfunc(T, c(0.500271751403524, 7362.93162373699,
                                46.9944253174359, -7.94134431636842e-2, 0, 0))
beta1 <- function(T) Qfunc(T, c(1.02168434322046, -4221.18387050915,
                                -13.2556426105993, 0, 0, 0))
beta2 <- function(T) Qfunc(T, c(-2.11840490018125, -10429.7675818801,
                                -69.1082762017204, 0.118455594720063, 0,
                                0.314460235765225))
cphi  <- function(T) Qfunc(T, c(-1.43154740436143e-2, 0, -0.125170094651518,
                                5.22666755731633e-4, 0, 0))

# ---------------- 滲透係數 / 水活性 / 活性係數 ----------------
phi_osm <- function(T, m) {
  I  <- IFAKT*m; sqI <- sqrt(I)
  x1 <- ALPHA1*sqI; x2 <- ALPHA2*sqI
  fphi <- -aphi(T)*sqI/(1 + 1.2*sqI)
  Bphi <- beta0(T) + beta1(T)*exp(-x1) + beta2(T)*exp(-x2)
  1 + ABSZ*fphi + 2*NUM*NUX/NU*m*Bphi + 2*(NUM*NUX)^1.5/NU*m^2*cphi(T)
}

aw <- function(T, m) exp(-phi_osm(T, m)*MW/1000*NU*m)   # RH_eq(%) = 100*aw

lngamma <- function(T, m) {
  I <- IFAKT*m; sqI <- sqrt(I)
  x1 <- ALPHA1*sqI; x2 <- ALPHA2*sqI
  fg <- -aphi(T)*(sqI/(1 + 1.2*sqI) + 2/1.2*log(1 + 1.2*sqI))
  Bg <- 2*beta0(T) +
        2*beta1(T)/x1^2*(1 - (1 + x1 - 0.5*x1^2)*exp(-x1)) +
        2*beta2(T)/x2^2*(1 - (1 + x2 - 0.5*x2^2)*exp(-x2))
  ABSZ*fg + 2*NUM*NUX/NU*m*Bg + 1.5*2*(NUM*NUX)^1.5/NU*m^2*cphi(T)
}

# ---------------- 溶解平衡常數 ln K ----------------
lnKf <- function(T, p) p[1] + p[2]*(1/T - 1/TR) + p[3]*log(T/TR) + p[4]*(T - TR)
lnK_V   <- function(T) lnKf(T, c(-0.707283559813531, -9860.32992397018, -33.9727082671554, 0))  # thenardite
lnK_III <- function(T) lnKf(T, c(-0.140054987730789, -9565.2275844851, -34.1939432311946, 0))
lnK_10  <- function(T) lnKf(T, c(-2.84554475448134, 0, 79.7760744619095, -0.16100201693158))    # mirabilite
lnK_7   <- function(T) lnKf(T, c(-1.61986647528533, -20733.6510183963, -48.9522629908378, 0))   # heptahydrate

# ---------------- 相邊界 ----------------
sat_molality <- function(T, lnK, n_hyd) {
  f <- function(m) log(4) + 3*log(m) + 3*lngamma(T, m) + n_hyd*log(aw(T, m)) - lnK(T)
  grid <- c(1e-4, seq(0.05, 12, by = 0.05))
  fv <- sapply(grid, f)
  i <- which(sign(fv[-1]) != sign(fv[-length(fv)]))[1]
  if (is.na(i)) return(NA_real_)
  uniroot(f, c(grid[i], grid[i+1]), tol = 1e-10)$root
}

rh_deliq        <- function(TC, lnK, n_hyd) {
  m <- sat_molality(TC + 273.15, lnK, n_hyd)
  if (is.na(m)) NA_real_ else 100*aw(TC + 273.15, m)
}
rh_mirabilite   <- function(TC) rh_deliq(TC, lnK_10, 10)
rh_thenardite   <- function(TC) rh_deliq(TC, lnK_V, 0)
rh_heptahydrate <- function(TC) rh_deliq(TC, lnK_7, 7)
rh_then_mir     <- function(TC) {                  # 水合/脫水邊界
  T <- TC + 273.15
  100*exp((lnK_10(T) - lnK_V(T))/10)
}

# 三相點 (兩條藍線交點)
Ttp <- uniroot(function(tc) rh_mirabilite(tc) - rh_then_mir(tc), c(25, 34))$root
RHtp <- rh_then_mir(Ttp)
cat(sprintf("三相點: T = %.2f C, RH = %.2f %%\n", Ttp, RHtp))
cat(sprintf("驗證: 20C mirabilite 潮解 RH = %.2f %% (文獻 95.6%%)\n", rh_mirabilite(20)))

# ---------------- 建立邊界內插表 ----------------
gT   <- seq(0, 50, by = 0.05)
gMIR <- sapply(gT, function(t) if (t <= Ttp) rh_mirabilite(t) else rh_thenardite(t))
gLOW <- sapply(gT, function(t) if (t <= Ttp) rh_then_mir(t)  else rh_thenardite(t))

classify <- function(TC, RH) {
  up  <- approx(gT, gMIR, xout = TC, rule = 2)$y
  low <- approx(gT, gLOW, xout = TC, rule = 2)$y
  ifelse(RH >= up, "solution", ifelse(RH < low, "thenardite", "mirabilite"))
}

# =====================================================================
# 讀取 IoT 資料 (依你原本的工作流程)
# =====================================================================
library(readxl); library(dplyr); library(lubridate)

xlsx <- "D:/IoT/T_vs_RH_phase_diagram/1黃土水_水牛群像/1_Huang_water_buffalo_20210809_20220922.xlsx"
sheet <- excel_sheets(xlsx)
IoT <- bind_rows(lapply(setNames(sheet, sheet),
                        function(x) read_xlsx(xlsx, sheet = x, skip = 7)),
                 .id = "Sheet")

# 解析中文 上午/下午 時間字串 (不依賴系統 locale)
parse_tw <- function(s) {
  pm <- grepl("下午", s)
  s2 <- gsub("[ ]*(上午|下午)[ ]*", " ", s)
  t  <- as.POSIXct(s2, format = "%Y/%m/%d %H:%M:%S", tz = "Asia/Taipei")
  h  <- as.integer(format(t, "%H"))
  t  <- t + ifelse(pm & h < 12, 12*3600, 0)      # 下午 1~11 點 +12h
  t  - ifelse(!pm & h == 12, 12*3600, 0)          # 上午 12 點 = 00 時
}
IoT$datetime <- parse_tw(IoT$DataRecordTime)

analyse <- function(Tvec, RHvec, dt, label) {
  ok <- !is.na(Tvec) & !is.na(RHvec) & !is.na(dt) & RHvec > 0 & RHvec <= 100
  d  <- data.frame(datetime = dt[ok], T = Tvec[ok], RH = RHvec[ok]) |>
        arrange(datetime)
  d$phase <- classify(d$T, d$RH)
  cat("\n===", label, "===\n"); print(table(d$phase))

  # 跨線事件: 相區改變的時刻
  ch <- which(d$phase[-1] != d$phase[-nrow(d)])
  ev <- data.frame(time_before = d$datetime[ch],   time_after = d$datetime[ch+1],
                   from_phase  = d$phase[ch],      to_phase   = d$phase[ch+1],
                   T_after = d$T[ch+1],            RH_after   = d$RH[ch+1])
  cat("跨線事件總數:", nrow(ev), "\n")
  cat("  潮解 (→solution):",  sum(ev$to_phase   == "solution"), "\n")
  cat("  再結晶 (solution→):", sum(ev$from_phase == "solution"), "\n")
  print(ev |> mutate(month = format(time_after, "%Y-%m")) |> count(month))
  write.csv(ev, paste0("crossing_events_", label, "_R.csv"), row.names = FALSE)

  # 相圖
  png(paste0("phase_diagram_", label, "_R.png"), width = 900, height = 620)
  cols <- c(solution = "red", mirabilite = "green3", thenardite = "orange")
  plot(d$T, d$RH, col = cols[d$phase], pch = 20, cex = .4,
       xlim = c(0, 50), ylim = c(15, 100),
       main = paste("2021-2022 Huang Water Buffalo —", label),
       xlab = expression(paste("Temperature ( ", degree, "C )")),
       ylab = "Relative Humidity (%)")
  tt <- seq(0, Ttp, length.out = 200); tt2 <- seq(Ttp, 50, length.out = 80)
  lines(tt,  sapply(tt,  rh_mirabilite), col = "blue", lwd = 2)
  lines(tt,  sapply(tt,  rh_then_mir),   col = "blue", lwd = 2, lty = 2)
  lines(tt2, sapply(tt2, rh_thenardite), col = "blue", lwd = 2)
  abline(h = 75, col = "orange", lwd = 2)
  text(40, 78, "NaCl"); text(8, 50, "Thenardite"); text(10, 82, "Mirabilite")
  text(35, 96, expression("Na"[2]*"SO"[4]*" solution"))
  legend("bottomright", legend = names(cols), col = cols, pch = 19, cex = .8)
  dev.off()
  invisible(ev)
}

ev_micro <- analyse(IoT$`微環境溫度`, IoT$`微環境濕度`, IoT$datetime, "micro")
ev_peri  <- analyse(IoT$`周邊溫度`,   IoT$`周邊濕度`,   IoT$datetime, "peripheral")