Pybind11_vectorize (#1)

* vectorize the function using pybind11/numpy.h, achieving handling numpy array, float and dtype

* add dewpoint, add constants.hpp, adding other functions WIP

Author: gitlinffff

src/calcmod.cpp

@@ -1,46 +1,12 @@
 #include "virtual_temperature.hpp"
 #include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
 #include <pybind11/stl.h>
 
 namespace py = pybind11;
 
-// Flexible dispatcher that supports scalar/list combinations
-std::vector<double> virtual_temperature_py(py::object temperature, py::object mixing_ratio, double epsilon) {
-    std::vector<double> temperature_vec;
-    std::vector<double> mixing_ratio_vec;
-
-    // Case 1: both are lists
-    if (py::isinstance<py::list>(temperature) && py::isinstance<py::list>(mixing_ratio)) {
-        temperature_vec = temperature.cast<std::vector<double>>();
-        mixing_ratio_vec = mixing_ratio.cast<std::vector<double>>();
-        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<py::float_>(temperature) && py::isinstance<py::list>(mixing_ratio)) {
-        mixing_ratio_vec = mixing_ratio.cast<std::vector<double>>();
-        temperature_vec = std::vector<double>(mixing_ratio_vec.size(), temperature.cast<double>());
-    }
-    // Case 3: temperature is list, mixing_ratio is float
-    else if (py::isinstance<py::list>(temperature) && py::isinstance<py::float_>(mixing_ratio)) {
-        temperature_vec = temperature.cast<std::vector<double>>();
-        mixing_ratio_vec = std::vector<double>(temperature_vec.size(), mixing_ratio.cast<double>());
-    }
-    // Case 4: both are floats
-    else if (py::isinstance<py::float_>(temperature) && py::isinstance<py::float_>(mixing_ratio)) {
-        temperature_vec = {temperature.cast<double>()};
-        mixing_ratio_vec = {mixing_ratio.cast<double>()};
-    }
-    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);
 }

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

src/virtual_temperature.cpp

@@ -1,19 +1,29 @@
+#include <cmath>
+#include "constants.hpp"
 #include "virtual_temperature.hpp"
 //#include <stdexcept>
 
-std::vector<double> VirtualTemperature(
-    const std::vector<double>& temperature,
-    const std::vector<double>& 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<double> 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));
 }

src/virtual_temperature.hpp

@@ -1,12 +1,10 @@
 #ifndef VIRTUAL_TEMPERATURE_HPP // if not defined
 #define VIRTUAL_TEMPERATURE_HPP // define the header file
 
-#include <vector>
+//#include <vector>
 
-// Compute virtual temperature: vector + vector
-std::vector<double> VirtualTemperature(
-    const std::vector<double>& temperature,
-    const std::vector<double>& mixing_ratio,
-    double epsilon);
+double DewPoint(double vapor_pressure);
+
+double VirtualTemperature(double temperature, double mixing_ratio, double epsilon);
 
 #endif // VIRTUAL_TEMPERATURE_HPP