diff --git a/src/include/libthermo.hpp b/src/include/libthermo.hpp index da2e937..93e081a 100644 --- a/src/include/libthermo.hpp +++ b/src/include/libthermo.hpp @@ -66,12 +66,12 @@ constexpr T potential_temperature(const T pressure, const T temperature) noexcep template constexpr T equivalent_potential_temperature( const T pressure, const T temperature, const T dewpoint -) noexcept; // theta_e +) noexcept; template 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& other) noexcept = default; - // move constructor - constexpr lcl(lcl&& other) noexcept = default; - // copy assignment operator - constexpr lcl& operator=(const lcl& other) noexcept = default; - // destructor - ~lcl() noexcept = default; constexpr lcl(const T pressure, const T temperature) noexcept : pressure(pressure), temperature(temperature){}; diff --git a/src/include/wind.hpp b/src/include/wind.hpp index 48270bb..615e760 100644 --- a/src/include/wind.hpp +++ b/src/include/wind.hpp @@ -1,35 +1,39 @@ #pragma once -#include #include -#include -#include -#include -#include #include namespace libthermo { template -struct WindVector { - T speed; - T direction; -}; +class wind_vector; template -struct WindComponents { - T u; - T v; -}; +class wind_components; template -constexpr WindComponents wind_components(const T direction, const T magnitude) noexcept; +constexpr T wind_direction(const T u, const T v) noexcept; template -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 -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& uv) noexcept; +}; + +template +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& uv) noexcept; +}; } // namespace libthermo diff --git a/src/lib/libthermo.cpp b/src/lib/libthermo.cpp index ad874c7..ca08125 100644 --- a/src/lib/libthermo.cpp +++ b/src/lib/libthermo.cpp @@ -58,7 +58,15 @@ template 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 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 constexpr T wet_bulb_potential_temperature( const T pressure, const T temperature, const T dewpoint diff --git a/src/lib/wind.cpp b/src/lib/wind.cpp index e9f6950..74bf153 100644 --- a/src/lib/wind.cpp +++ b/src/lib/wind.cpp @@ -2,33 +2,35 @@ namespace libthermo { -/* ........................................{ common }........................................... */ template -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 +constexpr T wind_direction(const wind_components& uv) noexcept { + return wind_direction(uv.u, uv.v); } + template constexpr T wind_magnitude(const T u, const T v) noexcept { return hypot(u, v); } template -constexpr WindComponents 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& uv) noexcept { + return hypot(uv.u, uv.v); +} - return {-u, -v}; +template +constexpr wind_vector::wind_vector(const wind_components& uv) noexcept : + direction(wind_direction(uv)), magnitude(wind_magnitude(uv)) { } +template +constexpr wind_components::wind_components(const wind_vector& dm) noexcept { + const T d = radians(dm.direction); + const T m = -dm.magnitude; + u = m * sin(d); + v = m * cos(d); +} + } // namespace libthermo \ No newline at end of file diff --git a/src/nzthermo/_C.pxd b/src/nzthermo/_C.pxd index 8245b37..94d3603 100644 --- a/src/nzthermo/_C.pxd +++ b/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: diff --git a/src/nzthermo/__init__.py b/src/nzthermo/__init__.py index ded7071..b1af963 100644 --- a/src/nzthermo/__init__.py +++ b/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, diff --git a/src/nzthermo/_core.pyi b/src/nzthermo/_core.pyi index 8dc5ceb..b4b1290 100644 --- a/src/nzthermo/_core.pyi +++ b/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]]( diff --git a/src/nzthermo/_core.pyx b/src/nzthermo/_core.pyx index efada56..1687c67 100644 --- a/src/nzthermo/_core.pyx +++ b/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 # ............................................................................................... # diff --git a/src/nzthermo/_ufunc.pyi b/src/nzthermo/_ufunc.pyi index 5e59ab1..10899de 100644 --- a/src/nzthermo/_ufunc.pyi +++ b/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 diff --git a/src/nzthermo/_ufunc.pyx b/src/nzthermo/_ufunc.pyx index 12773bc..5dfd3f5 100644 --- a/src/nzthermo/_ufunc.pyx +++ b/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 wnd.u, 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 uv.u, 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 dm.direction, dm.magnitude + + # 1x1 @cython.ufunc diff --git a/tests/ufunc_test.py b/tests/ufunc_test.py index c44a23c..16b3de5 100644 --- a/tests/ufunc_test.py +++ b/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):