Excel VBA to Write SWCC Equation

Question

I'm tring to make Satayanaga's best fit equation for unimodal 2017 to VBA, because I keep having problem with Excel's solver if I put the equation straight to the cell. Satyanaga's 2017 SWCC equation

I'm trying with this. But I'm not sure if it's right, because it such a complicated one for me. It'll be great if someone can confirm this or give a suggestion.

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 Double
    Dim beta As Double
    Dim term1 As Double, term2 As Double
    Dim erfcPart As Double
    Dim insideLog1 As Double, insideLog2 As Double
    
    ' --- Define ß ---
    If psi <= psi_a Then
        beta = 0
    Else
        beta = 1
    End If
    
    ' --- Check log arguments ---
    insideLog1 = 1 + psi / psi_r
    insideLog2 = 1 + (10 ^ 6) / psi_r
    
    If insideLog1 <= 0 Or insideLog2 <= 0 Then
        Satyanaga = CVErr(xlErrNum) ' return #NUM! error in Excel
        Exit Function
    End If
    
    ' --- First log fraction (always between 0 and 1) ---
    term1 = 1 - (Log(insideLog1) / Log(insideLog2))
    
    ' --- ERFC part ---
    If (psi_a - psi_m) = 0 Then
        erfcPart = 0
    Else
        On Error Resume Next
        erfcPart = Application.WorksheetFunction.ErfC(Log((psi_a - psi) / (psi_a - psi_m)) / s)
        If Err.Number <> 0 Then
            Satyanaga = CVErr(xlErrNum)
            Exit Function
        End If
        On Error GoTo 0
    End If
    
    ' --- Second bracket ---
    term2 = theta_r + (theta_s - theta_r) * (1 - beta * erfcPart)
    
    ' --- Final result ---
    Satyanaga = term1 * term2
End Function

Does the function return wrong values? Are the values general wrong or only on some value ranges of the input parameters? — Shrotter
– Shrotter, Commented Aug 18 at 5:50
Please edit the question and provide enough code to create a Minimal Reproducible Example we need to be able to copy your code and try it ourselves. — MT1
– MT1, Commented Aug 18 at 6:46
We need some representative test data to be able to test your code. Please edit your example input data into the question. davetallett26.github.io/excel-markdown.html — MT1
– MT1, Commented Aug 18 at 6:48
The code looks alright. What are you having problems with exactly? You could validate the logic like this: find some input and output that are known to be valid and check against your implementation — Michał Turczyn
– Michał Turczyn, Commented Aug 18 at 7:37
You need to provide a sample test harness around it that shows typical input parameters and range of phi. The equation is a little odd to my eyes. Use of erfc to blend in the more complex behaviour looks to me like an engineering hack rather than a physical equation. BTW when beta=0 the equation simplifies a lot and you should avoid evaluating the nasty complicated term Log((psi_a-psi)/(psi_a-psi_m)) (which is guaranteed to fail if it evaluates correctly when psi>psi_a). When debugging avoid on error resume next it makes code behaviour hard to understand. — Martin Brown
– Martin Brown, Commented Aug 18 at 8:17

Martin Brown · Accepted Answer · 2025-08-18 14:03:37Z

I think your problem is mainly with the peculiar way that the original Satyanaga equation has been kludged together by the author(s). If you try to take the log of a negative value your code will return an error value. This failure mode is guaranteed to happen for one of phi<=phi_a or psi>=psi_a with your function as originally written. Once a Nan always a Nan - even when you multiply through by zero (when beta = 0).

Log needs protecting from invalid negative arguments. You had it right in the first section although I would hope that pressures in this problem are always positive.

This is a quick rewrite that avoids the obvious pitfalls with log and assumes that the conditions for valid psi specified in the original equation definition are correct.

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 Double
   Dim term1 As Double, term2 As Double
   Dim erfcPart As Double
   Dim insideLog1 As Double, insideLog2 As Double

' --- Check log arguments ---
   insideLog1 = 1 + psi / psi_r
   insideLog2 = 1 + (10 ^ 6) / psi_r

   If insideLog1 <= 0 Or insideLog2 <= 0 Then
       Satyanaga = CVErr(xlErrNum) ' return #NUM! error in Excel
       Exit Function
   End If

' --- First log fraction (always between 0 and 1) ---
   term1 = 1 - (Log(insideLog1) / Log(insideLog2))
    
' --- Act on ß ---
   If psi<=psi_a Then '' beta = 0 so simplify accordingly
      term2 = theta_s
   Else
' --- ERFC part --- now only evaluated for valid conditions psi > psi_a
' NB this means numerator is negative so I presume denominator is too
      If (psi_a - psi_m) = 0 Then
         erfcPart = 0
      Else
         erfcPart = Application.WorksheetFunction.ErfC(Log((psi_a - psi) / (psi_a - psi_m)) / s)
      End If

      ' --- Second bracket ---
      term2 = theta_r + (theta_s - theta_r) * (1 - erfcPart)
   End If
' --- Final result ---
   Satyanaga = term1 * term2
End Function

You may still find that Excel solver doesn't like this strange function when asked to optimise it but at least this variant should be able to execute and return a numeric value without any errors for typical values of psi.

Collectives™ on Stack Overflow

Excel VBA to Write SWCC Equation

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related