diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 06d4d6098a..b32aeb41be 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,46 +1,12 @@ #include "virtual_temperature.hpp" #include +#include #include namespace py = pybind11; -// Flexible dispatcher that supports scalar/list combinations -std::vector virtual_temperature_py(py::object temperature, py::object mixing_ratio, double epsilon) { - std::vector temperature_vec; - std::vector mixing_ratio_vec; - - // Case 1: both are lists - if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - temperature_vec = temperature.cast>(); - mixing_ratio_vec = mixing_ratio.cast>(); - if (temperature_vec.size() != mixing_ratio_vec.size()) { - throw std::invalid_argument("Temperature and mixing ratio lists must be the same length."); - } - } - // Case 2: temperature is float, mixing_ratio is list - else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - mixing_ratio_vec = mixing_ratio.cast>(); - temperature_vec = std::vector(mixing_ratio_vec.size(), temperature.cast()); - } - // Case 3: temperature is list, mixing_ratio is float - else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - temperature_vec = temperature.cast>(); - mixing_ratio_vec = std::vector(temperature_vec.size(), mixing_ratio.cast()); - } - // Case 4: both are floats - else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - temperature_vec = {temperature.cast()}; - mixing_ratio_vec = {mixing_ratio.cast()}; - } - else { - throw std::invalid_argument("Inputs must be float or list."); - } - - return VirtualTemperature(temperature_vec, mixing_ratio_vec, epsilon); -} - int add(int i, int j) { - return i - j; + return i + j; } PYBIND11_MODULE(_calc_mod, m) { @@ -49,11 +15,11 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("add", &add, "Add two numbers"); // Unified binding with default epsilon - m.def("virtual_temperature", &virtual_temperature_py, - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622, - "Compute virtual temperature.\n" - "Accepts:\n" - " - two lists of equal length\n" - " - one scalar and one list\n" - "Defaults to epsilon = 0.622"); + m.def("dewpoint", py::vectorize(DewPoint), + "Calculate dew point from water vapor partial pressure.", + py::arg("vapor_pressure")); + + m.def("virtual_temperature", py::vectorize(VirtualTemperature), + "Calculate virtual temperature from temperature and mixing ratio.", + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622); } diff --git a/src/constants.hpp b/src/constants.hpp new file mode 100644 index 0000000000..8e9d596a61 --- /dev/null +++ b/src/constants.hpp @@ -0,0 +1,37 @@ +#ifndef CONSTANTS_HPP +#define CONSTANTS_HPP + +namespace metpy_constants { +// Gas constants (J / kg / K) +constexpr double R = 8.314462618; // Universal gas constant (J / mol / K) + +// Dry air +constexpr double Md = 28.96546e-3; // Molar mass of dry air (kg / mol) +constexpr double Rd = R / Md; // Dry air (J / kg / K) +constexpr double dry_air_spec_heat_ratio = 1.4; +constexpr double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) + +// Water +constexpr double Mw = 18.015268e-3; // Molar mass of water (kg / mol) +constexpr double Rv = R / Mw; // Water vapor (J / kg / K) +constexpr double wv_spec_heat_ratio = 1.33; +constexpr double Cp_v = Rv * wv_spec_heat_ratio / (wv_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) +constexpr double Cp_l = 4.2194e3; // Specific heat capacity of liquid water (J / kg / K) +constexpr double Lv = 2.50084e6; // Latent heat of vaporization of water (J / kg) +constexpr double T0 = 273.16; // Triple point of water (K) + +// General Meteorological constants +constexpr double epsilon = Mw / Md; // ≈ 0.622 + + +// Standard gravity +constexpr double g = 9.80665; // m / s^2 + +// Reference pressure +constexpr double P0 = 100000.0; // Pa (hPa = 1000) + +} + +#endif // CONSTANTS_HPP diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 592365cce0..0204a07066 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,19 +1,29 @@ +#include +#include "constants.hpp" #include "virtual_temperature.hpp" //#include -std::vector VirtualTemperature( - const std::vector& temperature, - const std::vector& mixing_ratio, - double epsilon -) { +double water_latent_heat_vaporization(double temperature) { + // Calculate the latent heat of vaporization of water in J/kg at a given temperature. + using namespace metpy_constants; + return Lv - (Cp_l - Cp_v) * (temperature - T0); +} + +double _saturation_vapor_pressure(double temperature) { + // Calculate saturation (equilibrium) water vapor (partial) pressure over liquid water. + // Constants for the Magnus-Tetens approximation + //const double a = 17.67; + //const double b = 243.5; - std::vector result(temperature.size()); - double T, w; - for (size_t i = 0; i < temperature.size(); ++i) { - T = temperature[i]; - w = mixing_ratio[i]; - result[i] = T * (w + epsilon) / (epsilon * (1 + w)); - } + // Calculate saturation vapor pressure using the Magnus-Tetens formula + //return 6.112 * exp((a * temperature) / (b + temperature)); +} + +double DewPoint(double vapor_pressure) { + double val = log(vapor_pressure / 6.112); + return 243.5 * val / (17.67 - val); +} - return result; +double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { + return temperature * (mixing_ratio + epsilon) / (epsilon * (1. + mixing_ratio)); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 2fc91564ef..4b871bebd0 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -1,12 +1,10 @@ #ifndef VIRTUAL_TEMPERATURE_HPP // if not defined #define VIRTUAL_TEMPERATURE_HPP // define the header file -#include +//#include -// Compute virtual temperature: vector + vector -std::vector VirtualTemperature( - const std::vector& temperature, - const std::vector& mixing_ratio, - double epsilon); +double DewPoint(double vapor_pressure); + +double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); #endif // VIRTUAL_TEMPERATURE_HPP