From 08664fa668fb6e13c8e1d7222be239d5fa7ff814 Mon Sep 17 00:00:00 2001 From: Patrick Grawehr Date: Fri, 7 Oct 2022 08:57:19 +0200 Subject: [PATCH] Add a few new WeatherHelper functions (#1941) --- .../Common/Iot/Device/Common/WeatherHelper.cs | 104 +++++++++++++++++- src/devices/Common/tests/WeatherTests.cs | 50 +++++++++ 2 files changed, 149 insertions(+), 5 deletions(-) diff --git a/src/devices/Common/Iot/Device/Common/WeatherHelper.cs b/src/devices/Common/Iot/Device/Common/WeatherHelper.cs index af26066c..571b3553 100644 --- a/src/devices/Common/Iot/Device/Common/WeatherHelper.cs +++ b/src/devices/Common/Iot/Device/Common/WeatherHelper.cs @@ -12,6 +12,21 @@ namespace Iot.Device.Common /// public static class WeatherHelper { + /// + /// Gas constant of dry Air, J / (kg * K) + /// + internal const double SpecificGasConstantOfAir = 287.058; + + /// + /// Gas constant of vapor, J / (kg * K) + /// + internal const double SpecificGasConstantOfVapor = 461.523; + + /// + /// Default atmospheric temperature gradient = 0.0065K/m (or 0.65K per 100m) + /// + internal const double DefaultTemperatureGradient = 0.0065; + /// /// The mean sea-level pressure (MSLP) is the average atmospheric pressure at mean sea level /// @@ -183,7 +198,7 @@ namespace Iot.Device.Common /// The altitude public static Length CalculateAltitude(Pressure pressure, Pressure seaLevelPressure, Temperature airTemperature) { - double meters = ((Math.Pow(seaLevelPressure.Pascals / pressure.Pascals, 1 / 5.255) - 1) * airTemperature.Kelvins) / 0.0065; + double meters = ((Math.Pow(seaLevelPressure.Pascals / pressure.Pascals, 1 / 5.255) - 1) * airTemperature.Kelvins) / DefaultTemperatureGradient; return Length.FromMeters(meters); } @@ -222,7 +237,7 @@ namespace Iot.Device.Common /// The estimated absolute sea-level pressure /// solved for sea level pressure public static Pressure CalculateSeaLevelPressure(Pressure pressure, Length altitude, Temperature airTemperature) => - Pressure.FromPascals(Math.Pow((((0.0065 * altitude.Meters) / airTemperature.Kelvins) + 1), 5.255) * pressure.Pascals); + Pressure.FromPascals(Math.Pow((((DefaultTemperatureGradient * altitude.Meters) / airTemperature.Kelvins) + 1), 5.255) * pressure.Pascals); /// /// Calculates the approximate absolute pressure from given sea-level pressure, altitude and air temperature. @@ -232,7 +247,7 @@ namespace Iot.Device.Common /// The air temperature at the point for which pressure is being calculated /// The estimated absolute pressure at the given altitude public static Pressure CalculatePressure(Pressure seaLevelPressure, Length altitude, Temperature airTemperature) => - Pressure.FromPascals(seaLevelPressure.Pascals / Math.Pow((((0.0065 * altitude.Meters) / airTemperature.Kelvins) + 1), 5.255)); + Pressure.FromPascals(seaLevelPressure.Pascals / Math.Pow((((DefaultTemperatureGradient * altitude.Meters) / airTemperature.Kelvins) + 1), 5.255)); /// /// Calculates the temperature gradient for the given pressure difference @@ -243,7 +258,7 @@ namespace Iot.Device.Common /// The standard temperature at the given altitude, when the given pressure difference is known /// solved for temperature public static Temperature CalculateTemperature(Pressure pressure, Pressure seaLevelPressure, Length altitude) => - Temperature.FromKelvins((0.0065 * altitude.Meters) / (Math.Pow(seaLevelPressure.Pascals / pressure.Pascals, 1 / 5.255) - 1)); + Temperature.FromKelvins((DefaultTemperatureGradient * altitude.Meters) / (Math.Pow(seaLevelPressure.Pascals / pressure.Pascals, 1 / 5.255) - 1)); /// /// Calculates the barometric pressure from a raw reading, using the reduction formula from the german met service. @@ -302,7 +317,7 @@ namespace Iot.Device.Common Length measurementAltitude) { double x = (9.80665 / (287.05 * ((measuredTemperature.Kelvins) + 0.12 * vaporPressure.Hectopascals + - (0.0065 * measurementAltitude.Meters) / 2))) * measurementAltitude.Meters; + (DefaultTemperatureGradient * measurementAltitude.Meters) / 2))) * measurementAltitude.Meters; double barometricPressure = measuredPressure.Hectopascals * Math.Exp(x); return Pressure.FromHectopascals(barometricPressure); } @@ -331,5 +346,84 @@ namespace Iot.Device.Common } #endregion + + /// + /// Simplified air density (not taking humidity into account) + /// + /// Measured air pressure + /// Measured temperature + /// Approximate standard air density + /// From https://de.wikipedia.org/wiki/Luftdichte + public static Density CalculateAirDensity(Pressure airPressure, Temperature temperature) + { + var result = airPressure.Pascals / (SpecificGasConstantOfAir * temperature.Kelvins); + return Density.FromKilogramsPerCubicMeter(result); + } + + /// + /// Calculates the air density + /// + /// Measured air pressure + /// Measured temperature + /// Measured relative humidity + /// Approximate standard air density at sea level + /// From https://de.wikipedia.org/wiki/Luftdichte + public static Density CalculateAirDensity(Pressure airPressure, Temperature temperature, RelativeHumidity humidity) + { + double rs = SpecificGasConstantOfAir; + double rd = SpecificGasConstantOfVapor; + var pd = CalculateSaturatedVaporPressureOverWater(temperature); + // It's still called "constant" even though it's not constant + double gasConstant = rs / (1 - (humidity.Percent / 100) * (pd.Pascals / airPressure.Pascals) * (1 - (rs / rd))); + var result = airPressure.Pascals / (gasConstant * temperature.Kelvins); + return Density.FromKilogramsPerCubicMeter(result); + } + + /// + /// Calculates the wind chill temperature - this is the perceived temperature in (heavy) winds at cold temperatures. + /// This is only useful at temperatures below about 20°C, above use instead. + /// Not suitable for wind speeds < 5 km/h. + /// + /// The measured air temperature + /// The wind speed (measured at 10m above ground) + /// The perceived temperature. Note that this is not a real temperature, and the skin will never really reach + /// this temperature. This is more an indication on how fast the skin will reach the air temperature. If the skin + /// reaches a temperature of about -5°C, frostbite might occur. + /// From https://de.wikipedia.org/wiki/Windchill + /// The wind speed is less than zero + public static Temperature CalculateWindchill(Temperature temperature, Speed windSpeed) + { + if (windSpeed < Speed.Zero) + { + throw new ArgumentOutOfRangeException(nameof(windSpeed), "The wind speed cannot be negative"); + } + + double va = temperature.DegreesCelsius; + if (windSpeed < Speed.FromKilometersPerHour(1)) + { + // otherwise, the result is unusable, because the second and third terms of the equation are 0, resulting in a constant offset from the input temperature + windSpeed = Speed.FromKilometersPerHour(1); + } + + double wct = 13.12 + 0.6215 * va + (0.3965 * va - 11.37) * Math.Pow(windSpeed.KilometersPerHour, 0.16); + return Temperature.FromDegreesCelsius(wct); + } + + /// + /// Calculates the wind force on an object. + /// + /// The denisty of the air, calculated using one of the overloads of + /// The speed of the wind + /// Pressure coefficient for the shape of the object. Use 1 for a rectangular object directly facing the wind + /// The Pressure the wind applies on the object + /// From https://de.wikipedia.org/wiki/Winddruck + public static Pressure CalculateWindForce(Density densityOfAir, Speed windSpeed, double pressureCoefficient = 1.0) + { + double v = windSpeed.MetersPerSecond; + double rho = densityOfAir.KilogramsPerCubicMeter; + + double wd = pressureCoefficient * rho / 2 * (v * v); + return Pressure.FromNewtonsPerSquareMeter(wd); + } } } diff --git a/src/devices/Common/tests/WeatherTests.cs b/src/devices/Common/tests/WeatherTests.cs index 59cc2466..6c2b06ca 100644 --- a/src/devices/Common/tests/WeatherTests.cs +++ b/src/devices/Common/tests/WeatherTests.cs @@ -185,5 +185,55 @@ namespace Iot.Device.Common.Tests Assert.Equal(outHumidityExpected, result.Percent, 3); } + + [Theory] + [InlineData(1013.25, 35, 1.1455)] + [InlineData(1013.25, 0, 1.2922)] + [InlineData(1013.25, -25, 1.4224)] + public void CalculateAirDensitySimple(double inPressure, double inTemp, double expected) + { + var result = WeatherHelper.CalculateAirDensity(Pressure.FromMillibars(inPressure), + Temperature.FromDegreesCelsius(inTemp)); + + Assert.Equal(expected, result.KilogramsPerCubicMeter, 4); + } + + [Theory] + [InlineData(1013.25, 35, 0, 1.1455)] + [InlineData(1013.25, 35, 50, 1.1335)] + [InlineData(1013.25, 35, 100, 1.1214)] + public void CalculateAirDensity(double inPressure, double inTemp, double inHumidity, double expected) + { + var result = WeatherHelper.CalculateAirDensity(Pressure.FromMillibars(inPressure), + Temperature.FromDegreesCelsius(inTemp), RelativeHumidity.FromPercent(inHumidity)); + + Assert.Equal(expected, result.KilogramsPerCubicMeter, 4); + } + + [Theory] + // When the wind is calm, the windchill temperature is higher than the actual temperature, because the air close to the skin heats up. + // Explained in https://de.wikipedia.org/wiki/Windchill + [InlineData(0, 0, 1.75)] + [InlineData(10, 5, 9.8)] + [InlineData(-10, 20, -17.9)] + [InlineData(-15, 40, -27.4)] + public void CalculateWindchill(double temperature, double windSpeed, double expected) + { + var result = WeatherHelper.CalculateWindchill(Temperature.FromDegreesCelsius(temperature), + Speed.FromKilometersPerHour(windSpeed)); + + Assert.Equal(expected, result.DegreesCelsius, 1); + } + + [Theory] + [InlineData(20, 38.5, 1013, 68.83)] + [InlineData(20, 88.9, 1013, 367.0)] + public void CalculateWindForce(double temperature, double windSpeed, double pressure, double expected) + { + Density airDensity = WeatherHelper.CalculateAirDensity(Pressure.FromHectopascals(pressure), + Temperature.FromDegreesCelsius(temperature)); + var density = WeatherHelper.CalculateWindForce(airDensity, Speed.FromKilometersPerHour(windSpeed), 1.0); + Assert.Equal(expected, density.NewtonsPerSquareMeter, 1); + } } } -- GitLab