⛈️

README.md

@@ -1,5 +1,33 @@
 # nzthermo
 
+This work has been heavily inspired by the excellent work and code base that has been developed by
+the `MetPy` team. The concept of `(N, Z)` is simply to solve for `N` profiles of `Z` levels. So
+regardless of what what you data looks like, if it can be reshaped to `(N, Z)` then it can be used
+with this library. Where possible all iterations of `N` are run in parallel using `OpenMP` and
+`Cython`.  Most of the root functions are written in `C++` and wrapped with `Cython` for use in
+`Python`.
+
+Where this code currently lacks in complete documentation it makes up for with the extensive and
+verbose usage of type annotations. For example, the `parcel_profile` function is defined as follows:
+
+```python
+def parcel_profile(
+    pressure: Pascal[np.ndarray[shape[Z], np.dtype[T]] | np.ndarray[shape[N, Z], np.dtype[T]]],
+    temperature: Kelvin[np.ndarray[shape[N], np.dtype[np.float_]]],
+    dewpoint: Kelvin[np.ndarray[shape[N], np.dtype[np.float_]]],
+    /,
+    *,
+    step: float = ...,
+    eps: float = ...,
+    max_iters: int = ...,
+) -> Kelvin[np.ndarray[shape[N, Z], np.dtype[T]]]: ...
+```
+
+Which make it quite clean that the pressure array is expected to be of shape `(Z,)` or `(N, Z)` and
+have the units of `Pascal`. The temperature and dewpoint arrays are expected to be of shape `(N,)`
+and have the units of `Kelvin`. The return value is expected to be of shape `(N, Z)` and have the
+units of `Kelvin`.
+
 ## Getting Started
 
 ```bash

pyproject.toml

@@ -74,7 +74,7 @@ exclude_lines = [
 
 [tool.black]
 target-version = ['py312']
-line-length = 119
+line-length = 99
 include = '(nzthermo|tests)\/.*(\.py|\.pyi)'
 force-exclude = '(nzthermo|tests)\/.*(\.pyx)'
 
@@ -91,14 +91,14 @@ combine_as_imports = true
 
 [tool.ruff]
 target-version = "py312"
-line-length = 119
+line-length = 99
 fix = true
 
-# [tool.ruff.lint]
-# ignore = [
-#     "E731", # do not assign a lambda expression, use a def
-#     "E402", # module level import not at top of file
-#     "E402", # module level import not at top of file
-#     "E501", # line too long
-#     "E741", # do not use variables named 'l', 'O', or 'I'
-# ]
+[tool.ruff.lint]
+ignore = [
+    "E731", # do not assign a lambda expression, use a def
+    "E402", # module level import not at top of file
+    "E402", # module level import not at top of file
+    "E501", # line too long
+    "E741", # do not use variables named 'l', 'O', or 'I'
+]

src/include/common.hpp

@@ -4,10 +4,18 @@
 
 namespace libthermo {
 
+#define LIMIT_ZERO(x) ((x) ? (x) : std::numeric_limits<T>::epsilon())
+
 template <typename T>
 concept floating = std::is_floating_point_v<T>;
 
 template <typename R, typename... Args>
 using Fn = R (*)(Args...);
 
+template <floating T>
+struct point {
+    T x;
+    T y;
+};
+
 }  // namespace libthermo

src/include/functional.hpp

@@ -7,6 +7,7 @@
 
 #include <common.hpp>
 #include <functional.hpp>
+
 namespace libthermo {
 
 typedef enum direction {
@@ -36,13 +37,9 @@ size_t search_sorted(
   const T x[], const T value, const size_t size, const bool inverted = false
 ) noexcept;
 
-template <floating T>
-constexpr T interpolate_z(const T x, const T xp[], const T fp[], const size_t size) noexcept;
-
 /**
  * see: https://github.com/numpy/numpy/blob/main/numpy/_core/src/npymath/npy_math_internal.h.src#L426
  */
-
 template <floating T>
 constexpr T heaviside(const T x, const T h0) noexcept;
 
@@ -85,7 +82,7 @@ constexpr T fixed_point(
  * Based on std::lower_bound, this iterates over an array using a binary search
  * to find the first element that does not satisfy the comparison condition. By
  * default, the comparison is std::less. Binary search expects data to be sorted
- * in ascending order -- for pressure level data, change the comparitor.
+ * in ascending order -- for pressure level data, change the comparator.
  *
  * We use a custom implementation of sharp::lower_bound rather than
  * std::lower_bound for a few reasons. First, we prefer raw pointer
@@ -136,6 +133,16 @@ template <floating T, typename C = std::less<T>>
 constexpr size_t upper_bound(const T array[], const int N, const T& value, const C cmp = C{});
 
 template <floating T>
-constexpr T trapezoidal(const Fn<T, T> fn, const T a, const T b, const T n) noexcept;
+constexpr T interpolate_1d(const T x, const T xp[], const T fp[], const size_t size) noexcept;
 
+template <floating T>
+constexpr point<T> intersect_1d(
+  const T x[],
+  const T a[],
+  const T b[],
+  const size_t size,
+  const bool log_x = false,
+  const bool increasing = true,
+  const bool bottom = true
+) noexcept;
 }  // namespace libthermo

src/include/libthermo.hpp

@@ -21,7 +21,7 @@ static constexpr double Rv = 461.52311572606084;  // `(J/kg*K)` - gas constant f
 static constexpr double Lv = 2501000.0;  // `(J/kg)` - latent heat of vaporization
 static constexpr double P0 = 100000.0;  // `(Pa)` - standard pressure at sea level
 static constexpr double Mw = 18.01528;  // `(g/mol)` - molecular weight of water
-static constexpr double Md = 28.96546e-3 * 1e3;  // `(g/mol)` - molecular weight of dry air
+static constexpr double Md = 28.96546;  // `(g/mol)` - molecular weight of dry air
 static constexpr double epsilon = Mw / Md;  // `Mw / Md` - molecular weight ratio
 static constexpr double kappa = Rd / Cp;  // `Rd / Cp`  - ratio of gas constants
 

src/lib/functional.cpp

@@ -4,6 +4,24 @@
 
 namespace libthermo {
 
+/**
+ * @brief The sign function returns -1 if x < 0, 0 if x==0, 1 if x > 0. nan is returned for nan inputs.
+ * 
+ * @tparam T 
+ * @param x 
+ * @return constexpr T 
+ */
+template <floating T>
+constexpr T sign(const T x) noexcept {
+    if (isnan(x))
+        return NAN;
+    else if (x == 0.)
+        return 0.;
+    else if (x < 0)
+        return -1.;
+    return 1.;
+}
+
 template <floating T>
 constexpr bool monotonic(const T x[], const size_t size, direction direction) noexcept {
     if (direction == direction::increasing) {
@@ -31,14 +49,14 @@ constexpr T radians(const T degrees) noexcept {
 
 template <floating T>
 constexpr T norm(const T x, const T x0, const T x1) noexcept {
-    return (x - x0) / (x1 - x0);
+    return (x - x0) / LIMIT_ZERO(x1 - x0);
 }
 
 template <floating T>
 constexpr T linear_interpolate(
   const T x, const T x0, const T x1, const T y0, const T y1
 ) noexcept {
-    return y0 + (x - x0) * (y1 - y0) / (x1 - x0);
+    return y0 + (x - x0) * (y1 - y0) / LIMIT_ZERO(x1 - x0);
 }
 
 template <floating T, typename C>
@@ -73,15 +91,6 @@ size_t search_sorted(const T x[], const T value, const size_t size, const bool i
     return upper_bound(x, size, value, std::less_equal());
 }
 
-template <floating T>
-constexpr T interpolate_z(const T x, const T xp[], const T fp[], const size_t size) noexcept {
-    const size_t i = lower_bound(xp, size, x, std::greater_equal());
-    if (i == 0)
-        return fp[0];
-
-    return linear_interpolate(x, xp[i - 1], xp[i], fp[i - 1], fp[i]);
-}
-
 template <floating T>
 constexpr T heaviside(const T x, const T h0) noexcept {
     if (isnan(x))
@@ -130,12 +139,9 @@ constexpr T fixed_point(
         p1 = fn(p0, x0, args...);
         p2 = fn(p1, x0, args...);
         delta = p2 - 2.0 * p1 + p0;
-        if (delta)
-            p2 = p0 - pow(p1 - p0, 2) / delta; /* delta squared */
 
-        err = p2;
-        if (p0)
-            err = fabs((p2 - p0) / p0); /* absolute relative error */
+        p2 = delta ? p0 - pow(p1 - p0, 2) / delta : p2; /* delta squared */
+        err = p0 ? fabs((p2 - p0) / p0) : p2; /* absolute relative error */
 
         if (err < eps)
             return p2;
@@ -147,36 +153,71 @@ constexpr T fixed_point(
 }
 
 template <floating T>
-constexpr T trapezoidal(const Fn<T, T> fn, const T a, const T b, const T n) noexcept {
-    // Grid spacing
-    T h = (b - a) / n;
-
-    // Computing sum of first and last terms
-    // in above formula
-    T s = fn(a) + fn(b);
-
-    // Adding middle terms in above formula
-    for (size_t i = 1; i < n; i++)
-        s += 2 * fn(a + i * h);
+constexpr T interpolate_1d(const T x, const T xp[], const T fp[], const size_t size) noexcept {
+    const size_t i = lower_bound(xp, size, x, std::greater_equal());
+    if (i == 0)
+        return fp[0];
 
-    // h/2 indicates (b-a)/2n. Multiplying h/2
-    // with s.
-    return (h / 2) * s;
+    return linear_interpolate(x, xp[i - 1], xp[i], fp[i - 1], fp[i]);
 }
 
 template <floating T>
-constexpr T* find_intersections(T x[], T a[], T b[], size_t size, bool log_x) {
-    T* intersections = new T[size];
+constexpr point<T> intersect_1d(
+  const T x[],
+  const T a[],
+  const T b[],
+  const size_t size,
+  const bool log_x,
+  const bool increasing,
+  const bool bottom
+) noexcept {
+    const T sign_check = increasing ? 1. : -1.;
+    T x0, x1, a0, a1, b0, b1, delta_y0, delta_y1, x_intercept, y_intercept;
+    size_t i0, i1;
+    if (bottom) {
+        for (i0 = 0, i1 = 1; i0 < size - 1; i0++, i1++) {
+            T s0 = sign(a[i0] - b[i0]);
+            T s1 = sign(a[i1] - b[i1]);
+
+            if ((s0 != s1) && (s1 == sign_check))
+                break;
+        }
+        if (i1 == size)
+            return {NAN, NAN};
+    } else {
+        for (i0 = size - 1, i1 = size; i0 > 0; i0--, i1--) {
+            T s0 = sign(a[i0] - b[i0]);
+            T s1 = sign(a[i1] - b[i1]);
+
+            if ((s0 != s1) && (s1 == sign_check))
+                break;
+        }
+        if (i0 == 0)
+            return {NAN, NAN};
+    }
+
+    x0 = x[i0];
+    x1 = x[i1];
+    if (log_x) {
+        x0 = log(x0);
+        x1 = log(x1);
+    }
+
+    a0 = a[i0];
+    a1 = a[i1];
+    b0 = b[i0];
+    b1 = b[i1];
+
+    delta_y0 = a0 - b0;
+    delta_y1 = a1 - b1;
 
-    size_t nearest_idx;
-    size_t sign_change_idx;
-    //     difference = a - b
+    x_intercept = (delta_y1 * x0 - delta_y0 * x1) / LIMIT_ZERO(delta_y1 - delta_y0);
+    y_intercept = ((x_intercept - x0) / LIMIT_ZERO(x1 - x0)) * (a1 - a0) + a0;
 
-    // # Determine the point just before the intersection of the lines
-    // # Will return multiple points for multiple intersections
-    // sign_change_idx, = np.nonzero(np.diff(np.sign(difference)))
+    if (log_x)
+        x_intercept = exp(x_intercept);
 
-    return intersections;
+    return {x_intercept, y_intercept};
 };
 
 }  // namespace libthermo

src/lib/libthermo.cpp

@@ -146,75 +146,6 @@ constexpr T wet_bulb_temperature(
     return moist_lapse(x.pressure, pressure, x.temperature, step);
 }
 
-// template <floating T>
-// constexpr T cape_cin(
-//   const T pressure[], const T temperature[], const T dewpoint[], const size_t size
-// ) noexcept {
-//     return 0.0;
-// }
-
-// template <floating T>
-// constexpr T downdraft_cape(
-//   const T pressure[], const T temperature[], const T dewpoint[], const size_t size
-// ) noexcept {
-//     T p, t, td, delta, p_top, t_top, td_top, wb_top, trace, ln, delta_start;
-//     size_t start, stop;
-//     start = 0;
-//     stop = size - 1;
-//     for (size_t i = 0; i < size; i++) {
-//         if (pressure[i] >= 7e4) {  // start a 7000. Pa
-//             continue;
-
-//         } else if (start == 0) {
-//             start = i;
-//         }
-
-//         if (pressure[i] <= 5e4) {  // stop at 5000. Pa
-//             stop = i;
-//             break;
-//         }
-//     }
-//     p_top = pressure[stop];  //  # (N,)
-//     t_top = temperature[stop];  // # (N,)
-//     td_top = dewpoint[stop];  //];  # (N,)
-
-//     wb_top = wet_bulb_temperature(p_top, t_top, td_top, .1, 1000.0, 50);  // # (N,)
-
-//     p = pressure[start];
-//     t = temperature[start];
-//     td = dewpoint[start];
-
-//     trace = moist_lapse(p, wb_top, p_top, 1000.0);
-//     delta_start = virtual_temperature(p, t, td) - virtual_temperature(p, trace, trace);
-//     T ln_start = std::log(p);
-
-//     T dcape = 0.0;
-//     for (size_t i = start + 1; i < stop; i++) {
-//         p = pressure[i];
-//         t = temperature[i];
-//         td = dewpoint[i];
-//         trace = moist_lapse(p, wb_top, p_top, 1000.0);
-//         dcape -= virtual_temperature(p, t, td) - virtual_temperature(p, trace, trace);
-
-//         // ln = std::log(p);
-//         // // for (size_t j = 1; j < 5; j++)
-//         // //     delta += 2 * virtual_temperature(pressure[j], t, td) -
-//         // //       virtual_temperature(p, trace, trace);  //fn(a + i * h);
-//         // // dcape += trapezoidal(delta, ln, (T)5.0);
-//         // // dcape += (delta - delta_start) * (ln + ln_start) / 2.0;
-//         // // // next iter
-//         // ln_start = ln;
-//         // delta_start = delta;
-//     }
-//     // dcape = -(Rd * F.nantrapz(delta, logp, axis=1))
-//     return -Rd * dcape;
-// }
-
-/* ........................................{ sharp  }...........................................
-
-
-   ........................................{ sharp  }........................................... */
-
 template <floating T>
 constexpr T wobus(T temperature) {
     T pol;