将Satyanaga 2017单峰SWCC拟合方程转为Excel VBA的验证请求

代码正确性验证 #

对照Satyanaga 2017年提出的单峰土壤水分特征曲线（SWCC）拟合方程，你的VBA代码核心逻辑完全符合方程定义：

1. β参数定义：当ψ ≤ ψₐ时β=0，否则β=1，与方程规则一致。
2. 第一部分项（term1）：使用VBA自然对数Log()计算1 - ln(1+ψ/ψᵣ)/ln(1+10⁶/ψᵣ)，匹配方程要求（VBA的Log对应数学中的ln）。
3. 互补误差函数部分（erfcPart）：调用Excel工作表函数ErfC()实现方程中的互补误差函数，参数ln((ψₐ-ψ)/(ψₐ-ψₘ))/s的计算逻辑正确。
4. 第二部分项（term2）：θᵣ + (θₛ-θᵣ)(1 - β·erfcPart)的计算完全贴合方程结构。
5. 错误处理：提前检查对数参数有效性、避免除以0错误、捕获ErfC计算异常并返回#NUM!，覆盖了主要异常场景。

整体代码可以正确实现Satyanaga 2017年的SWCC方程计算。

改进建议 #

针对代码的可读性、健壮性和规范性，给出以下优化点：

1. 修正函数返回值类型 #

原函数声明为Double，但返回了错误值CVErr(xlErrNum)，严格来说Double类型无法存储错误值，建议改为Variant类型，避免潜在类型冲突：

Function Satyanaga(psi As Double, theta_s As Double, theta_r As Double, _
                   psi_a As Double, psi_m As Double, s As Double, psi_r As Double) As Variant

2. 细化异常判断，减少错误捕获依赖 #

原代码通过On Error Resume Next捕获Log异常，可提前判断(ψₐ-ψ)/(ψₐ-ψₘ)是否大于0，主动避免对数定义域错误：

ratio = (psi_a - psi) / (psi_a - psi_m)
If ratio <= 0 Then
    Satyanaga = CVErr(xlErrNum)
    Exit Function
End If

3. 常量与参数优化 #

• 将方程中的10^6定义为命名常量，提升代码可读性和可维护性：

  Const PSI_MAX As Double = 10 ^ 6
  

• 参数psi_r与theta_r命名易混淆，建议重命名为psi_entry这类更清晰的名称，避免使用时出错。

4. 简化代码写法 #

用IIf函数简化β参数的判断逻辑，让代码更简洁：

beta = IIf(psi <= psi_a, 0, 1)

改进后的完整代码 #

Function Satyanaga(psi As Double, theta_s As Double, theta_r As Double, _
                   psi_a As Double, psi_m As Double, s As Double, psi_entry As Double) As Variant
    Const PSI_MAX As Double = 10 ^ 6 ' 方程中用于归一化的最大基质势常量
    
    Dim beta As Double
    Dim term1 As Double, term2 As Double
    Dim erfcPart As Double
    Dim insideLog1 As Double, insideLog2 As Double
    Dim ratio As Double
    
    ' 定义β参数
    beta = IIf(psi <= psi_a, 0, 1)
    
    ' 检查对数参数有效性，避免定义域错误
    insideLog1 = 1 + psi / psi_entry
    insideLog2 = 1 + PSI_MAX / psi_entry
    If insideLog1 <= 0 Or insideLog2 <= 0 Then
        Satyanaga = CVErr(xlErrNum)
        Exit Function
    End If
    
    ' 计算第一部分项
    term1 = 1 - (Log(insideLog1) / Log(insideLog2))
    
    ' 计算互补误差函数部分
    If (psi_a - psi_m) = 0 Then
        erfcPart = 0
    Else
        ratio = (psi_a - psi) / (psi_a - psi_m)
        ' 提前检查对数输入的有效性
        If ratio <= 0 Then
            Satyanaga = CVErr(xlErrNum)
            Exit Function
        End If
        ' 计算ErfC，保留错误捕获处理极端情况
        On Error Resume Next
        erfcPart = Application.WorksheetFunction.ErfC(Log(ratio) / s)
        If Err.Number <> 0 Then
            Satyanaga = CVErr(xlErrNum)
            Exit Function
        End If
        On Error GoTo 0
    End If
    
    ' 计算第二部分项
    term2 = theta_r + (theta_s - theta_r) * (1 - beta * erfcPart)
    
    ' 返回最终结果
    Satyanaga = term1 * term2
End Function

内容的提问来源于stack exchange，提问作者Una Qunn