Revert "use corrective term in SI formulae for GetSatVapPres to erase…

… the gap at 0 C (#26)" (#29) This reverts commit 223ef02. This issue will be handled by a new PR addressing the discontinuity in the NR method in GetTDewPointFromVapPres.
psychrometrics · Mar 19, 2019 · 91dd95f · 91dd95f
1 parent 223ef02
commit 91dd95f
Show file tree

Hide file tree

Showing 5 changed files with 10 additions and 43 deletions.
diff --git a/src/c/psychrolib.c b/src/c/psychrolib.c
@@ -647,16 +647,11 @@ 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())
   {
@@ -679,10 +674,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) + CORRECTIVE_TERM_SI;
+          + 2.0747825E-09 * pow(T, 3) - 9.484024E-13 * pow(T, 4) + 4.1635019 * log(T);
     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) - CORRECTIVE_TERM_SI;
+        - 1.4452093E-08 * pow(T, 3) + 6.5459673 * log(T);
     else
       return INVALID;             // TDryBulb is out of range [-100, 200]
   }

diff --git a/src/fortran/psychrolib.f90 b/src/fortran/psychrolib.f90
@@ -878,11 +878,6 @@ 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]
@@ -892,8 +887,6 @@ 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
@@ -919,10 +912,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) + CORRECTIVE_TERM_SI
+                  + 2.0747825E-09 * T**3 - 9.484024E-13 * T**4 + 4.1635019 * log(T)
         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) - CORRECTIVE_TERM_SI
+                  - 1.4452093E-08 * T**3 + 6.5459673 * log(T)
         end if
       end if
 

diff --git a/src/js/psychrolib.js b/src/js/psychrolib.js
@@ -605,15 +605,10 @@ 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())
     {
@@ -638,10 +633,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) + CORRECTIVE_TERM_SI;
+                + 2.0747825E-09 * pow(T, 3) - 9.484024E-13 * pow(T, 4) + 4.1635019 * log(T);
       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) - CORRECTIVE_TERM_SI;
+                - 1.4452093E-08 * pow(T, 3) + 6.5459673 * log(T);
       else
         return INVALID;             // TDryBulb is out of range [-100, 200]
     }

diff --git a/src/python/psychrolib.py b/src/python/psychrolib.py
@@ -871,14 +871,7 @@ 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")
@@ -899,12 +892,10 @@ 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) \
-                  + CORRECTIVE_TERM_SI
+                  + 2.0747825E-09 * math.pow(T, 3) - 9.484024E-13 * math.pow(T, 4) + 4.1635019 * math.log(T)
         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) \
-                  - CORRECTIVE_TERM_SI
+                  - 1.4452093E-08 * math.pow(T, 3) + 6.5459673 * math.log(T)
 
     SatVapPres = math.exp(LnPws)
     return SatVapPres

diff --git a/src/vba/psychrolib.bas b/src/vba/psychrolib.bas
@@ -1153,15 +1153,8 @@ 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
 
@@ -1191,10 +1184,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) + CORRECTIVE_TERM_SI
+              + 2.0747825E-09 * T ^ 3 - 9.484024E-13 * T ^ 4 + 4.1635019 * Log(T)
     Else
         LnPws = -5800.2206 / T + 1.3914993 - 0.048640239 * T + 0.000041764768 * T ^ 2 _
-              - 0.000000014452093 * T ^ 3 + 6.5459673 * Log(T) - CORRECTIVE_TERM_SI
+              - 0.000000014452093 * T ^ 3 + 6.5459673 * Log(T)
     End If
   End If