💨wind vector stuff

src/include/libthermo.hpp

@@ -66,12 +66,12 @@ constexpr T potential_temperature(const T pressure, const T temperature) noexcep
 template <floating T>
 constexpr T equivalent_potential_temperature(
   const T pressure, const T temperature, const T dewpoint
-) noexcept;  // theta_e
+) noexcept;
 
 template <floating T>
 constexpr T wet_bulb_potential_temperature(
   const T pressure, const T temperature, const T dewpoint
-) noexcept;  // theta_w
+) noexcept;
 
 /* ........................................{ ode    }........................................... */
 
@@ -98,16 +98,7 @@ class lcl {
     // constructor with the function arguments.
   public:
     T pressure, temperature;
-    // default constructor
     constexpr lcl() noexcept = default;
-    // copy constructor
-    constexpr lcl(const lcl<T>& other) noexcept = default;
-    // move constructor
-    constexpr lcl(lcl<T>&& other) noexcept = default;
-    // copy assignment operator
-    constexpr lcl<T>& operator=(const lcl<T>& other) noexcept = default;
-    // destructor
-    ~lcl() noexcept = default;
     constexpr lcl(const T pressure, const T temperature) noexcept :
         pressure(pressure), temperature(temperature){};
 

src/include/wind.hpp

@@ -1,35 +1,39 @@
 #pragma once
 
-#include <algorithm>
 #include <cmath>
-#include <functional>
-#include <vector>
-#include <array>
 
-#include <functional.hpp>
 #include <common.hpp>
 
 namespace libthermo {
 
 template <floating T>
-struct WindVector {
-    T speed;
-    T direction;
-};
+class wind_vector;
 
 template <floating T>
-struct WindComponents {
-    T u;
-    T v;
-};
+class wind_components;
 
 template <floating T>
-constexpr WindComponents<T> wind_components(const T direction, const T magnitude) noexcept;
+constexpr T wind_direction(const T u, const T v) noexcept;
 
 template <floating T>
-constexpr T wind_direction(const T u, const T v, const bool from = false) noexcept;
+constexpr T wind_magnitude(const T u, const T v) noexcept;
 
 template <floating T>
-constexpr T wind_magnitude(const T u, const T v) noexcept;
+class wind_components {
+  public:
+    T u, v;
+    constexpr wind_components() noexcept = default;
+    constexpr wind_components(const T u, const T v) noexcept : u(u), v(v) {}
+    constexpr explicit wind_components(const wind_vector<T>& uv) noexcept;
+};
+
+template <floating T>
+class wind_vector {
+  public:
+    T direction, magnitude;
+    constexpr wind_vector() noexcept = default;
+    constexpr wind_vector(const T d, const T m) noexcept : direction(d), magnitude(m) {}
+    constexpr explicit wind_vector(const wind_components<T>& uv) noexcept;
+};
 
 }  // namespace libthermo

src/lib/libthermo.cpp

@@ -58,7 +58,15 @@ template <floating T>
 constexpr T potential_temperature(const T pressure, const T temperature) noexcept {
     return temperature / exner_function(pressure);
 }
-
+/**
+ * @brief theta_e
+ * 
+ * @tparam T 
+ * @param pressure 
+ * @param temperature 
+ * @param dewpoint 
+ * @return constexpr T 
+ */
 template <floating T>
 constexpr T equivalent_potential_temperature(
   const T pressure, const T temperature, const T dewpoint
@@ -71,7 +79,15 @@ constexpr T equivalent_potential_temperature(
 
     return th_l * exp(r * (1 + 0.448 * r) * (3036.0 / t_l - 1.78));
 }
-
+/**
+ * @brief theta_w
+ * 
+ * @tparam T 
+ * @param pressure 
+ * @param temperature 
+ * @param dewpoint 
+ * @return constexpr T 
+ */
 template <floating T>
 constexpr T wet_bulb_potential_temperature(
   const T pressure, const T temperature, const T dewpoint

src/lib/wind.cpp

@@ -2,33 +2,35 @@
 
 namespace libthermo {
 
-/* ........................................{ common }........................................... */
 template <floating T>
-constexpr T wind_direction(const T u, const T v, const bool from) noexcept {
-    T wdir = degrees(atan2(u, v));
-    if (from)
-        wdir = fmod(wdir - 180.0, 360.0);
-
-    if (wdir <= 0)
-        wdir = 360.0;
-    if (u == 0 && v == 0)
-        wdir = 0.0;
-
-    return wdir;
+constexpr T wind_direction(const T u, const T v) noexcept {
+    return fmod(degrees(atan2(u, v)) + 180.0, 360.0);
+}
+template <floating T>
+constexpr T wind_direction(const wind_components<T>& uv) noexcept {
+    return wind_direction(uv.u, uv.v);
 }
+
 template <floating T>
 constexpr T wind_magnitude(const T u, const T v) noexcept {
     return hypot(u, v);
 }
 
 template <floating T>
-constexpr WindComponents<T> wind_components(const T direction, const T magnitude) noexcept {
-    if (direction < 0 || direction > 360)
-        return {NAN, NAN};
-
-    const T u = magnitude * sin(radians(direction));
-    const T v = magnitude * cos(radians(direction));
+constexpr T wind_magnitude(const wind_components<T>& uv) noexcept {
+    return hypot(uv.u, uv.v);
+}
 
-    return {-u, -v};
+template <floating T>
+constexpr wind_vector<T>::wind_vector(const wind_components<T>& uv) noexcept :
+    direction(wind_direction(uv)), magnitude(wind_magnitude(uv)) {
 }
+template <floating T>
+constexpr wind_components<T>::wind_components(const wind_vector<T>& dm) noexcept {
+    const T d = radians(dm.direction);
+    const T m = -dm.magnitude;
+    u = m * sin(d);
+    v = m * cos(d);
+}
+
 }  // namespace libthermo

src/nzthermo/_C.pxd

@@ -36,14 +36,20 @@ cdef inline size_t search_pressure(Float[:] pressure, Float value) noexcept nogi
 
 
 cdef extern from "wind.cpp" namespace "libthermo" nogil:
-    cdef cppclass WindComponents[T]:
-        T u
-        T v
+    cdef cppclass wind_components[T]:
+        T u, v
+        wind_components() noexcept
+        wind_components(T d, T m) noexcept
+        wind_components(wind_vector[T] dm) noexcept
+
+    cdef cppclass wind_vector[T]:
+        T direction, magnitude
+        wind_vector() noexcept
+        wind_vector(T d, T s) noexcept
+        wind_vector(wind_components[T] uv)  noexcept
 
     T wind_direction[T](T u, T v) noexcept
-    T wind_direction[T](T u, T v, T from_) noexcept
     T wind_magnitude[T](T u, T v) noexcept
-    WindComponents[T] wind_components[T](T direction, T magnitude) noexcept
 
 
 cdef extern from "libthermo.cpp" namespace "libthermo" nogil:

src/nzthermo/__init__.py

@@ -5,6 +5,17 @@
     "__version__",
     # ._core
     "OPENMP_ENABLED",
+    "E0",
+    "P0",
+    "T0",
+    "Cp",
+    "Lv",
+    "Md",
+    "Mw",
+    "Rd",
+    "Rv",
+    "epsilon",
+    "kappa",
     "moist_lapse",
     "parcel_profile",
     "parcel_profile_with_lcl",
@@ -20,15 +31,16 @@
     "saturation_vapor_pressure",
     "wet_bulb_temperature",
     "wet_bulb_potential_temperature",
-    "wind_direction",
-    "wind_components",
-    "wind_magnitude",
     "dewpoint_from_specific_humidity",
     "lfc",
     "parcel_profile",
     "mixing_ratio",
     "vapor_pressure",
     "virtual_temperature",
+    "wind_components",
+    "wind_direction",
+    "wind_magnitude",
+    "wind_vector",
     # .core
     "el",
     "cape_cin",
@@ -46,7 +58,23 @@
     __version__ = "undefined"
 
 from . import functional
-from ._core import OPENMP_ENABLED, moist_lapse, parcel_profile, parcel_profile_with_lcl
+from ._core import (
+    E0,
+    OPENMP_ENABLED,
+    P0,
+    T0,
+    Cp,
+    Lv,
+    Md,
+    Mw,
+    Rd,
+    Rv,
+    epsilon,
+    kappa,
+    moist_lapse,
+    parcel_profile,
+    parcel_profile_with_lcl,
+)
 from ._ufunc import (
     delta_t,
     dewpoint,
@@ -65,6 +93,7 @@
     wind_components,
     wind_direction,
     wind_magnitude,
+    wind_vector,
 )
 from .core import (
     cape_cin,

src/nzthermo/_core.pyi

@@ -15,7 +15,19 @@ _DualArrayLike = (
 )
 _float = TypeVar("_float", np.float32, np.float64)
 
-OPENMP_ENABLED: bool
+OPENMP_ENABLED: bool = ...
+
+T0: float = ...
+E0: float = ...
+Cp: float = ...
+Rd: float = ...
+Rv: float = ...
+Lv: float = ...
+P0: float = ...
+Mw: float = ...
+Md: float = ...
+epsilon: float = ...
+kappa: float = ...
 
 @overload
 def moist_lapse[T: np.floating[Any]](

src/nzthermo/_core.pyx

@@ -71,6 +71,32 @@ ctypedef enum ProfileStrategy:
 
 OPENMP_ENABLED = bool(OPENMP)
 
+
+T0 = C.T0
+"""`(J/kg*K)` - freezing point in kelvin"""
+E0 = C.E0
+"""`(Pa)` - vapor pressure at T0"""
+Cp = C.Cp
+"""`(J/kg*K)` - specific heat of dry air"""
+Rd = C.Rd
+"""`(J/kg*K)` - gas constant for dry air"""
+Rv = C.Rv
+"""`(J/kg*K)` - gas constant for water vapor"""
+Lv = C.Lv
+"""`(J/kg)` - latent heat of vaporization"""
+P0 = C.P0
+"""`(Pa)` - standard pressure at sea level"""
+Mw = C.Mw
+"""`(g/mol)` - molecular weight of water"""
+Md = C.Md
+"""`(g/mol)` - molecular weight of dry air"""
+epsilon = C.epsilon
+"""`Mw / Md` - molecular weight ratio"""
+kappa = C.kappa
+"""`Rd / Cp`  - ratio of gas constants"""
+
+
+
 # ............................................................................................... #
 # helpers
 # ............................................................................................... #

src/nzthermo/_ufunc.pyi

@@ -25,6 +25,8 @@ def wind_direction(u: float, v: float) -> float: ...
 @_ufunc2x1
 def wind_magnitude(u: float, v: float) -> float: ...
 @_ufunc2x2
+def wind_vector(u: float, v: float) -> tuple[float, float]: ...
+@_ufunc2x2
 def wind_components(direction: float, speed: float) -> tuple[float, float]: ...
 
 # 1x1

src/nzthermo/_ufunc.pyx

@@ -70,19 +70,25 @@ cdef T wind_direction(T u, T v) noexcept nogil:
 cdef T wind_magnitude(T u, T v) noexcept nogil:
     return C.wind_magnitude(u, v)
 
-
+# TODO: the signature....
+# cdef (T, T) wind_components(T direction, T magnitude) noexcept nogil:...
+# 
+# Is unsupported by the ufunc signature. So the option are:
+# - maintain gil
+# - cast to double 
+# - write the template in C
 @cython.ufunc
-cdef (double, double) wind_components(T direction, T magnitude) noexcept nogil:
-    # TODO: the signature....
-    # cdef (T, T) wind_components(T direction, T magnitude) noexcept nogil:...
-    # 
-    # Is unsupported by the ufunc signature. So the option are:
-    # - maintain gil
-    # - cast to double 
-    # - write the template in C
-    cdef C.WindComponents[T] wnd = C.wind_components(direction, magnitude)
-    return <double>wnd.u, <double>wnd.v
+cdef (double, double) wind_components(T d, T m) noexcept nogil: 
+    cdef C.wind_components[T] uv = C.wind_components[T](C.wind_vector[T](d, m))
+    return <double>uv.u, <double>uv.v
+
 
+@cython.ufunc
+cdef (double, double) wind_vector(T u, T v) noexcept nogil:
+    cdef C.wind_vector[T] dm = C.wind_vector[T](C.wind_components[T](u, v))
+    return <double>dm.direction, <double>dm.magnitude
+
+
 
 # 1x1
 @cython.ufunc

tests/ufunc_test.py

@@ -12,12 +12,13 @@
     wet_bulb_potential_temperature,
     wet_bulb_temperature,
     wind_components,
+    wind_vector,
 )
 
 Pa = units.pascal
 K = units.kelvin
 dimensionless = units.dimensionless
-WIND_DIRECTIONS = np.array([0, 90, 180, 270, 360])
+WIND_DIRECTIONS = np.array([0, 90, 180, 270, 0])
 WIND_MAGNITUDES = np.array([10, 20, 30, 40, 50])
 
 
@@ -43,6 +44,9 @@ def test_wind_components() -> None:
         ],
     )
 
+    u, v = wind_components(WIND_DIRECTIONS, WIND_MAGNITUDES)
+    assert_allclose([WIND_DIRECTIONS, WIND_MAGNITUDES], wind_vector(u, v), atol=1e-9)
+
 
 @pytest.mark.parametrize("dtype", [np.float64, np.float32])
 def test_wet_bulb_temperature(dtype):