diff --git a/project.clj b/project.clj new file mode 100644 index 0000000..326d174 --- /dev/null +++ b/project.clj @@ -0,0 +1,7 @@ +(defproject humidity "0.1.0-SNAPSHOT" + :description "FIXME: write description" + :url "http://example.com/FIXME" + :license {:name "EPL-2.0 OR GPL-2.0-or-later WITH Classpath-exception-2.0" + :url "https://www.eclipse.org/legal/epl-2.0/"} + :dependencies [[org.clojure/clojure "1.11.1"]] + :repl-options {:init-ns humidity.core}) diff --git a/resources/absolute-humidity-eqn.gif b/resources/absolute-humidity-eqn.gif new file mode 100644 index 0000000..e620841 Binary files /dev/null and b/resources/absolute-humidity-eqn.gif differ diff --git a/resources/katzenberger.txt b/resources/katzenberger.txt new file mode 100644 index 0000000..76c60dc --- /dev/null +++ b/resources/katzenberger.txt @@ -0,0 +1,49 @@ +# I'm guessing this is Python + +# Saturated water vapor pressure in hPa, Magnus Formula +# Parameters from Sonntag1990, for −45 °C ≤ T ≤ 60 °C (error ±0.35 °C). +def saturatedVaporPressureMagnusSonntag1990(temperatureCelsius): + a = 6.112 + b = 17.62 + c = 243.12 + + saturatedVaporPressure = a * math.exp((b * temperatureCelsius) / (c + temperatureCelsius)) + return saturatedVaporPressure + + +# Saturated water vapor pressure in hPa, Tetens Formula +def saturatedVaporPressureTetens(temperatureCelsius): + a = 6.1078 + b = 17.27 + c = 237.3 + + saturatedVaporPressure = a * math.exp((b * temperatureCelsius) / (c + temperatureCelsius)) + return saturatedVaporPressure + + +# Saturated water vapor pressure in hPa, Buck 1996 Formula +def saturatedVaporPressureBuck1996(temperatureCelsius): + if temperatureCelsius < 0: + return 6.1115 * math.exp( + (23.06 - temperatureCelsius / 333.7) * (temperatureCelsius / (279.82 + temperatureCelsius))) + else: + return 6.1121 * math.exp( + (18.678 - temperatureCelsius / 234.5) * (temperatureCelsius / (257.14 + temperatureCelsius))) + + +# Actual water vapor pressure in hPa +def vaporPressure(relativeHumidity, temperatureCelsius): + vaporPressure = relativeHumidity / 100.0 * saturatedVaporPressureBuck1996(temperatureCelsius) + return vaporPressure + + +# Absolute humidity in g/m³ +def absoluteHumidity(relativeHumidity, temperatureCelsius): + molarMassOfWaterVapor = 18.01528 + universalGasConstant = 8314.46261815324 + zeroCelsiusInKelvin = 273.15 + + absoluteHumidity = 10 ** 5 * molarMassOfWaterVapor / universalGasConstant * vaporPressure(relativeHumidity, + temperatureCelsius) / ( + temperatureCelsius + zeroCelsiusInKelvin) + return absoluteHumidity \ No newline at end of file diff --git a/resources/saturation-vapour-pressure-eqn.png b/resources/saturation-vapour-pressure-eqn.png new file mode 100644 index 0000000..d2dd39d Binary files /dev/null and b/resources/saturation-vapour-pressure-eqn.png differ diff --git a/resources/wet-bulb-temp-eqn.png b/resources/wet-bulb-temp-eqn.png new file mode 100644 index 0000000..df55be6 Binary files /dev/null and b/resources/wet-bulb-temp-eqn.png differ diff --git a/src/humidity/absolute.clj b/src/humidity/absolute.clj new file mode 100644 index 0000000..4bd1367 --- /dev/null +++ b/src/humidity/absolute.clj @@ -0,0 +1,38 @@ +(ns ^{:doc "Function to calculate wet bulb temperature (in ° C) given temperature in + ° C and relative-humidity percentage (0...100). + + See [authority](https://carnotcycle.wordpress.com/2012/08/04/how-to-convert-relative-humidity-to-absolute-humidity/)"} + humidity.absolute + (:require + [humidity.constants :refer [Rw]])) + +;; Absolute Humidity (grams/m3) = 6.112 × e^[(17.67 × T)/(T+243.5)] × rh × 2.1674 +;; ----- +;; (273.15+T) + +(def ^:const magic-numbers + "The magic number constants from the derived expression. It would give me + more comfort if there were some agreement between these magic numbers and + those used in the wet bulb function, but there is not." + [0 6.112 17.67 243.5 2.1674 273.15]) + +(def C magic-numbers) + +;; (defn abs-humidity +;; "Returns absolute humidity in grammes per cubic metre given this +;; `temperature` in ° C and `relative-hunidity` percentage (0...100). +;; ![function for absolute humidity given temperature and relative humidity](https://carnotcycle.wordpress.com/wp-content/uploads/2014/09/ah3.gif)" +;; [temperature relative-humidity] +;; (when (empty? (filter false? (map number? [temperature relative-humidity]))) +;; (let [T temperature +;; RHP relative-humidity] +;; (/ (* (C 1) (pow (* (/ (C 2) T) (+ T (C 3))) e) RHP (C 4)) +;; (+ celsius-offset temperature))))) + +;; AH=(RH × Ps)​​/(Rw ​× T × 100) + +(defn abs-humidity [& {:keys [relative-humidity saturation-pressure temperature-kelvin]}] + (when (and relative-humidity saturation-pressure temperature-kelvin) + (/ (* relative-humidity saturation-pressure) (* Rw temperature-kelvin 100))) + ) + diff --git a/src/humidity/constants.clj b/src/humidity/constants.clj new file mode 100644 index 0000000..94bcead --- /dev/null +++ b/src/humidity/constants.clj @@ -0,0 +1,22 @@ +(ns humidity.constants) + + +(def ^:const celsius-offset + "Value of 0° Celsius in ° Kelvin" + 273.15) + +(def ^:const e + "Base of natural logarithms" + 2.71828) + +(def ^:const Pc + "Critical pressure for water, in MPa (Megapascal)" + 22.064) + +(def ^:const Rw + "Specific gas constant for water vapour (J/(kg⋅K))" + 461.5) + +(def ^:const Tc + "Critical temperature for water, ° Kelvin" + 647.096) \ No newline at end of file diff --git a/src/humidity/core.clj b/src/humidity/core.clj new file mode 100644 index 0000000..893bd93 --- /dev/null +++ b/src/humidity/core.clj @@ -0,0 +1,49 @@ +(ns humidity.core + (:require [humidity.absolute :refer [abs-humidity]] + [humidity.constants :refer [celsius-offset]] + [humidity.saturation-vp :as saturation-vp] + [humidity.utils :refer [?assoc]] + [humidity.wet-bulb :refer [wet-bulb-temperature]])) + +(defn resolve-humidity + "Resolve humidity equations given the data in this packet." + [& {:keys [temperature temperature-celsius temperature-kelvin relative-humidity + pressure-kilopascals pressure-millibars] :as packet}] + (println packet) + (if temperature (resolve-humidity (assoc (dissoc packet :temperature) :temperature-celsius temperature)) + (?assoc packet + :temperature-celsius (when + (and temperature-kelvin (not temperature-celsius)) + (- temperature-kelvin celsius-offset)) + :temperature-kelvin (when + (and temperature-celsius (not temperature-kelvin)) + (+ temperature-celsius celsius-offset)) + :wet-bulb-temperature (when + (and temperature-celsius relative-humidity) + (wet-bulb-temperature temperature-celsius relative-humidity)) + :absolute-humidity (when (and temperature-celsius relative-humidity) + (abs-humidity temperature-celsius relative-humidity)) + :pressure-kilopascals (when-not pressure-kilopascals + (if pressure-millibars (/ pressure-millibars 10) 100)) + :pressure-millibars (when-not pressure-millibars + (if pressure-kilopascals (* pressure-kilopascals 10) 1000)) + :saturation-pressure (saturation-vp/saturation-vp packet)))) + + +;; the example in the paper gives this as 13.7 +(wet-bulb-temperature 20 50) +;; by definition, this must equal 35 +(wet-bulb-temperature 35 100) +;; from published graphs, I interpolatethis in the range 33...35 +(wet-bulb-temperature 40 70) + +;; (resolve-humidity {:relative-humidity 100 :temperature 35}) + +(defn resolution + "Repeatedly run `resolve-humidity`, q.v., until we can infer no further + information." + [packet] + (loop [p packet] + (let [p' (resolve-humidity p)] + (if (= (keys p) (keys p')) p + (recur p'))))) \ No newline at end of file diff --git a/src/humidity/relative.clj b/src/humidity/relative.clj new file mode 100644 index 0000000..b45c850 --- /dev/null +++ b/src/humidity/relative.clj @@ -0,0 +1 @@ +(ns humidity.relative) \ No newline at end of file diff --git a/src/humidity/saturation_vp.clj b/src/humidity/saturation_vp.clj new file mode 100644 index 0000000..b7dd8d4 --- /dev/null +++ b/src/humidity/saturation_vp.clj @@ -0,0 +1,25 @@ +(ns ^{:doc "Saturation vapour pressure function" + :authority "https://en.wikipedia.org/wiki/Tetens_equation" + :author "simon"} + humidity.saturation-vp + (:require [clojure.math :refer [pow]] + [humidity.constants :refer [celsius-offset e]])) + + +(defn saturation-vp-tetens + "Calculate and return saturation vapour pressure in hPa using Tetens' + formula, given temerature in ° Celsius." + [temperature-celsius] + (let [a 6.1078 ;; yet more magic numbers... + b 17.27 + c 237.3] + (* a (pow e (/ (* b temperature-celsius) (+ c temperature-celsius)))))) + +(defn saturation-vp + "Calculate and return saturation vapour pressure, is possible, + given the data in this packet." + [packet] + (let [Tc (packet :temperature-celsius) + Tk (packet :temperature-kelvin)] + (cond Tc (saturation-vp-tetens Tc) + Tk (saturation-vp-tetens (- Tk celsius-offset))))) diff --git a/src/humidity/utils.clj b/src/humidity/utils.clj new file mode 100644 index 0000000..1fef5ca --- /dev/null +++ b/src/humidity/utils.clj @@ -0,0 +1,27 @@ +(ns ^{:doc "Utility functions for humidity calculation"} humidity.utils + (:require [clojure.math :refer [pow]])) + +(def ^:const celsius-offset + "Value of 0° Celsius in ° Kelvin" + 273.15) + +(def ^:const e + "Base of natural logarithms" + 2.71828) + + +(defmacro expt + "Clojure.math `pow` function does not allow ratios as exponents; + this is equivalent but coerces the exponent argument to float." + [base exponent] + `(pow ~base (float ~exponent))) + +(defn ?assoc + "Same as assoc, but skip the assoc if v is nil. + Borrowed from https://stackoverflow.com/questions/16356888/assoc-if-in-clojure" + [m & kvs] + (->> kvs + (partition 2) + (filter second) + (map vec) + (into m))) \ No newline at end of file diff --git a/src/humidity/wet_bulb.clj b/src/humidity/wet_bulb.clj new file mode 100644 index 0000000..5443d0a --- /dev/null +++ b/src/humidity/wet_bulb.clj @@ -0,0 +1,26 @@ +(ns ^{:doc "Function to calculate wet bulb temperature (in ° C) given temperature in + ° C and relative-humidity percentage (0...100)"} + humidity.wet-bulb + (:require [clojure.math :refer [atan]] + [humidity.utils :refer [expt]])) + +(def ^:const magic-numbers + "The magic numbers (empirical constants) from the derived expression. It + bothers me that all these formulae depend on 'empirical constants', *and* + that there is no commonality between the empirical constants of different + formulae." + [0 0.151977 8.313659 1.676331 0.00391838 0.023101 4.686035]) + +(def C magic-numbers) + +(defn wet-bulb-temperature + "Returns wet bulb temperature (in ° C) given this `temperature` in + ° C and `relative-humidity` percentage (0...100)" + [temperature relative-humidity] + (let [T temperature + RHP relative-humidity] + (+ (* T (atan (* (C 1) (expt (+ RHP (C 2)) 1/2)))) + (atan (+ T RHP)) + (- 0 (atan (- RHP (C 3)))) + (* (C 4) (expt RHP 3/2) (atan (* (C 5) RHP))) + (- 0 (C 6))))) diff --git a/test/humidity/wet_bulb_test.clj b/test/humidity/wet_bulb_test.clj new file mode 100644 index 0000000..6dc802b --- /dev/null +++ b/test/humidity/wet_bulb_test.clj @@ -0,0 +1,47 @@ +(ns humidity.wet-bulb-test + (:require [clojure.test :refer [deftest is testing]] + [humidity.wet-bulb :refer [wet-bulb-temperature]])) + + +(deftest wet-bulb-test + (testing "Example from paper" + (let [expected 13.7 + temperature 20 + rel-humidity 50 + actual (wet-bulb-temperature temperature rel-humidity)] + (println (format "Expected: %s; actual: %s" expected actual)) + (is (< (abs (- expected actual)) 0.001)))) + + (testing "Interpolations from graph" + (let [min-expected 31 + max-expected 33 + temperature 37 + rel-humidity 70 + actual (wet-bulb-temperature temperature rel-humidity)] + (println (format "Expected: %s...%s; actual: %s" min-expected max-expected actual)) + (is (< min-expected actual max-expected))) + + (let [min-expected 34 + max-expected 36 + temperature 39.7 + rel-humidity 70 + actual (wet-bulb-temperature temperature rel-humidity)] + (println (format "Expected: %s...%s; actual: %s" min-expected max-expected actual)) + (is (< min-expected actual max-expected))) + + + (let [min-expected 31 + max-expected 34 + temperature 40 + rel-humidity 60 + actual (wet-bulb-temperature temperature rel-humidity)] + (println (format "Expected: %s...%s; actual: %s" min-expected max-expected actual)) + (is (< min-expected actual max-expected))) + + (let [min-expected 32 + max-expected 35 + temperature 51 + rel-humidity 30 + actual (wet-bulb-temperature temperature rel-humidity)] + (println (format "Expected: %s...%s; actual: %s" min-expected max-expected actual)) + (is (< min-expected actual max-expected))))) \ No newline at end of file