diff --git a/src/c/psychrolib.c b/src/c/psychrolib.c
index 5e768efc..9230135e 100644
--- a/src/c/psychrolib.c
+++ b/src/c/psychrolib.c
@@ -53,16 +53,24 @@
  * Global constants
  *****************************************************************************************************/
 
-#define R_DA_IP 53.350            // Universal gas constant for dry air (IP version) in ft∙lbf/lb_da/R
-                                  // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1
-#define R_DA_SI 287.042           // Universal gas constant for dry air (SI version) in J/kg_da/K
-                                  // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1
-#define INVALID -99999            // Invalid value
+#define R_DA_IP 53.350                  // Universal gas constant for dry air (IP version) in ft∙lbf/lb_da/R.
+                                        // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1.
+#define R_DA_SI 287.042                 // Universal gas constant for dry air (SI version) in J/kg_da/K.
+                                        // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1.
+#define INVALID -99999                  // Invalid value.
 
-#define MAX_ITER_COUNT 100        // Maximum number of iterations before exiting while loops.
+#define MAX_ITER_COUNT 100              // Maximum number of iterations before exiting while loops.
 
-#define MIN_HUM_RATIO 1e-7        // Minimum acceptable humidity ratio used/returned by any functions.
-                                  // Any value above 0 or below the MIN_HUM_RATIO will be reset to this value.
+#define MIN_HUM_RATIO 1e-7              // Minimum acceptable humidity ratio used/returned by any functions.
+                                        // Any value above 0 or below the MIN_HUM_RATIO will be reset to this value.
+
+#define FREEZING_POINT_WATER_IP 32.0    // Freezing point of water in Fahrenheit.
+
+#define FREEZING_POINT_WATER_SI 0.0     // Freezing point of water in Celsius.
+
+#define TRIPLE_POINT_WATER_IP 32.018    // Triple point of water in Fahrenheit.
+
+#define TRIPLE_POINT_WATER_SI 0.01      // Triple point of water in Celsius.
 
 
 /******************************************************************************************************
@@ -303,7 +311,7 @@ double dLnPws_        // (o) Derivative of natural log of vapor pressure of satu
   {
     T = GetTRankineFromTFahrenheit(TDryBulb);
 
-    if (TDryBulb <= 32.)
+    if (TDryBulb <= TRIPLE_POINT_WATER_IP)
       dLnPws = 1.0214165E+04 / pow(T, 2) - 5.3765794E-03 + 2 * 1.9202377E-07 * T
                + 2 * 3.5575832E-10 * pow(T, 2) - 4 * 9.0344688E-14 * pow(T, 3) + 4.1635019 / T;
     else
@@ -314,7 +322,7 @@ double dLnPws_        // (o) Derivative of natural log of vapor pressure of satu
   {
     T = GetTKelvinFromTCelsius(TDryBulb);
 
-    if (TDryBulb <= 0.)
+    if (TDryBulb <= TRIPLE_POINT_WATER_SI)
       dLnPws = 5.6745359E+03 / pow(T, 2) - 9.677843E-03 + 2 * 6.2215701E-07 * T
                + 3 * 2.0747825E-09 * pow(T, 2) - 4 * 9.484024E-13 * pow(T, 3) + 4.1635019 / T;
     else
@@ -342,40 +350,21 @@ double GetTDewPointFromVapPres  // (o) Dew Point temperature in °F [IP] or °C
 {
   // Bounds function of the system of units
   double BOUNDS[2];              // Domain of validity of the equations
-  double T_WATER_FREEZE, T_WATER_FREEZE_LOW, T_WATER_FREEZE_HIGH, PWS_FREEZE_LOW, PWS_FREEZE_HIGH;
 
   if (isIP())
   {
     BOUNDS[0] = -148.;
     BOUNDS[1] = 392.;
-    T_WATER_FREEZE = 32.;
   }
   else
   {
     BOUNDS[0] = -100.;
     BOUNDS[1] = 200.;
-    T_WATER_FREEZE = 0.;
   }
 
   // Bounds outside which a solution cannot be found
-  ASSERT (VapPres >= GetSatVapPres(BOUNDS[0]) && VapPres <= GetSatVapPres(BOUNDS[1]), "Partial pressure of water vapor is outside range of validity of equations")
-
-  // Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing
-  T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.;          // Temperature just below freezing
-  T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE / 10.;         // Temperature just above freezing
-  PWS_FREEZE_LOW = GetSatVapPres(T_WATER_FREEZE_LOW);
-  PWS_FREEZE_HIGH = GetSatVapPres(T_WATER_FREEZE_HIGH);
-
-  // Restrict iteration to either left or right part of the saturation vapor pressure curve
-  // to avoid iterating back and forth across the discontinuity of the curve at the freezing point
-  // When the partial pressure of water vapor is within the discontinuity of GetSatVapPres,
-  // simply return the freezing point of water.
-  if (VapPres < PWS_FREEZE_LOW)
-    BOUNDS[1] = T_WATER_FREEZE_LOW;
-  else if (VapPres > PWS_FREEZE_LOW)
-    BOUNDS[0] = T_WATER_FREEZE_HIGH;
-  else
-    return T_WATER_FREEZE;
+  ASSERT (VapPres >= GetSatVapPres(BOUNDS[0]) && VapPres <= GetSatVapPres(BOUNDS[1]),
+          "Partial pressure of water vapor is outside range of validity of equations")
 
   // We use NR to approximate the solution.
   // First guess
@@ -484,7 +473,7 @@ double GetHumRatioFromTWetBulb  // (o) Humidity Ratio in lb_H₂O lb_Air⁻¹ [I
 
   if (isIP())
   {
-    if (TWetBulb >= 32.)
+    if (TWetBulb >= FREEZING_POINT_WATER_IP)
       HumRatio = ((1093. - 0.556 * TWetBulb) * Wsstar - 0.240 * (TDryBulb - TWetBulb))
       / (1093. + 0.444 * TDryBulb - TWetBulb);
     else
@@ -493,7 +482,7 @@ double GetHumRatioFromTWetBulb  // (o) Humidity Ratio in lb_H₂O lb_Air⁻¹ [I
   }
   else
   {
-    if (TWetBulb >= 0.)
+    if (TWetBulb >= FREEZING_POINT_WATER_SI)
       HumRatio = ((2501. - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb))
          / (2501. + 1.86 * TDryBulb - 4.186 * TWetBulb);
     else
@@ -731,6 +720,12 @@ 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
+// Important note: the ASHRAE formulae are defined above and below the freezing point but have
+// a discontinuity at the freezing point. This is a small inaccuracy on ASHRAE's part: the formulae
+// should be defined above and below the triple point of water (not the feezing point) in which case
+// the discontinuity vanishes. It is essential to use the triple point of water otherwise function
+// GetTDewPointFromVapPres, which inverts the present function, does not converge properly around
+// the freezing point.
 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]
   )
@@ -742,28 +737,26 @@ double GetSatVapPres            // (o) Vapor Pressure of saturated air in Psi [I
     ASSERT(TDryBulb >= -148. && TDryBulb <= 392., "Dry bulb temperature is outside range [-148, 392]")
 
     T = GetTRankineFromTFahrenheit(TDryBulb);
-    if (TDryBulb >= -148. && TDryBulb <= 32.)
+
+    if (TDryBulb <= TRIPLE_POINT_WATER_IP)
       LnPws = (-1.0214165E+04 / T - 4.8932428 - 5.3765794E-03 * T + 1.9202377E-07 * T * T
         + 3.5575832E-10 * pow(T, 3) - 9.0344688E-14 * pow(T, 4) + 4.1635019 * log(T));
-    else if (TDryBulb > 32. && TDryBulb <= 392.)
+    else
       LnPws = -1.0440397E+04 / T - 1.1294650E+01 - 2.7022355E-02 * T + 1.2890360E-05 * T * T
       - 2.4780681E-09 * pow(T, 3) + 6.5459673 * log(T);
-    else
-      return INVALID;             // TDryBulb is out of range [-148, 392]
   }
   else
   {
     ASSERT(TDryBulb >= -100. && TDryBulb <= 200., "Dry bulb temperature is outside range [-100, 200]")
 
     T = GetTKelvinFromTCelsius(TDryBulb);
-    if (TDryBulb >= -100. && TDryBulb <= 0.)
+
+    if (TDryBulb <= TRIPLE_POINT_WATER_SI)
       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);
-    else if (TDryBulb > 0. && TDryBulb <= 200.)
+    else
       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);
-    else
-      return INVALID;             // TDryBulb is out of range [-100, 200]
   }
 
   return exp(LnPws);
diff --git a/src/fortran/psychrolib.f90 b/src/fortran/psychrolib.f90
index c4300c6e..e5824fd5 100644
--- a/src/fortran/psychrolib.f90
+++ b/src/fortran/psychrolib.f90
@@ -96,19 +96,19 @@ module psychrolib
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   real, parameter ::  R_DA_IP = 53.350
-    !+ Universal gas constant for dry air (IP version) in ft lb_Force lb_DryAir⁻¹ R⁻¹
+    !+ Universal gas constant for dry air (IP version) in ft lb_Force lb_DryAir⁻¹ R⁻¹.
     !+ Reference:
-    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1
+    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1.
 
   real, parameter ::  R_DA_SI = 287.042
-    !+ Universal gas constant for dry air (SI version) in J kg_DryAir⁻¹ K⁻¹
+    !+ Universal gas constant for dry air (SI version) in J kg_DryAir⁻¹ K⁻¹.
     !+ Reference:
-    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1
+    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1.
 
   integer, parameter :: IP = 1
   integer, parameter :: SI = 2
 
-  integer  :: PSYCHROLIB_UNITS = 0 ! 0 = undefined
+  integer  :: PSYCHROLIB_UNITS = 0 ! 0 = undefined.
     !+ Unit system to use.
 
   real ::  PSYCHROLIB_TOLERANCE = 1.0
@@ -121,6 +121,18 @@ module psychrolib
     !+ Minimum acceptable humidity ratio used/returned by any functions.
     !+ Any value above 0 or below the MIN_HUM_RATIO will be reset to this value.
 
+  real, parameter  :: FREEZING_POINT_WATER_IP = 32.0
+    !+ float: Freezing point of water in Fahrenheit.
+
+  real, parameter  :: FREEZING_POINT_WATER_SI = 0.0
+    !+ float: Freezing point of water in Celsius.
+
+  real, parameter  :: TRIPLE_POINT_WATER_IP = 32.018
+    !+ float: Triple point of water in Fahrenheit.
+
+  real, parameter  :: TRIPLE_POINT_WATER_SI = 0.01
+    !+ float: Triple point of water in Celsius.
+
 
   contains
 
@@ -408,7 +420,7 @@ function dLnPws_(TDryBulb) result(dLnPws)
 
       T = GetTRankineFromTFahrenheit(TDryBulb)
 
-      if (TDryBulb <= 32.) then
+      if (TDryBulb <= TRIPLE_POINT_WATER_IP) then
         dLnPws = 1.0214165E+04 / T**2 - 5.3765794E-03 + 2 * 1.9202377E-07 * T &
                  + 2 * 3.5575832E-10 * T**2 - 4 * 9.0344688E-14 * T**3 + 4.1635019 / T
       else
@@ -420,7 +432,7 @@ function dLnPws_(TDryBulb) result(dLnPws)
 
       T = GetTKelvinFromTCelsius(TDryBulb)
 
-      if (TDryBulb <= 0) then
+      if (TDryBulb <= TRIPLE_POINT_WATER_SI) then
         dLnPws = 5.6745359E+03 / T**2 - 9.677843E-03 + 2 * 6.2215701E-07 * T &
                  + 3 * 2.0747825E-09 * T**2 - 4 * 9.484024E-13 * T**3 + 4.1635019 / T
       else
@@ -462,17 +474,14 @@ function GetTDewPointFromVapPres(TDryBulb, VapPres) result(TDewPoint)
       !+ Valid temperature range in °F [IP] or °C [SI]
     integer             :: index
       !+ Index used in the calculation
-     real ::  T_WATER_FREEZE, T_WATER_FREEZE_LOW, T_WATER_FREEZE_HIGH, PWS_FREEZE_LOW, PWS_FREEZE_HIGH
 
     ! Bounds and step size as a function of the system of units
     if (isIP()) then
         BOUNDS(1) = -148.0
         BOUNDS(2) =  392.0
-        T_WATER_FREEZE = 32.
     else
         BOUNDS(1) = -100.0
         BOUNDS(2) =  200.0
-        T_WATER_FREEZE = 0.
     end if
 
     ! Validity check -- bounds outside which a solution cannot be found
@@ -480,26 +489,6 @@ function GetTDewPointFromVapPres(TDryBulb, VapPres) result(TDewPoint)
       error stop "Error: partial pressure of water vapor is outside range of validity of equations"
     end if
 
-    ! Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing
-    T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.        ! Temperature just below freezing
-    T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.        ! Temperature just below freezing
-    T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE             ! Temperature just above freezing
-    PWS_FREEZE_LOW = GetSatVapPres(T_WATER_FREEZE_LOW)
-    PWS_FREEZE_HIGH = GetSatVapPres(T_WATER_FREEZE_HIGH)
-
-    ! Restrict iteration to either left or right part of the saturation vapor pressure curve
-    ! to avoid iterating back and forth across the discontinuity of the curve at the freezing point
-    ! When the partial pressure of water vapor is within the discontinuity of GetSatVapPres,
-    ! simply return the freezing point of water.
-    if (VapPres < PWS_FREEZE_LOW) then
-        BOUNDS(2) = T_WATER_FREEZE_LOW
-    else if (VapPres > PWS_FREEZE_HIGH) then
-        BOUNDS(1) = T_WATER_FREEZE_HIGH
-    else
-        TDewPoint = T_WATER_FREEZE
-        return
-    end if
-
     ! We use NR to approximate the solution.
     TDewPoint = TDryBulb
     lnVP = log(VapPres)
@@ -635,7 +624,7 @@ function GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) result(HumRatio)
     Wsstar = GetSatHumRatio(TWetBulb, Pressure)
 
     if (isIP()) then
-      if (TWetBulb >= 32.0) then
+      if (TWetBulb >= FREEZING_POINT_WATER_IP) then
         HumRatio = ((1093.0 - 0.556 * TWetBulb) * Wsstar - 0.240 * (TDryBulb - TWetBulb))   &
                    / (1093.0 + 0.444 * TDryBulb - TWetBulb)
       else
@@ -643,7 +632,7 @@ function GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) result(HumRatio)
                    / (1220.0 + 0.444 * TDryBulb - 0.48 * TWetBulb)
       end if
     else
-      if (TWetBulb >= 0.0) then
+      if (TWetBulb >= FREEZING_POINT_WATER_SI) then
         HumRatio = ((2501.0 - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb))   &
                    / (2501.0 + 1.86 * TDryBulb - 4.186 * TWetBulb)
       else
@@ -970,6 +959,12 @@ function GetSatVapPres(TDryBulb) result(SatVapPres)
     !+ Return saturation vapor pressure given dry-bulb temperature.
     !+ Reference:
     !+ ASHRAE Handbook - Fundamentals (2017) ch. 1  eqn 5
+    !+ Important note: the ASHRAE formulae are defined above and below the freezing point but have
+    !+ a discontinuity at the freezing point. This is a small inaccuracy on ASHRAE's part: the formulae
+    !+ should be defined above and below the triple point of water (not the feezing point) in which case 
+    !+ the discontinuity vanishes. It is essential to use the triple point of water otherwise function
+    !+ GetTDewPointFromVapPres, which inverts the present function, does not converge properly around
+    !+ the freezing point.
 
     real, intent(in)  ::  TDryBulb
       !+ Dry-bulb temperature in °F [IP] or °C [SI]
@@ -987,7 +982,7 @@ function GetSatVapPres(TDryBulb) result(SatVapPres)
 
       T = GetTRankineFromTFahrenheit(TDryBulb)
 
-      if (TDryBulb <= 32.) then
+      if (TDryBulb <= TRIPLE_POINT_WATER_IP) then
         LnPws = (-1.0214165E+04 / T - 4.8932428 - 5.3765794E-03 * T + 1.9202377E-07 * T**2    &
                 + 3.5575832E-10 * T**3 - 9.0344688E-14 * T**4 + 4.1635019 * log(T))
       else
@@ -1002,7 +997,7 @@ function GetSatVapPres(TDryBulb) result(SatVapPres)
 
         T = GetTKelvinFromTCelsius(TDryBulb)
 
-        if (TDryBulb <= 0) then
+        if (TDryBulb <= TRIPLE_POINT_WATER_SI) 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)
         else
diff --git a/src/js/psychrolib.js b/src/js/psychrolib.js
index fd1b3f98..81e6f0b1 100644
--- a/src/js/psychrolib.js
+++ b/src/js/psychrolib.js
@@ -57,16 +57,25 @@ function Psychrometrics() {
    * Global constants
    *****************************************************************************************************/
 
-  var R_DA_IP = 53.350;         // Universal gas constant for dry air (IP version) in ft lb_Force lb_DryAir⁻¹ R⁻¹
-                                // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1
-  var R_DA_SI = 287.042;        // Universal gas constant for dry air (SI version) in J kg_DryAir⁻¹ K⁻¹
-                                // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1
-  var INVALID = -99999;         // Invalid value (dimensionless)
+  var R_DA_IP = 53.350;               // Universal gas constant for dry air (IP version) in ft lb_Force lb_DryAir⁻¹ R⁻¹.
+                                      // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1.
+  var R_DA_SI = 287.042;              // Universal gas constant for dry air (SI version) in J kg_DryAir⁻¹ K⁻¹.
+                                      // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1.
+  var INVALID = -99999;               // Invalid value (dimensionless).
 
-  var MAX_ITER_COUNT = 100      // Maximum number of iterations before exiting while loops.
+  var MAX_ITER_COUNT = 100            // Maximum number of iterations before exiting while loops.
+
+  var MIN_HUM_RATIO = 1e-7            // Minimum acceptable humidity ratio used/returned by any functions.
+                                      // Any value above 0 or below the MIN_HUM_RATIO will be reset to this value.
+
+  var FREEZING_POINT_WATER_IP = 32.0  // Freezing point of water in Fahrenheit.
+
+  var FREEZING_POINT_WATER_SI = 0.0   // Freezing point of water in Celsius.
+
+  var TRIPLE_POINT_WATER_IP = 32.018  // Triple point of water in Fahrenheit.
+
+  var TRIPLE_POINT_WATER_SI = 0.01    // Triple point of water in Celsius.
 
-  var MIN_HUM_RATIO = 1e-7      // Minimum acceptable humidity ratio used/returned by any functions.
-                                // Any value above 0 or below the MIN_HUM_RATIO will be reset to this value.
 
 
   /******************************************************************************************************
@@ -268,7 +277,7 @@ function Psychrometrics() {
     {
       T = this.GetTRankineFromTFahrenheit(TDryBulb);
 
-      if (TDryBulb <= 32.)
+      if (TDryBulb <= TRIPLE_POINT_WATER_IP)
         dLnPws = 1.0214165E+04 / pow(T, 2) - 5.3765794E-03 + 2 * 1.9202377E-07 * T
                  + 2 * 3.5575832E-10 * pow(T, 2) - 4 * 9.0344688E-14 * pow(T, 3) + 4.1635019 / T;
       else
@@ -279,7 +288,7 @@ function Psychrometrics() {
     {
       T = this.GetTKelvinFromTCelsius(TDryBulb);
 
-      if (TDryBulb <= 0.)
+      if (TDryBulb <= TRIPLE_POINT_WATER_SI)
         dLnPws = 5.6745359E+03 / pow(T, 2) - 9.677843E-03 + 2 * 6.2215701E-07 * T
                  + 3 * 2.0747825E-09 * pow(T, 2) - 4 * 9.484024E-13 * pow(T, 3) + 4.1635019 / T;
       else
@@ -306,40 +315,20 @@ function Psychrometrics() {
     ) {
    // Bounds function of the system of units
   var BOUNDS              // Domain of validity of the equations
-  var T_WATER_FREEZE, T_WATER_FREEZE_LOW, T_WATER_FREEZE_HIGH, PWS_FREEZE_LOW, PWS_FREEZE_HIGH;
 
   if (this.isIP())
   {
     BOUNDS = [-148., 392.];   // Domain of validity of the equations
-    T_WATER_FREEZE = 32.;
   }
   else
   {
     BOUNDS = [-100., 200.];   // Domain of validity of the equations
-    T_WATER_FREEZE = 0.;
   }
 
   // Bounds outside which a solution cannot be found
   if (VapPres < this.GetSatVapPres(BOUNDS[0]) || VapPres > this.GetSatVapPres(BOUNDS[1]))
     throw new Error("Partial pressure of water vapor is outside range of validity of equations");
 
-  // Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing
-  T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.;          // Temperature just below freezing
-  T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE / 10.;         // Temperature just above freezing
-  PWS_FREEZE_LOW = this.GetSatVapPres(T_WATER_FREEZE_LOW);
-  PWS_FREEZE_HIGH = this.GetSatVapPres(T_WATER_FREEZE_HIGH);
-
-  // Restrict iteration to either left or right part of the saturation vapor pressure curve
-  // to avoid iterating back and forth across the discontinuity of the curve at the freezing point
-  // When the partial pressure of water vapor is within the discontinuity of GetSatVapPres,
-  // simply return the freezing point of water.
-  if (VapPres < PWS_FREEZE_LOW)
-    BOUNDS[1] = T_WATER_FREEZE_LOW;
-  else if (VapPres > PWS_FREEZE_HIGH)
-    BOUNDS[0] = T_WATER_FREEZE_HIGH;
-  else
-    return T_WATER_FREEZE;
-
   // We use NR to approximate the solution.
   // First guess
   var TDewPoint = TDryBulb;      // Calculated value of dew point temperatures, solved for iteratively in °F [IP] or °C [SI]
@@ -447,7 +436,7 @@ function Psychrometrics() {
 
       if (this.isIP())
       {
-        if (TWetBulb >= 32.)
+        if (TWetBulb >= FREEZING_POINT_WATER_IP)
           HumRatio = ((1093. - 0.556 * TWetBulb) * Wsstar - 0.240 * (TDryBulb - TWetBulb))
           / (1093. + 0.444 * TDryBulb - TWetBulb);
         else
@@ -456,7 +445,7 @@ function Psychrometrics() {
       }
       else
       {
-        if (TWetBulb >= 0.)
+        if (TWetBulb >= FREEZING_POINT_WATER_SI)
           HumRatio = ((2501. - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb))
              / (2501. + 1.86 * TDryBulb - 4.186 * TWetBulb);
         else
@@ -687,6 +676,12 @@ function Psychrometrics() {
 
   // Return saturation vapor pressure given dry-bulb temperature.
   // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 & 6
+  // Important note: the ASHRAE formulae are defined above and below the freezing point but have
+  // a discontinuity at the freezing point. This is a small inaccuracy on ASHRAE's part: the formulae
+  // should be defined above and below the triple point of water (not the feezing point) in which case
+  // the discontinuity vanishes. It is essential to use the triple point of water otherwise function
+  // GetTDewPointFromVapPres, which inverts the present function, does not converge properly around
+  // the freezing point.
   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]
     ) {
@@ -698,14 +693,12 @@ function Psychrometrics() {
         throw new Error("Dry bulb temperature is outside range [-148, 392]");
 
       T = this.GetTRankineFromTFahrenheit(TDryBulb);
-      if (TDryBulb >= -148. && TDryBulb <= 32.)
+      if (TDryBulb <= TRIPLE_POINT_WATER_IP)
         LnPws = (-1.0214165E+04 / T - 4.8932428 - 5.3765794E-03 * T + 1.9202377E-07 * T * T
                 + 3.5575832E-10 * pow(T, 3) - 9.0344688E-14 * pow(T, 4) + 4.1635019 * log(T));
-      else if (TDryBulb > 32. && TDryBulb <= 392.)
+      else
         LnPws = -1.0440397E+04 / T - 1.1294650E+01 - 2.7022355E-02 * T + 1.2890360E-05 * T * T
                 - 2.4780681E-09 * pow(T, 3) + 6.5459673 * log(T);
-      else
-        return INVALID;             // TDryBulb is out of range [-148, 392]
     }
     else
     {
@@ -713,14 +706,12 @@ function Psychrometrics() {
         throw new Error("Dry bulb temperature is outside range [-100, 200]");
 
       T = this.GetTKelvinFromTCelsius(TDryBulb);
-      if (TDryBulb >= -100. && TDryBulb <= 0.)
+      if (TDryBulb <= TRIPLE_POINT_WATER_SI)
         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);
-      else if (TDryBulb > 0. && TDryBulb <= 200.)
+      else
         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);
-      else
-        return INVALID;             // TDryBulb is out of range [-100, 200]
     }
 
     return exp(LnPws);
diff --git a/src/python/psychrolib.py b/src/python/psychrolib.py
index 4cd247fe..2beb8898 100644
--- a/src/python/psychrolib.py
+++ b/src/python/psychrolib.py
@@ -84,6 +84,25 @@
 
 """
 
+FREEZING_POINT_WATER_IP = 32.0
+"""float: Freezing point of water in Fahrenheit.
+
+"""
+
+FREEZING_POINT_WATER_SI = 0.0
+"""float: Freezing point of water in Celsius.
+
+"""
+
+TRIPLE_POINT_WATER_IP = 32.018
+"""float: Triple point of water in Fahrenheit.
+
+"""
+
+TRIPLE_POINT_WATER_SI = 0.01
+"""float: Triple point of water in Celsius.
+
+"""
 
 #######################################################################################################
 # Helper functions
@@ -403,7 +422,7 @@ def dLnPws_(TDryBulb: float) -> float:
     """
     if isIP():
         T = GetTRankineFromTFahrenheit(TDryBulb)
-        if TDryBulb <= 32.:
+        if TDryBulb <= TRIPLE_POINT_WATER_IP:
             dLnPws = 1.0214165E+04 / math.pow(T, 2) - 5.3765794E-03 + 2 * 1.9202377E-07 * T \
                   + 2 * 3.5575832E-10 * math.pow(T, 2) - 4 * 9.0344688E-14 * math.pow(T, 3) + 4.1635019 / T
         else:
@@ -411,7 +430,7 @@ def dLnPws_(TDryBulb: float) -> float:
                   - 3 * 2.4780681E-09 * math.pow(T, 2) + 6.5459673 / T
     else:
         T = GetTKelvinFromTCelsius(TDryBulb)
-        if TDryBulb <= 0.:
+        if TDryBulb <= TRIPLE_POINT_WATER_SI:
             dLnPws = 5.6745359E+03 / math.pow(T, 2) - 9.677843E-03 + 2 * 6.2215701E-07 * T \
                   + 3 * 2.0747825E-09 * math.pow(T, 2) - 4 * 9.484024E-13 * math.pow(T, 3) + 4.1635019 / T
         else:
@@ -441,38 +460,19 @@ def GetTDewPointFromVapPres(TDryBulb: float, VapPres: float) -> float:
         narrower range of validity.
         The Newton-Raphson (NR) method is used on the logarithm of water vapour
         pressure as a function of temperature, which is a very smooth function
-        Convergence is usually achieved in 3 to 5 iterations.
+        Convergence is usually achieved in 3 to 5 iterations. 
         TDryBulb is not really needed here, just used for convenience.
 
     """
     if isIP():
         BOUNDS = [-148, 392]
-        T_WATER_FREEZE = 32.
     else:
         BOUNDS = [-100, 200]
-        T_WATER_FREEZE = 0.
 
     # Validity check -- bounds outside which a solution cannot be found
     if VapPres < GetSatVapPres(BOUNDS[0]) or VapPres > GetSatVapPres(BOUNDS[1]):
         raise ValueError("Partial pressure of water vapor is outside range of validity of equations")
 
-    # Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing
-    T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.          # Temperature just below freezing
-    T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE / 10.         # Temperature just above freezing
-    PWS_FREEZE_LOW = GetSatVapPres(T_WATER_FREEZE_LOW)
-    PWS_FREEZE_HIGH = GetSatVapPres(T_WATER_FREEZE_HIGH)
-
-    # Restrict iteration to either left or right part of the saturation vapor pressure curve
-    # to avoid iterating back and forth across the discontinuity of the curve at the freezing point
-    # When the partial pressure of water vapor is within the discontinuity of GetSatVapPres,
-    # simply return the freezing point of water.
-    if (VapPres < PWS_FREEZE_LOW):
-        BOUNDS[1] = T_WATER_FREEZE_LOW
-    elif (VapPres > PWS_FREEZE_HIGH):
-        BOUNDS[0] = T_WATER_FREEZE_HIGH
-    else:
-        return T_WATER_FREEZE
-
     # We use NR to approximate the solution.
     # First guess
     TDewPoint = TDryBulb        # Calculated value of dew point temperatures, solved for iteratively
@@ -596,14 +596,14 @@ def GetHumRatioFromTWetBulb(TDryBulb: float, TWetBulb: float, Pressure: float) -
     Wsstar = GetSatHumRatio(TWetBulb, Pressure)
 
     if isIP():
-       if TWetBulb >= 32:
+       if TWetBulb >= FREEZING_POINT_WATER_IP:
            HumRatio = ((1093 - 0.556 * TWetBulb) * Wsstar - 0.240 * (TDryBulb - TWetBulb)) \
                     / (1093 + 0.444 * TDryBulb - TWetBulb)
        else:
            HumRatio = ((1220 - 0.04 * TWetBulb) * Wsstar - 0.240 * (TDryBulb - TWetBulb)) \
                     / (1220 + 0.444 * TDryBulb - 0.48*TWetBulb)
     else:
-       if TWetBulb >= 0:
+       if TWetBulb >= FREEZING_POINT_WATER_SI:
            HumRatio = ((2501. - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) \
                     / (2501. + 1.86 * TDryBulb - 4.186 * TWetBulb)
        else:
@@ -949,6 +949,12 @@ def GetSatVapPres(TDryBulb: float) -> float:
 
     Reference:
         ASHRAE Handbook - Fundamentals (2017) ch. 1  eqn 5 & 6
+        Important note: the ASHRAE formulae are defined above and below the freezing point but have
+        a discontinuity at the freezing point. This is a small inaccuracy on ASHRAE's part: the formulae
+        should be defined above and below the triple point of water (not the feezing point) in which case
+        the discontinuity vanishes. It is essential to use the triple point of water otherwise function
+        GetTDewPointFromVapPres, which inverts the present function, does not converge properly around
+        the freezing point.
 
     """
     if isIP():
@@ -957,7 +963,7 @@ def GetSatVapPres(TDryBulb: float) -> float:
 
         T = GetTRankineFromTFahrenheit(TDryBulb)
 
-        if (TDryBulb <= 32.):
+        if (TDryBulb <= TRIPLE_POINT_WATER_IP):
             LnPws = (-1.0214165E+04 / T - 4.8932428 - 5.3765794E-03 * T + 1.9202377E-07 * T**2 \
                   + 3.5575832E-10 * math.pow(T, 3) - 9.0344688E-14 * math.pow(T, 4) + 4.1635019 * math.log(T))
         else:
@@ -969,7 +975,7 @@ def GetSatVapPres(TDryBulb: float) -> float:
 
         T = GetTKelvinFromTCelsius(TDryBulb)
 
-        if (TDryBulb <= 0):
+        if (TDryBulb <= TRIPLE_POINT_WATER_SI):
             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)
         else:
diff --git a/src/vba/psychrolib.bas b/src/vba/psychrolib.bas
index 9c2ff810..bf075608 100644
--- a/src/vba/psychrolib.bas
+++ b/src/vba/psychrolib.bas
@@ -62,12 +62,22 @@ End Enum
 ' Global constants
 '******************************************************************************************************
 
-Private Const R_DA_IP = 53.35               ' Universal gas constant for dry air (IP version) in ft lbf/lb_DryAir/R
-Private Const R_DA_SI = 287.042             ' Universal gas constant for dry air (SI version) in J/kg_DryAir/K
-Private Const MAX_ITER_COUNT = 100         ' Maximum number of iterations before exiting while loops.
-Private Const MIN_HUM_RATIO = 1e-7          ' Minimum acceptable humidity ratio used/returned by any functions.
-                                            ' Any value above 0 or below the MIN_HUM_RATIO will be reset to this value.
+Private Const R_DA_IP = 53.35                 ' Universal gas constant for dry air (IP version) in ft lbf/lb_DryAir/R.
 
+Private Const R_DA_SI = 287.042               ' Universal gas constant for dry air (SI version) in J/kg_DryAir/K.
+
+Private Const MAX_ITER_COUNT = 100            ' Maximum number of iterations before exiting while loops.
+
+Private Const MIN_HUM_RATIO = 1e-7            ' Minimum acceptable humidity ratio used/returned by any functions.
+                                              ' Any value above 0 or below the MIN_HUM_RATIO will be reset to this value.
+
+Private Const FREEZING_POINT_WATER_IP = 32.0  ' Freezing point of water, in °F
+
+Private Const FREEZING_POINT_WATER_SI = 0.0   ' Freezing point of water, in °C
+
+Private Const TRIPLE_POINT_WATER_IP = 32.018  ' Triple point of water, in °F
+
+Private Const TRIPLE_POINT_WATER_SI = 0.01    ' Triple point of water, in °C
 
 '******************************************************************************************************
 ' Helper functions
@@ -502,7 +512,7 @@ Private Function dLnPws_(TDryBulb As Variant) As Variant
   Dim T As Variant
   If (isIP()) Then
     T = GetTRankineFromTFahrenheit(TDryBulb)
-    If (TDryBulb <= 32.) Then
+    If (TDryBulb <= TRIPLE_POINT_WATER_IP) Then
       dLnPws_ = 10214.165 / T ^ 2 - 0.0053765794 + 2 * 0.00000019202377 * T _
              + 2 * 3.5575832E-10 * T ^ 2 - 4 * 9.0344688E-14 * T ^ 3 + 4.1635019 / T
     Else
@@ -511,7 +521,7 @@ Private Function dLnPws_(TDryBulb As Variant) As Variant
     End If
   Else
     T = GetTKelvinFromTCelsius(TDryBulb)
-    If (TDryBulb <= 0.) Then
+    If (TDryBulb <= TRIPLE_POINT_WATER_SI) Then
       dLnPws_ = 5674.5359 / T ^ 2 - 0.009677843 + 2 * 0.00000062215701 * T _
              + 3 * 2.0747825E-09 * T ^ 2 - 4 * 9.484024E-13 * T ^ 3 + 4.1635019 / T
     Else
@@ -546,18 +556,14 @@ Function GetTDewPointFromVapPres(ByVal TDryBulb As Variant, ByVal VapPres As Var
 '        TDryBulb is not really needed here, just used for convenience.
 '
   Dim BOUNDS(2) As Variant
-  Dim T_WATER_FREEZE As Variant, T_WATER_FREEZE_LOW As Variant, T_WATER_FREEZE_HIGH As Variant
-  Dim PWS_FREEZE_LOW As Variant, PWS_FREEZE_HIGH As Variant
   Dim PSYCHROLIB_TOLERANCE As Variant
 
   If (isIP()) Then
     BOUNDS(1) = -148.
     BOUNDS(2) = 392.
-    T_WATER_FREEZE = 32.
   Else
     BOUNDS(1) = -100.
     BOUNDS(2) = 200.
-    T_WATER_FREEZE = 0.
   End If
 
   On Error GoTo ErrHandler
@@ -568,24 +574,6 @@ Function GetTDewPointFromVapPres(ByVal TDryBulb As Variant, ByVal VapPres As Var
   End If
 
   PSYCHROLIB_TOLERANCE = GetTol()
-  ' Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing
-  T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.        ' Temperature just below freezing
-  T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE / 10.        ' Temperature just above freezing
-  PWS_FREEZE_LOW = GetSatVapPres(T_WATER_FREEZE_LOW)
-  PWS_FREEZE_HIGH = GetSatVapPres(T_WATER_FREEZE_HIGH)
-
-  ' Restrict iteration to either left or right part of the saturation vapor pressure curve
-  ' to avoid iterating back and forth across the discontinuity of the curve at the freezing point
-  ' When the partial pressure of water vapor is within the discontinuity of GetSatVapPres,
-  ' simply return the freezing point of water.
-  If (VapPres < PWS_FREEZE_LOW) Then
-    BOUNDS(2) = T_WATER_FREEZE_LOW
-  ElseIf (VapPres > PWS_FREEZE_HIGH) Then
-    BOUNDS(1) = T_WATER_FREEZE_HIGH
-  Else
-    GetTDewPointFromVapPres = T_WATER_FREEZE
-    Exit Function
-  End If
 
   Dim TDewPoint As Variant
   Dim lnVP As Variant
@@ -600,7 +588,7 @@ Function GetTDewPointFromVapPres(ByVal TDryBulb As Variant, ByVal VapPres As Var
   TDewPoint = TDryBulb        ' Calculated value of dew point temperatures, solved for iteratively
   lnVP = Log(VapPres)         ' Partial pressure of water vapor in moist air
 
-
+  ' Iteration
   Do
     TDewPoint_iter = TDewPoint   ' Value of Tdp used in NR calculation
     lnVP_iter = Log(GetSatVapPres(TDewPoint_iter))
@@ -752,13 +740,13 @@ Function GetHumRatioFromTWetBulb(ByVal TDryBulb As Variant, ByVal TWetBulb As Va
   End If
 
   If isIP() Then
-    If (TWetBulb >= 32) Then
+    If (TWetBulb >= FREEZING_POINT_WATER_IP) Then
       HumRatio = ((1093 - 0.556 * TWetBulb) * Wsstar - 0.24 * (TDryBulb - TWetBulb)) / (1093 + 0.444 * TDryBulb - TWetBulb)
     Else
       HumRatio = ((1220 - 0.04 * TWetBulb) * Wsstar - 0.24 * (TDryBulb - TWetBulb)) / (1220 + 0.444 * TDryBulb - 0.48 * TWetBulb)
     End If
   Else
-    If (TWetBulb >= 0) Then
+    If (TWetBulb >= FREEZING_POINT_WATER_SI) Then
       HumRatio = ((2501 - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) / (2501 + 1.86 * TDryBulb - 4.186 * TWetBulb)
     Else
       HumRatio = ((2830 - 0.24 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) / (2830 + 1.86 * TDryBulb - 2.1 * TWetBulb)
@@ -1226,6 +1214,12 @@ Function GetSatVapPres(ByVal TDryBulb As Variant) As Variant
 '
 ' Reference:
 '        ASHRAE Handbook - Fundamentals (2017) ch. 1  eqn 5 & 6
+'        Important note: the ASHRAE formulae are defined above and below the freezing point but have
+'        a discontinuity at the freezing point. This is a small inaccuracy on ASHRAE's part: the formulae
+'        should be defined above and below the triple point of water (not the feezing point) in which case
+'        the discontinuity vanishes. It is essential to use the triple point of water otherwise function
+'        GetTDewPointFromVapPres, which inverts the present function, does not converge properly around
+'        the freezing point.
 '
   Dim LnPws As Variant, T As Variant
 
@@ -1239,7 +1233,7 @@ Function GetSatVapPres(ByVal TDryBulb As Variant) As Variant
 
     T = GetTRankineFromTFahrenheit(TDryBulb)
 
-    If (TDryBulb <= 32) Then
+    If (TDryBulb <= TRIPLE_POINT_WATER_IP) Then
       LnPws = (-10214.165 / T - 4.8932428 - 0.0053765794 * T + 0.00000019202377 * T ^ 2 _
             + 3.5575832E-10 * T ^ 3 - 9.0344688E-14 * T ^ 4 + 4.1635019 * Log(T))
     Else
@@ -1255,7 +1249,7 @@ Function GetSatVapPres(ByVal TDryBulb As Variant) As Variant
 
     T = GetTKelvinFromTCelsius(TDryBulb)
 
-    If (TDryBulb <= 0) Then
+    If (TDryBulb <= TRIPLE_POINT_WATER_SI) 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)
     Else
diff --git a/tests/vba/test_psychrolib_ip.xlsm b/tests/vba/test_psychrolib_ip.xlsm
index 8b604d28..2dcb24f4 100755
Binary files a/tests/vba/test_psychrolib_ip.xlsm and b/tests/vba/test_psychrolib_ip.xlsm differ
diff --git a/tests/vba/test_psychrolib_si.bas b/tests/vba/test_psychrolib_si.bas
index 660a064b..28054cfb 100644
--- a/tests/vba/test_psychrolib_si.bas
+++ b/tests/vba/test_psychrolib_si.bas
@@ -61,7 +61,7 @@ End Sub
 
 ' Test of helper functions
 Sub test_GetTKelvinFromTCelsius()
-  Call TestExpression("", GetTKelvinFromTCelsius(20), 293.15, 0.000001)
+  Call TestExpression("GetTKelvinFromTCelsius", GetTKelvinFromTCelsius(20), 293.15, 0.000001)
 End Sub
 
 
@@ -287,4 +287,3 @@ Sub test_GetTDewPointFromVapPres_convergence()
   Call TestExpression("GetTDewPointFromVapPres convergence test", 1, 1, abst:=0.1)
 End Sub
 
-
diff --git a/tests/vba/test_psychrolib_si.xlsm b/tests/vba/test_psychrolib_si.xlsm
index cfbe0674..e279747c 100755
Binary files a/tests/vba/test_psychrolib_si.xlsm and b/tests/vba/test_psychrolib_si.xlsm differ