Merge pull request #23 from leaver2000/development

⛈️clean up, testing, documentation

src/include/libthermo.hpp

@@ -95,7 +95,7 @@ class lcl {
       const size_t max_iters = DEFAULT_ITERS
     ) noexcept;
     constexpr T wet_bulb_temperature(const T pressure, const T step = DEFAULT_STEP) noexcept;
-    constexpr size_t index(const T pressure[], const size_t size) noexcept;
+    constexpr size_t index_on(const T pressure[], const size_t size) noexcept;
 };
 
 template <floating T>

src/lib/functional.cpp

@@ -158,7 +158,7 @@ constexpr T rk2(Fn<T, T, T> fn, T x0, T x1, T y, T step /* 1000.0 (Pa) */) noexc
     T k1, delta, abs_delta;
     size_t N = 1;
 
-    abs_delta = fabs(delta = (x1 - x0));
+    abs_delta = std::fabs(delta = (x1 - x0));
     if (abs_delta > step) {
         N = (size_t)ceil(abs_delta / step);
         delta = delta / (T)N;
@@ -206,9 +206,9 @@ constexpr T fixed_point(
         p2 = fn(p1, x0, args...);
         delta = p2 - 2.0 * p1 + p0;
 
-        p2 = delta ? p0 - pow(p1 - p0, 2) / delta : p2; /* delta squared */
+        p2 = delta ? p0 - std::pow(p1 - p0, 2) / delta : p2; /* delta squared */
 
-        if ((p0 ? fabs((p2 - p0) / p0) : p2 /* absolute relative error */) < eps)
+        if ((p0 ? std::fabs((p2 - p0) / p0) : p2 /* absolute relative error */) < eps)
             return p2;
 
         p0 = p2;
@@ -289,8 +289,8 @@ constexpr point<T> intersect_1d(
     x0 = x[i0];
     x1 = x[i1];
     if (log_x) {
-        x0 = log(x0);
-        x1 = log(x1);
+        x0 = std::log(x0);
+        x1 = std::log(x1);
     }
 
     a0 = a[i0];
@@ -305,7 +305,7 @@ constexpr point<T> intersect_1d(
     y_intercept = ((x_intercept - x0) / LIMIT_ZERO(x1 - x0)) * (a1 - a0) + a0;
 
     if (log_x)
-        x_intercept = exp(x_intercept);
+        x_intercept = std::exp(x_intercept);
 
     return {x_intercept, y_intercept};
 };

src/lib/libthermo.cpp

@@ -2,17 +2,15 @@
 
 namespace libthermo {
 
-// helper functions
-
 /**
  * @brief given a pressure array of decreasing values, find the index of the pressure level that
  * corresponds to the given value. This function is optimized for evenly distributed pressure
  * and will typically find the index in O(1) time.
  * 
  * @tparam T 
- * @param levels 
- * @param value 
- * @param size 
+ * @param levels        pressure levels array
+ * @param value         value to find
+ * @param size          size of the array
  * @return constexpr size_t 
  */
 template <floating T>
@@ -92,34 +90,38 @@ constexpr T potential_temperature(const T pressure, const T temperature) noexcep
 }
 
 /**
- * @brief theta_e
+ * @brief Equivalent potential temperature, `theta-e (𝜃𝑒)`, is a quantity that is conserved during
+ * changes to an air parcel's pressure (that is, during vertical motions in the atmosphere), even
+ * if water vapor condenses during that pressure change. It is therefore more conserved than the
+ * ordinary potential temperature, which remains constant only for unsaturated vertical motions
+ * (pressure changes).
  * 
  * @tparam T 
- * @param pressure 
- * @param temperature 
- * @param dewpoint 
- * @return constexpr T 
+ * @param pressure          `(Pa)`
+ * @param temperature       `(K)`
+ * @param dewpoint          `(K)`
+ * @return `theta-e` -      `(K)`
  */
 template <floating T>
 constexpr T equivalent_potential_temperature(
   const T pressure, const T temperature, const T dewpoint
 ) noexcept {
     const T r = saturation_mixing_ratio(pressure, dewpoint);
     const T e = saturation_vapor_pressure(dewpoint);
-    const T t_l = 56 + 1.0 / (1.0 / (dewpoint - 56) + log(temperature / dewpoint) / 800.0);
+    const T t_l = 56. + 1. / (1. / (dewpoint - 56.) + log(temperature / dewpoint) / 800.);
     const T th_l =
       potential_temperature(pressure - e, temperature) * pow(temperature / t_l, 0.28 * r);
 
-    return th_l * exp(r * (1 + 0.448 * r) * (3036.0 / t_l - 1.78));
+    return th_l * exp(r * (1. + 0.448 * r) * (3036. / t_l - 1.78));
 }
 
 /**
  * @brief theta_w
  * 
  * @tparam T 
- * @param pressure 
- * @param temperature 
- * @param dewpoint 
+ * @param pressure          `(Pa)`
+ * @param temperature       `(K)`
+ * @param dewpoint          `(K)`
  * @return constexpr T 
  */
 template <floating T>
@@ -188,7 +190,7 @@ constexpr T lcl<T>::wet_bulb_temperature(const T pressure, const T step) noexcep
 }
 
 template <floating T>
-constexpr size_t lcl<T>::index(const T pressure[], const size_t size) noexcept {
+constexpr size_t lcl<T>::index_on(const T pressure[], const size_t size) noexcept {
     return index_pressure(pressure, this->pressure, size);
 }
 

src/nzthermo/_C.pxd

@@ -54,7 +54,7 @@ cdef extern from "libthermo.cpp" namespace "libthermo" nogil:
         lcl() noexcept
         lcl(T pressure, T temperature) noexcept
         lcl(T pressure, T temperature, T dewpoint) noexcept
-        size_t index(T* levels, size_t size) noexcept
+        size_t index_on(T* levels, size_t size) noexcept
 
     # 1x1
     T saturation_vapor_pressure[T](T temperature) noexcept

src/nzthermo/_core.pyi

@@ -14,7 +14,6 @@ _Dtype_t = TypeVar("_Dtype_t", bound=np.generic)
 _DualArrayLike = (
     SupportsArray[_Dtype_t] | NestedSequence[SupportsArray[_Dtype_t]] | _T | NestedSequence[_T]
 )
-_float = TypeVar("_float", np.float32, np.float64)
 
 OPENMP_ENABLED: bool = ...
 

src/nzthermo/_core.pyx

@@ -103,16 +103,6 @@ kappa = C.kappa
 # ............................................................................................... #
 # helpers
 # ............................................................................................... #
-cdef cvarray nzarray((size_t, size_t) shape, size_t size):
-    return cvarray(
-        shape, 
-        itemsize=size, 
-        format='f' if size == sizeof(float) else 'd', 
-        mode='c',
-        allocate_buffer=True,
-    )
-
-
 def broadcast_and_nanmask(
     np.ndarray pressure not None,
     np.ndarray temperature not None, 
@@ -284,7 +274,7 @@ cdef T[:, :] moist_lapse_2d(
         T[:, :] out
 
     N, Z = temperature.shape[0], pressure.shape[1]
-    out = nzarray((N, Z), pressure.itemsize)
+    out = np.empty((N, Z), dtype=np.dtype(f"f{pressure.itemsize}"))
     with nogil, parallel():
         if BROADCAST is mode:
             for i in prange(N, schedule='dynamic'):
@@ -466,7 +456,7 @@ cdef void parcel_profile_1d(
 
     # [dry ascent]
     # stop the dry ascent at the LCL
-    stop = lcl.index(&pressure[0], Z)
+    stop = lcl.index_on(&pressure[0], Z)
     for i in prange(0, stop, schedule='dynamic'): # parallelize the dry ascent
         out[i] = C.dry_lapse(pressure[i], p0, temperature)
 
@@ -549,8 +539,7 @@ cdef void parcel_profile_with_lcl_1d(
     lcl = C.lcl[T](p0, t0, td0)
 
     # [dry ascent] .. parcel temperature from the surface up to the LCL ..
-
-    stop = lcl.index(&pressure[0], Z)
+    stop = lcl.index_on(&pressure[0], Z)
 
     if stop:
         ep[:stop] = pressure[:stop]
@@ -774,7 +763,7 @@ cdef intersect_2d(
         T[:, :] out 
 
     N, Z = y0.shape[0], y1.shape[1]
-    out = nzarray((2, N), x.itemsize)
+    out = np.empty((2, N), f'f{x.itemsize}')
     with nogil, parallel():
         if BROADCAST is mode:
             for i in prange(N, schedule='dynamic'):

src/nzthermo/_ufunc.pyi

@@ -13,9 +13,9 @@ _DType_T_co = TypeVar("_DType_T_co", bound=np.dtype[Any])
 
 class pressure_vector(np.ndarray[_S, _DType_T_co]):
     @overload
-    def __new__(cls, array: np.ndarray[_S, _DType_T_co]) -> pressure_vector[_S, _DType_T_co]: ...
+    def __new__(cls, __x: np.ndarray[_S, _DType_T_co], /) -> pressure_vector[_S, _DType_T_co]: ...
     @overload
-    def __new__(cls, array: ArrayLike) -> pressure_vector[Any, np.dtype[Any]]: ...
+    def __new__(cls, __x: ArrayLike, /) -> pressure_vector[Any, np.dtype[Any]]: ...
     def is_above(
         self, bottom: ArrayLike, close: bool = ...
     ) -> np.ndarray[_S, np.dtype[np.bool_]]: ...