渗流理论核心组件 (Perkp) 代码解析

perkp.py 是整个仿真系统的“动力引擎”。它根据化学键断裂状态（配位数、桥键概率等），通过组合数学计算出煤结构的崩溃程度及产物分布。

🔬 物理数学背景：渗流理论

该代码实现了基于 Bethe 晶格（Bethe Lattice）的渗流统计：

$p$ ($l + c$)：桥键存活的总概率。
$p^*$：无穷集群（Infinite Network）存在的临界概率解。
$Q_n$：大小为 $n$ 的片段在统计意义上的出现频率。

⚙️ 函数逻辑实现：`perkp(...)`

这是一个高性能的数值求解函数，涵盖了从非线性迭代到组合数学计算的全过程。

1. $p^*$ 的牛顿迭代求解

当系统超过渗流阈值时，代码通过牛顿迭代法求解非线性方程，以确定可溶物（Sol fraction）的比例。

# 核心修复：高精度牛顿迭代
const_term = p * (1.0 - p) ** sig_minus_1
for _ in range(25):
    om_ps = 1.0 - pstar
    f = pstar * om_ps ** sig_minus_1 - const_term
    fp = om_ps ** sig_minus_1 - pstar * sig_minus_1 * om_ps ** sig_minus_2
    ppstar = pstar - f / fp # 更新迭代值
    if abs(1.0 - ppstar / pstar) <= 1e-4: break

2. 质量平衡逻辑（对齐 MATLAB）

代码严格遵循质量守恒定律，将煤的总量划分为：

fgas：桥键断裂产生的轻质气体。
ftart_out：理论上可转化为焦油的单体总量。
fchar：保持在交联网络中的焦炭部分。

3. 基于 Gamma 函数的分布计算

为了计算焦油的分子量分布，代码使用 lgamma 函数处理组合数，避免了大数阶乘导致的溢出：

# 使用对数 Gamma 避免大数溢出
ln_fgam = lgamma(xm + 1.0) - lgamma(yk + 1.0) - lgamma(xm - yk + 1.0)
qn = (sig_plus_1 / xm * exp(ln_fgam)) * (p ** yk) * (one_minus_p ** tn) / n

ft_raw: 存储了从单体到阶齐聚物的产率分布。
mt_raw: 计算每一阶片段对应的物理分子量。

⚡ 高性能优化特性

@njit(cache=True): 由于涉及大量的指数运算和 lgamma 调用，JIT 编译将此函数的运行速度提升了数十倍。
reshape(-1, 1): 输出结果自动转化为列向量，无缝对接后端的 NumPy 矩阵运算，避免了维度坍缩产生的 Bug。
预计算常数: 如 sig_minus_1, sig_plus_1 等在循环外计算，减少冗余浮点运算。

📋 关键参数对照表

变量名	代码实现	物理/数学意义
`sig` ()	`sig`	晶格配位数 (Coordination Number)
`y[0]` ()	`l_val`	桥键 (Bridges) 比例
`y[2]` ()	`c_val`	侧链 (Cobs) 比例
`ma` ()	`ma`	单体单元分子量
`rba`	`rba`	桥键与单体的质量比参数