The VBA listing below calculates some of the steam table functions using the 1967 ASME formulation. Note that it calculates in only one direction>> Pressure and Temperature (psia, degF) to the steam properties. It will not calculate "in the other direction." You must iterate (or use Excel's Solver) to do that.
Hope it helps.
Dim qDataCache(100), Cnst(100)
'____________________________ Specific Volume _______________________________ Function fnVsh(P, F) If Cnst(1) = 0 Then Initialize B = P / 3208.234759 'Beta (ratio of P/Pcrit) Th = (F + 459.67) / 1165.14 'Theta (ratio of degR/degRcrit) Chi = 4.260321148 * Th / B ' Initialize Chi (ratio of v/vCrit)
Kount1 = 7 For M = 1 To 5 Call Sum2(Sum, Kount1, 0) Chi = Chi  Sum * M * B ^ (M  1) Next
Kount1 = 36: Kount2 = 59 For M = 6 To 8 Call Sum2(Sum, Kount1, 0) Top = Sum * (M  2) * B ^ (1  M) Call Sum2(Sum, Kount2, 0) Bot = (Sum + B ^ (2  M)) ^ 2 Chi = Chi  Top / Bot Next
Call Sum3(Sum, Bl, 0) Chi = Chi + 11 * Sum * (B / Bl) ^ 10 fnVsh = Chi * 0.0507785289 End Function '____________________________________________________________________________
'______________________________ Entropy _____________________________________ Function fnSsh(P, F) ' Private B, Sig, Sum, Kount1, Kount2, Kount3, Kount4, Top, Bot, Bl, M If Cnst(1) = 0 Then Initialize B = P / 3208.234759 'Beta (ratio of P/Pcrit) Th = (F + 459.67) / 1165.14 'Theta (ratio of degR/degRcrit) Sig = 4.260321148 * Log(B) + Cnst(1) * Log(Th) ' Initialize Sig ' which is Sigma Call Sum1(Sum, 0) ' s/sCrit Sig = Sig  Sum
Kount1 = 7 For M = 1 To 5 Call Sum2(Sum, Kount1, 1) Sig = Sig  Sum * bLil * B ^ M Next
Kount1 = 59: Kount2 = 36: Kount3 = 59: Kount4 = 36 For M = 6 To 8 Call Sum2(Top, Kount1, 1) Call Sum2(Sum, Kount2, 0) Top = Top + Sum Call Sum2(Sum, Kount3, 0) Bot = Sum + B ^ (2  M) Call Sum2(Sum, Kount4, 1) Sig = Sig  bLil * (Sum  Top / Bot) / Bot Next
Call Sum3(Sum, Bl, 1) Sig = Sig + B * Sum * (B / Bl) ^ 10 fnSsh = Sig * 0.02587358228 End Function '____________________________________________________________________________
'_________________________________ Enthalpy _________________________________ Function fnHsh(P, F) 'Private B, Eps, Sum, Kount1, Kount2, Kount3, Kount4, Top, Bot, Buv, Bl, M If Cnst(1) = 0 Then Initialize B = P / 3208.234759 'Beta (ratio of P/Pcrit) Th = (F + 459.67) / 1165.14 'Theta (ratio of degR/degRcrit) Eps = Cnst(1) * Th ' Initialize Eps which is Epselon (h/hCrit) Call Sum1(Sum, 1) Eps = Eps  Sum
Kount1 = 7: Kount2 = 7 For M = 1 To 5 Call Sum2(Sum, Kount1, 0) Top = Sum Call Sum2(Sum, Kount2, 1) Eps = Eps  (Top + bLil * Th * Sum) * B ^ M Next
Kount1 = 59: Kount2 = 59: Kount3 = 36: Kount4 = 36 For M = 6 To 8 Call Sum2(Top, Kount1, 1) Call Sum2(Sum, Kount2, 0) Bot = Sum + B ^ (2  M) Call Sum2(Buv, Kount3, 0) Eps = Eps  (Buv + bLil * Th * Sum  Buv * bLil * Th * Top / Bot) / Bot Next
Call Sum3(Sum, Bl, 2) Eps = Eps + Sum * B * (B / Bl) ^ 10 fnHsh = Eps * 30.14634566 End Function '____________________________________________________________________________
'______________________________ Sat Pressure ________________________________ Function fnPSat(F) If F < 32 Then Beep: qPrint "Bad Data": Exit Function B = P / 3208.234759 'Beta (ratio of P/Pcrit) Th = (F + 459.67) / 1165.14 'Theta (ratio of degR/degRcrit) temp = 1  Th Top = temp * (temp * (temp * (temp * (temp * 118.9646225 + 64.23285504) _  168.1706546)  26.08023696)  7.691234564) Bot = 1 + temp * (temp * 20.9750676 + 4.16711732) fnPSat = 3208.234759 * Exp(Top / (Th * Bot)  temp / (1000000000# * temp ^ 2 + 6)) End Function '____________________________________________________________________________
'_______________________________ Sat Temperature ____________________________ Function fnTSat(P) 'Fit by Bize Not ASME Steam Tables If P < 0.0886 Then Beep: qPrint "Bad Data : exit function" temp = 268.633  0.0054 * P + 2.4417 * Sqr(P) + 32.558 * Sqr(Sqr(P)) _ + 50.338 * Sqr(Sqr(Sqr(P))) + 285.04 * Sqr(Sqr(Sqr(Sqr(P)))) If P > 1400 Then Del = P  1400 temp = temp  Del * (0.000106 + Del * (0.000000343 + Del * 0.00000000013)) End If fnTSat = temp End Function '____________________________________________________________________________
Sub Initialize() ' Reads ASME constants into Cnst() and bLil
'B sub 0,v ===> Cnst(1) thru (6) qData 16.83599274, 28.56067796, 54.38923329 qData 0.4330662834, 0.6547711697, 0.08565182058
'B sub u,v and z sub u,v where u = 1 to 5 '========> Cnst(7) thru (35) (zero is flag to end loop) qData 0.06670375918, 13, 1.388983801, 3, 0 qData 0.08390104328, 18, 0.02614670893, 2, 0.03373439453, 1, 0 qData 0.4520918904, 18, 0.1069036614, 10, 0 qData 0.5975336707, 25, 0.08847535804, 14, 0 qData 0.5958051609, 32, 0.5159303373, 28, 0.2075021122, 24, 0
'B sub u,v and z sub u,v where u = 6 to 8 '========> Cnst(36) thru (50) (zero is flag to end loop) qData 0.1190610271, 12, 0.09867174132, 11, 0 qData 0.1683998803, 24, 0.05809438001, 18, 0 qData 0.006552390126, 24, 0.0005710218649, 14, 0
'B sub 9,v =======> Cnst(51) thru (57) bLil = Cnst(58) qData 193.6587558, 1388.522425, 4126.607219, 6508.211677 qData 5745.984054, 2693.088365, 523.5718623, 0.7633333333 ' bLil
'b sub u,lambda and x sub u,lambda '=====> Cnst(59) thru (69) (zero is flag to end loop) qData 0.4006073948, 14, 0 qData 0.08636081627, 19, 0 qData 0.8532322921, 54, 0.3460208861, 27, 0
For I = 1 To 69 qRead Cnst(I) Next I bLil = Cnst(58) End Sub
Sub NonDim(Psia, degF, Beta, Theta) ' Changes Psia and degF into Nondimensional parameters ' Beta is ratio of pressure/pressure at critical point of steam ' Theta is ratio of degR/degR at critical point
Beta = Psia / 3208.234759 Theta = (degF + 459.67) / 1165.14
End Sub
Sub Sum1(Sum1, iFlag) 'Sum1 is sum of (N1*Beta(0,v)*Theta^N2) for v = 1 to 5 ' where if ' iFlag = 1 then N1 = v1 ' N2 = v2 ' else N1 = v2 ' N2 = v1 'Private LilV, N1, N2 Sum1 = 0 For LilV = 1 To 5 If iFlag Then N1 = LilV  2: N2 = LilV  1 Else N1 = LilV  1: N2 = LilV  2 End If Sum1 = Sum1 + N1 * Cnst(LilV + 1) * Th ^ N2 Next
End Sub
Sub Sum2(Sum2, K, iFlag) 'Private Prod
Sum2 = 0 While Cnst(K) <> 0 Prod = Cnst(K) * Exp(Cnst(K + 1) * bLil * (1  Th)) If iFlag Then Prod = Prod * Cnst(K + 1) Sum2 = Sum2 + Prod K = K + 2 Wend K = K + 1 End Sub
Sub Sum3(Sum3, Bl, iFlag) 'Private LilV, Prod, temp, L1, L2, Blprime L1 = 34.17061978: L2 = 19.31380707 Bl = 15.74373327 + Th * (L1 + Th * L2) Blprime = L1 + 2 * L2 * Th
Sum3 = 0 For LilV = 0 To 6 Prod = Cnst(51 + LilV) * Exp(LilV * bLil * (1  Th)) If iFlag <> 0 Then temp = 10 * Blprime / Bl + LilV * bLil If iFlag = 1 Then Prod = Prod * temp If iFlag = 2 Then Prod = Prod * (1 + Th * temp) End If Sum3 = Sum3 + Prod Next End Sub
