From 223ef028370c13f8e9d55dfeed43fcfb290b9f45 Mon Sep 17 00:00:00 2001 From: dmey Date: Sat, 16 Mar 2019 12:33:53 +0000 Subject: [PATCH] use corrective term in SI formulae for GetSatVapPres to erase the gap at 0 C (#26) --- src/c/psychrolib.c | 9 +++++++-- src/fortran/psychrolib.f90 | 11 +++++++++-- src/js/psychrolib.js | 9 +++++++-- src/python/psychrolib.py | 13 +++++++++++-- src/vba/psychrolib.bas | 11 +++++++++-- 5 files changed, 43 insertions(+), 10 deletions(-) diff --git a/src/c/psychrolib.c b/src/c/psychrolib.c index 2060e561..e18cea14 100644 --- a/src/c/psychrolib.c +++ b/src/c/psychrolib.c @@ -647,11 +647,16 @@ double GetHumRatioFromEnthalpyAndTDryBulb // (o) Humidity ratio in lb_H₂O lb_ // Return saturation vapor pressure given dry-bulb temperature. // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 & 6 +// Notes: the SI formulae show a discontinuity at 0 C. In rare cases this discontinuity creates issues +// in GetTDewPointFromVapPres. To avoid the problem, a small corrective term is added/subtracted +// to the ASHRAE formulae to make the formulae continuous at 0 C. The effect on the results is +// negligible (0.005%), well below the accuracy of the formulae double GetSatVapPres // (o) Vapor Pressure of saturated air in Psi [IP] or Pa [SI] ( double TDryBulb // (i) Dry bulb temperature in °F [IP] or °C [SI] ) { double LnPws, T; + double CORRECTIVE_TERM_SI = 4.851e-05; // Small corrective term to make the function continuous at 0 C. if (isIP()) { @@ -674,10 +679,10 @@ double GetSatVapPres // (o) Vapor Pressure of saturated air in Psi [I T = GetTKelvinFromTCelsius(TDryBulb); if (TDryBulb >= -100. && TDryBulb <= 0.) LnPws = -5.6745359E+03 / T + 6.3925247 - 9.677843E-03 * T + 6.2215701E-07 * T * T - + 2.0747825E-09 * pow(T, 3) - 9.484024E-13 * pow(T, 4) + 4.1635019 * log(T); + + 2.0747825E-09 * pow(T, 3) - 9.484024E-13 * pow(T, 4) + 4.1635019 * log(T) + CORRECTIVE_TERM_SI; else if (TDryBulb > 0. && TDryBulb <= 200.) LnPws = -5.8002206E+03 / T + 1.3914993 - 4.8640239E-02 * T + 4.1764768E-05 * T * T - - 1.4452093E-08 * pow(T, 3) + 6.5459673 * log(T); + - 1.4452093E-08 * pow(T, 3) + 6.5459673 * log(T) - CORRECTIVE_TERM_SI; else return INVALID; // TDryBulb is out of range [-100, 200] } diff --git a/src/fortran/psychrolib.f90 b/src/fortran/psychrolib.f90 index f7293bc0..ce19e588 100644 --- a/src/fortran/psychrolib.f90 +++ b/src/fortran/psychrolib.f90 @@ -878,6 +878,11 @@ function GetSatVapPres(TDryBulb) result(SatVapPres) !+ Return saturation vapor pressure given dry-bulb temperature. !+ Reference: !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 5 + !+ Notes: + !+ The SI formulae show a discontinuity at 0 C. In rare cases this discontinuity creates issues + !+ in GetTDewPointFromVapPres. To avoid the problem, a small corrective term is added/subtracted + !+ to the ASHRAE formulae to make the formulae continuous at 0 C. The effect on the results is + !+ negligible (0.005%), well below the accuracy of the formulae real, intent(in) :: TDryBulb !+ Dry-bulb temperature in °F [IP] or °C [SI] @@ -887,6 +892,8 @@ function GetSatVapPres(TDryBulb) result(SatVapPres) !+ Log of Vapor Pressure of saturated air (dimensionless) real :: T !+ Dry bulb temperature in R [IP] or K [SI] + real, parameter :: CORRECTIVE_TERM_SI = 4.851e-05 + !+ Small corrective term to make the function continuous at 0 C. if (isIP()) then if (TDryBulb < -148.0 .or. TDryBulb > 392.0) then @@ -912,10 +919,10 @@ function GetSatVapPres(TDryBulb) result(SatVapPres) if (TDryBulb <= 0) then LnPws = -5.6745359E+03 / T + 6.3925247 - 9.677843E-03 * T + 6.2215701E-07 * T**2 & - + 2.0747825E-09 * T**3 - 9.484024E-13 * T**4 + 4.1635019 * log(T) + + 2.0747825E-09 * T**3 - 9.484024E-13 * T**4 + 4.1635019 * log(T) + CORRECTIVE_TERM_SI else LnPws = -5.8002206E+03 / T + 1.3914993 - 4.8640239E-02 * T + 4.1764768E-05 * T**2 & - - 1.4452093E-08 * T**3 + 6.5459673 * log(T) + - 1.4452093E-08 * T**3 + 6.5459673 * log(T) - CORRECTIVE_TERM_SI end if end if diff --git a/src/js/psychrolib.js b/src/js/psychrolib.js index 757dc622..37679a4b 100644 --- a/src/js/psychrolib.js +++ b/src/js/psychrolib.js @@ -605,10 +605,15 @@ function Psychrometrics() { // Return saturation vapor pressure given dry-bulb temperature. // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 & 6 + // Notes: the SI formulae show a discontinuity at 0 C. In rare cases this discontinuity creates issues + // in GetTDewPointFromVapPres. To avoid the problem, a small corrective term is added/subtracted + // to the ASHRAE formulae to make the formulae continuous at 0 C. The effect on the results is + // negligible (0.005%), well below the accuracy of the formulae this.GetSatVapPres = function // (o) Vapor Pressure of saturated air in Psi [IP] or Pa [SI] ( TDryBulb // (i) Dry bulb temperature in °F [IP] or °C [SI] ) { var LnPws, T; + var CORRECTIVE_TERM_SI = 4.851e-05; // Small corrective term to make the function continuous at 0 C. if (this.isIP()) { @@ -633,10 +638,10 @@ function Psychrometrics() { T = this.GetTKelvinFromTCelsius(TDryBulb); if (TDryBulb >= -100. && TDryBulb <= 0.) LnPws = -5.6745359E+03 / T + 6.3925247 - 9.677843E-03 * T + 6.2215701E-07 * T * T - + 2.0747825E-09 * pow(T, 3) - 9.484024E-13 * pow(T, 4) + 4.1635019 * log(T); + + 2.0747825E-09 * pow(T, 3) - 9.484024E-13 * pow(T, 4) + 4.1635019 * log(T) + CORRECTIVE_TERM_SI; else if (TDryBulb > 0. && TDryBulb <= 200.) LnPws = -5.8002206E+03 / T + 1.3914993 - 4.8640239E-02 * T + 4.1764768E-05 * T * T - - 1.4452093E-08 * pow(T, 3) + 6.5459673 * log(T); + - 1.4452093E-08 * pow(T, 3) + 6.5459673 * log(T) - CORRECTIVE_TERM_SI; else return INVALID; // TDryBulb is out of range [-100, 200] } diff --git a/src/python/psychrolib.py b/src/python/psychrolib.py index 0a9ff266..59e041f1 100644 --- a/src/python/psychrolib.py +++ b/src/python/psychrolib.py @@ -871,7 +871,14 @@ def GetSatVapPres(TDryBulb: float) -> float: Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 5 & 6 + Notes: + The SI formulae show a discontinuity at 0 C. In rare cases this discontinuity creates issues + in GetTDewPointFromVapPres. To avoid the problem, a small corrective term is added/subtracted + to the ASHRAE formulae to make the formulae continuous at 0 C. The effect on the results is + negligible (0.005%), well below the accuracy of the formulae """ + CORRECTIVE_TERM_SI = 4.851e-05 # small corrective term to make the function continuous at 0 C. + if isIP(): if (TDryBulb < -148 or TDryBulb > 392): raise ValueError("Dry bulb temperature must be in range [-148, 392]°F") @@ -892,10 +899,12 @@ def GetSatVapPres(TDryBulb: float) -> float: if (TDryBulb <= 0): LnPws = -5.6745359E+03 / T + 6.3925247 - 9.677843E-03 * T + 6.2215701E-07 * T**2 \ - + 2.0747825E-09 * math.pow(T, 3) - 9.484024E-13 * math.pow(T, 4) + 4.1635019 * math.log(T) + + 2.0747825E-09 * math.pow(T, 3) - 9.484024E-13 * math.pow(T, 4) + 4.1635019 * math.log(T) \ + + CORRECTIVE_TERM_SI else: LnPws = -5.8002206E+03 / T + 1.3914993 - 4.8640239E-02 * T + 4.1764768E-05 * T**2 \ - - 1.4452093E-08 * math.pow(T, 3) + 6.5459673 * math.log(T) + - 1.4452093E-08 * math.pow(T, 3) + 6.5459673 * math.log(T) \ + - CORRECTIVE_TERM_SI SatVapPres = math.exp(LnPws) return SatVapPres diff --git a/src/vba/psychrolib.bas b/src/vba/psychrolib.bas index c75ffd62..89dc64f1 100644 --- a/src/vba/psychrolib.bas +++ b/src/vba/psychrolib.bas @@ -1153,8 +1153,15 @@ Function GetSatVapPres(ByVal TDryBulb As Variant) As Variant ' ' Reference: ' ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 5 & 6 +' +' Notes: +' The SI formulae show a discontinuity at 0 C. In rare cases this discontinuity creates issues +' in GetTDewPointFromVapPres. To avoid the problem, a small corrective term is added/subtracted +' to the ASHRAE formulae to make the formulae continuous at 0 C. The effect on the results is +' negligible (0.005%), well below the accuracy of the formulae ' Dim LnPws As Variant, T As Variant + Const CORRECTIVE_TERM_SI = 4.851e-05 ' Small corrective term to make the function continuous at 0 C. On Error GoTo ErrHandler @@ -1184,10 +1191,10 @@ Function GetSatVapPres(ByVal TDryBulb As Variant) As Variant If (TDryBulb <= 0) Then LnPws = -5674.5359 / T + 6.3925247 - 0.009677843 * T + 0.00000062215701 * T ^ 2 _ - + 2.0747825E-09 * T ^ 3 - 9.484024E-13 * T ^ 4 + 4.1635019 * Log(T) + + 2.0747825E-09 * T ^ 3 - 9.484024E-13 * T ^ 4 + 4.1635019 * Log(T) + CORRECTIVE_TERM_SI Else LnPws = -5800.2206 / T + 1.3914993 - 0.048640239 * T + 0.000041764768 * T ^ 2 _ - - 0.000000014452093 * T ^ 3 + 6.5459673 * Log(T) + - 0.000000014452093 * T ^ 3 + 6.5459673 * Log(T) - CORRECTIVE_TERM_SI End If End If