diff --git a/PySDM/physics/terminal_velocity/__init__.py b/PySDM/physics/terminal_velocity/__init__.py
index 2e012dc474..3434f7ad3a 100644
--- a/PySDM/physics/terminal_velocity/__init__.py
+++ b/PySDM/physics/terminal_velocity/__init__.py
@@ -1,4 +1,6 @@
-"""terminal velocity formulae"""
+"""
+terminal velocity formulae
+"""
from .rogers_yau import RogersYau
from .gunn_kinzer_1949 import GunnKinzer1949
diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py
new file mode 100644
index 0000000000..2d0f59fb30
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py
@@ -0,0 +1,8 @@
+"""
+The Physics of Falling Raindrops in Diverse Planetary Atmospheres
+[Loftus, K., & Wordsworth, R. D. 2021 (Journal of Geophysical Research: Planets, 126)](https://doi.org/10.1029/2020JE006653)
+
+"""
+
+from .settings import Settings
+from .simulation import Simulation
diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb
new file mode 100644
index 0000000000..6a27e3ffa6
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb
@@ -0,0 +1,14723 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1e8d983b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "!pip install joblib"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "d8f644e9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "from numba import njit\n",
+ "from joblib import Parallel, delayed\n",
+ "from open_atmos_jupyter_utils import show_plot\n",
+ "\n",
+ "from PySDM import Formulae\n",
+ "from PySDM.physics import si\n",
+ "from scipy.optimize import fsolve\n",
+ "from PySDM_examples.Loftus_and_Wordsworth_2021 import Settings\n",
+ "from PySDM_examples.Loftus_and_Wordsworth_2021.planet import (\n",
+ " Planet,\n",
+ " EarthLike,\n",
+ " Earth,\n",
+ " EarlyMars,\n",
+ " Jupiter,\n",
+ " Saturn,\n",
+ " K2_18B,\n",
+ ")\n",
+ "from PySDM_examples.Loftus_and_Wordsworth_2021 import Simulation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "92a3a574",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "formulae = Formulae(\n",
+ " ventilation=\"PruppacherAndRasmussen1979\",\n",
+ " saturation_vapour_pressure=\"AugustRocheMagnus\",\n",
+ " diffusion_coordinate=\"WaterMassLogarithm\",\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "41f0ed6d",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "300.0"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "new_Earth = EarthLike()\n",
+ "new_Earth.T_STP"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "d3a8f4b9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "radius_array = np.logspace(-4.5, -2.5, 50) * si.m\n",
+ "RH_array = np.linspace(0.25, 0.99, 50)\n",
+ "output_matrix = np.full((len(RH_array), len(radius_array)), np.nan)\n",
+ "const = formulae.constants\n",
+ "\n",
+ "\n",
+ "@njit()\n",
+ "def mix(dry, vap, ratio):\n",
+ " return (dry + ratio * vap) / (1 + ratio)\n",
+ "\n",
+ "\n",
+ "def compute_one_RH(i, RH):\n",
+ " \"\"\"\n",
+ " Compute one row of the output_matrix for a given RH.\n",
+ " Returns a 1D numpy array of length len(radius_array).\n",
+ " \"\"\"\n",
+ " new_Earth.RH_zref = RH\n",
+ "\n",
+ " pvs = formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP)\n",
+ " initial_water_vapour_mixing_ratio = const.eps / (\n",
+ " new_Earth.p_STP / new_Earth.RH_zref / pvs - 1\n",
+ " )\n",
+ "\n",
+ " Rair = mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio)\n",
+ " c_p = mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio)\n",
+ "\n",
+ " def f(x):\n",
+ " return initial_water_vapour_mixing_ratio / (\n",
+ " initial_water_vapour_mixing_ratio + const.eps\n",
+ " ) * new_Earth.p_STP * (x / new_Earth.T_STP) ** (\n",
+ " c_p / Rair\n",
+ " ) - formulae.saturation_vapour_pressure.pvs_water(\n",
+ " x\n",
+ " )\n",
+ "\n",
+ " tdews = fsolve(f, [150, 300])\n",
+ " Tcloud = np.max(tdews)\n",
+ " Zcloud = (new_Earth.T_STP - Tcloud) * c_p / new_Earth.g_std\n",
+ " thstd = formulae.trivia.th_std(new_Earth.p_STP, new_Earth.T_STP)\n",
+ "\n",
+ " pcloud = formulae.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio(\n",
+ " new_Earth.p_STP, thstd, initial_water_vapour_mixing_ratio, Zcloud\n",
+ " )\n",
+ "\n",
+ " np.testing.assert_approx_equal(\n",
+ " actual=pcloud\n",
+ " * (\n",
+ " initial_water_vapour_mixing_ratio\n",
+ " / (initial_water_vapour_mixing_ratio + const.eps)\n",
+ " )\n",
+ " / formulae.saturation_vapour_pressure.pvs_water(Tcloud),\n",
+ " desired=1,\n",
+ " significant=4,\n",
+ " )\n",
+ "\n",
+ " output = None\n",
+ " row_data = np.full(len(radius_array), np.nan)\n",
+ " for j, r in enumerate(radius_array[::-1]):\n",
+ " settings = Settings(\n",
+ " planet=new_Earth,\n",
+ " r_wet=r,\n",
+ " mass_of_dry_air=1e5 * si.kg,\n",
+ " initial_water_vapour_mixing_ratio=initial_water_vapour_mixing_ratio,\n",
+ " pcloud=pcloud,\n",
+ " Zcloud=Zcloud,\n",
+ " Tcloud=Tcloud,\n",
+ " formulae=formulae,\n",
+ " )\n",
+ " simulation = Simulation(settings)\n",
+ " try:\n",
+ " output = simulation.run()\n",
+ " if output[\"z\"][-1] > 0:\n",
+ " row_data[j] = np.nan\n",
+ " break\n",
+ " else:\n",
+ " row_data[j] = 1 - (output[\"r\"][-1] / (r * 1e6))\n",
+ " except Exception as _:\n",
+ " break\n",
+ "\n",
+ " return i, row_data, output"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "4352de81",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "all_rows = Parallel(n_jobs=os.cpu_count())(\n",
+ " delayed(compute_one_RH)(i, RH) for i, RH in enumerate(RH_array[::-1])\n",
+ ")\n",
+ "\n",
+ "last_output = None\n",
+ "for i, row_data, output in all_rows:\n",
+ " output_matrix[i] = row_data\n",
+ " last_output = output"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "8e6027d8",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "53496ecd340048b79ff90faa9da71b6e",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\".\\\\tmp_nham5pj.pdf
\"), HTML(val…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "row_data = output_matrix[18, ::-1] # Reverse the row for plotting\n",
+ "fig, ax = plt.subplots(\n",
+ " 2,\n",
+ " 1,\n",
+ " figsize=(6, 6),\n",
+ " sharex=True,\n",
+ " gridspec_kw={\"height_ratios\": [1, 3]},\n",
+ " constrained_layout=True,\n",
+ ")\n",
+ "ax[0].plot(radius_array, row_data, label=f\"Surface RH = {RH_array[-18]:.2f} %\")\n",
+ "ax[0].set_ylabel(\"Mass evaporated (%)\")\n",
+ "ax[0].legend()\n",
+ "\n",
+ "h = ax[1].contourf(\n",
+ " radius_array,\n",
+ " RH_array[::-1],\n",
+ " output_matrix[:, ::-1],\n",
+ " cmap=\"gray_r\",\n",
+ " levels=np.linspace(0, 1, 100),\n",
+ ")\n",
+ "ax[1].set_xscale(\"log\")\n",
+ "\n",
+ "# Add labels\n",
+ "ax[1].set_xlabel(\"Radius (µm)\")\n",
+ "ax[1].set_ylabel(\"Surface RH (%)\")\n",
+ "\n",
+ "cbar = fig.colorbar(h, ax=ax, shrink=0.4)\n",
+ "cbar.set_label(\"fraction Mass evaporated\")\n",
+ "contour_levels = [0.1] # Define the level for the contour\n",
+ "ax[1].contour(\n",
+ " radius_array,\n",
+ " RH_array[::-1],\n",
+ " output_matrix[:, ::-1],\n",
+ " levels=contour_levels,\n",
+ " colors=\"red\",\n",
+ " linewidths=1.5,\n",
+ ")\n",
+ "show_plot()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "369fa24a",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "0baa5c73262744e4890c68dcfaa37b13",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\".\\\\tmps5jrn8sn.pdf
\"), HTML(val…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "fig, axs = plt.subplots(1, 3, figsize=(10, 8), sharey=True)\n",
+ "axs[0].plot(last_output[\"r\"], last_output[\"z\"], label=\"Radius\")\n",
+ "axs[0].set_ylabel(\"Height (m)\")\n",
+ "axs[0].set_xlabel(\"Radius (um)\")\n",
+ "axs[0].legend()\n",
+ "axs[1].plot(last_output[\"S\"], last_output[\"z\"], label=\"Supersaturation\")\n",
+ "axs[1].set_xlabel(\"Supersaturation\")\n",
+ "axs[1].set_ylabel(\"Height (m)\")\n",
+ "axs[1].legend()\n",
+ "axs[2].plot(last_output[\"t\"], last_output[\"z\"], label=\"Height\")\n",
+ "axs[2].set_ylabel(\"Height (m)\")\n",
+ "axs[2].set_xlabel(\"Time (s)\")\n",
+ "axs[2].legend()\n",
+ "show_plot()"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": ".venv",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.13.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py
new file mode 100644
index 0000000000..7ac5eb0dac
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py
@@ -0,0 +1,62 @@
+from PySDM.environments.parcel import Parcel
+
+
+class AlienParcel(Parcel):
+ def __init__(
+ self,
+ dt,
+ mass_of_dry_air: float,
+ pcloud: float,
+ initial_water_vapour_mixing_ratio: float,
+ Tcloud: float,
+ w: float = 0,
+ zcloud: float = 0,
+ mixed_phase=False,
+ ):
+ super().__init__(
+ dt=dt,
+ mass_of_dry_air=mass_of_dry_air,
+ p0=pcloud,
+ initial_water_vapour_mixing_ratio=initial_water_vapour_mixing_ratio,
+ T0=Tcloud,
+ w=w,
+ z0=zcloud,
+ mixed_phase=mixed_phase,
+ variables=None,
+ )
+
+ def advance_parcel_vars(self):
+ """
+ Compute new values of displacement, dry-air density and volume,
+ and write them to self._tmp and self.mesh.dv
+ """
+ dt = self.particulator.dt
+ formulae = self.particulator.formulae
+ T = self["T"][0]
+ p = self["p"][0]
+
+ dz_dt = -self.particulator.attributes["terminal velocity"].to_ndarray()[0]
+ water_vapour_mixing_ratio = (
+ self["water_vapour_mixing_ratio"][0]
+ - self.delta_liquid_water_mixing_ratio / 2
+ )
+
+ drho_dz = formulae.hydrostatics.drho_dz(
+ p=p,
+ T=T,
+ water_vapour_mixing_ratio=water_vapour_mixing_ratio,
+ lv=formulae.latent_heat_vapourisation.lv(T),
+ d_liquid_water_mixing_ratio__dz=(
+ self.delta_liquid_water_mixing_ratio / dz_dt / dt
+ ),
+ )
+ drhod_dz = drho_dz
+
+ self.particulator.backend.explicit_euler(self._tmp["z"], dt, dz_dt)
+ self.particulator.backend.explicit_euler(
+ self._tmp["rhod"], dt, dz_dt * drhod_dz
+ )
+
+ self.mesh.dv = formulae.trivia.volume_of_density_mass(
+ (self._tmp["rhod"][0] + self["rhod"][0]) / 2, self.mass_of_dry_air
+ )
diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py
new file mode 100644
index 0000000000..0e4211b397
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py
@@ -0,0 +1,118 @@
+from PySDM.physics.constants import si
+from dataclasses import dataclass
+from typing import Optional, Dict, Any
+
+
+@dataclass
+class Planet:
+ g_std: float
+ T_STP: float
+ p_STP: float
+ RH_zref: float
+ dry_molar_conc_H2: float
+ dry_molar_conc_He: float
+ dry_molar_conc_N2: float
+ dry_molar_conc_O2: float
+ dry_molar_conc_CO2: float
+ H_LCL: Optional[float] = None
+
+ def to_dict(self) -> Dict[str, Any]:
+ return vars(self)
+
+
+@dataclass
+class EarthLike(Planet):
+ g_std: float = 9.82 * si.metre / si.second**2
+ T_STP: float = 300 * si.kelvin
+ p_STP: float = 1.01325 * 1e6 * si.pascal
+ RH_zref: float = 0.75
+ dry_molar_conc_H2: float = 0
+ dry_molar_conc_He: float = 0
+ dry_molar_conc_N2: float = 1
+ dry_molar_conc_O2: float = 0
+ dry_molar_conc_CO2: float = 0
+ H_LCL: float = 8.97 * si.kilometre
+
+
+@dataclass
+class Earth(Planet):
+ g_std: float = 9.82 * si.metre / si.second**2
+ T_STP: float = 290 * si.kelvin
+ p_STP: float = 1.01325 * 1e6 * si.pascal
+ RH_zref: float = 0.75
+ dry_molar_conc_H2: float = 0
+ dry_molar_conc_He: float = 0
+ dry_molar_conc_N2: float = 0.8
+ dry_molar_conc_O2: float = 0.2
+ dry_molar_conc_CO2: float = 0
+ H_LCL: float = 8.41 * si.kilometre
+
+
+@dataclass
+class EarlyMars(Planet):
+ g_std: float = 3.71 * si.metre / si.second**2
+ T_STP: float = 290 * si.kelvin
+ p_STP: float = 2 * 1e6 * si.pascal
+ RH_zref: float = 0.75
+ dry_molar_conc_H2: float = 0
+ dry_molar_conc_He: float = 0
+ dry_molar_conc_N2: float = 0
+ dry_molar_conc_O2: float = 0
+ dry_molar_conc_CO2: float = 1
+ H_LCL: float = 14.5 * si.kilometre
+
+
+@dataclass
+class Jupiter(Planet):
+ g_std: float = 24.84 * si.metre / si.second**2
+ T_STP: float = 274 * si.kelvin
+ p_STP: float = 4.85 * 1e6 * si.pascal
+ RH_zref: float = 1
+ dry_molar_conc_H2: float = 0.864
+ dry_molar_conc_He: float = 0.136
+ dry_molar_conc_N2: float = 0
+ dry_molar_conc_O2: float = 0
+ dry_molar_conc_CO2: float = 0
+ H_LCL: float = 39.8 * si.kilometre
+
+
+@dataclass
+class Saturn(Planet):
+ g_std: float = 10.47 * si.metre / si.second**2
+ T_STP: float = 284 * si.kelvin
+ p_STP: float = 10.4 * 1e6 * si.pascal
+ RH_zref: float = 1
+ dry_molar_conc_H2: float = 0.88
+ dry_molar_conc_He: float = 0.12
+ dry_molar_conc_N2: float = 0
+ dry_molar_conc_O2: float = 0
+ dry_molar_conc_CO2: float = 0
+ H_LCL: float = 99.2 * si.kilometre
+
+
+@dataclass
+class K2_18B(Planet):
+ g_std: float = 12.44 * si.metre / si.second**2
+ T_STP: float = 275 * si.kelvin
+ p_STP: float = 0.1 * 1e6 * si.pascal
+ RH_zref: float = 1
+ dry_molar_conc_H2: float = 0.9
+ dry_molar_conc_He: float = 0.1
+ dry_molar_conc_N2: float = 0
+ dry_molar_conc_O2: float = 0
+ dry_molar_conc_CO2: float = 0
+ H_LCL: float = 56.6 * si.kilometre
+
+
+@dataclass
+class CompositeTest(Planet):
+ g_std: float = 9.82 * si.metre / si.second**2
+ T_STP: float = 275 * si.kelvin
+ p_STP: float = 0.75 * 1e6 * si.pascal
+ RH_zref: float = 1
+ dry_molar_conc_H2: float = 0.1
+ dry_molar_conc_He: float = 0.1
+ dry_molar_conc_N2: float = 0.1
+ dry_molar_conc_O2: float = 0.1
+ dry_molar_conc_CO2: float = 0.1
+ H_LCL: Optional[float] = None
diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py
new file mode 100644
index 0000000000..b0a81752a7
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py
@@ -0,0 +1,46 @@
+# Planetary Properties, Loftus and Wordsworth 2021 Table 1
+
+from pystrict import strict
+
+from PySDM import Formulae
+from PySDM.dynamics import condensation
+from PySDM.physics.constants import si
+from PySDM_examples.Loftus_and_Wordsworth_2021.planet import Planet
+
+
+@strict
+class Settings:
+ def __init__(
+ self,
+ r_wet: float,
+ mass_of_dry_air: float,
+ planet: Planet,
+ initial_water_vapour_mixing_ratio: float,
+ pcloud: float,
+ Zcloud: float,
+ Tcloud: float,
+ formulae: Formulae = None,
+ ):
+ self.formulae = formulae or Formulae(
+ saturation_vapour_pressure="AugustRocheMagnus",
+ diffusion_coordinate="WaterMassLogarithm",
+ )
+
+ self.initial_water_vapour_mixing_ratio = initial_water_vapour_mixing_ratio
+ self.p0 = planet.p_STP
+ self.RH0 = planet.RH_zref
+ self.kappa = 0.2
+ self.T0 = planet.T_STP
+ self.z_half = 150 * si.metres
+ self.dt = 1 * si.second
+ self.pcloud = pcloud
+ self.Zcloud = Zcloud
+ self.Tcloud = Tcloud
+
+ self.r_wet = r_wet
+ self.mass_of_dry_air = mass_of_dry_air
+ self.n_output = 500
+
+ self.rtol_x = 0.5 * (condensation.DEFAULTS.rtol_x)
+ self.rtol_thd = condensation.DEFAULTS.rtol_thd
+ self.dt_cond_range = condensation.DEFAULTS.cond_range
diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py
new file mode 100644
index 0000000000..071dde8bf2
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py
@@ -0,0 +1,82 @@
+import numpy as np
+
+import PySDM.products as PySDM_products
+from PySDM.backends import CPU
+from PySDM.builder import Builder
+from PySDM.dynamics import AmbientThermodynamics, Condensation
+from PySDM.physics import constants as const
+from PySDM_examples.Loftus_and_Wordsworth_2021.parcel import AlienParcel
+
+
+class Simulation:
+ def __init__(self, settings, backend=CPU):
+ builder = Builder(
+ backend=backend(
+ formulae=settings.formulae,
+ **(
+ {"override_jit_flags": {"parallel": False}}
+ if backend == CPU
+ else {}
+ ),
+ ),
+ n_sd=1,
+ environment=AlienParcel(
+ dt=settings.dt,
+ mass_of_dry_air=settings.mass_of_dry_air,
+ pcloud=settings.pcloud,
+ zcloud=settings.Zcloud,
+ initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
+ Tcloud=settings.Tcloud,
+ ),
+ )
+
+ builder.add_dynamic(AmbientThermodynamics())
+ builder.add_dynamic(
+ Condensation(
+ rtol_x=settings.rtol_x,
+ rtol_thd=settings.rtol_thd,
+ dt_cond_range=settings.dt_cond_range,
+ )
+ )
+ builder.request_attribute("terminal velocity")
+
+ attributes = {}
+ r_dry = 1e-10
+ attributes["dry volume"] = settings.formulae.trivia.volume(radius=r_dry)
+ attributes["kappa times dry volume"] = attributes["dry volume"] * settings.kappa
+ attributes["multiplicity"] = np.array([1], dtype=np.int64)
+ r_wet = settings.r_wet
+ attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet)
+ products = [
+ PySDM_products.MeanRadius(name="radius_m1", unit="um"),
+ PySDM_products.ParcelDisplacement(name="z"),
+ PySDM_products.AmbientRelativeHumidity(name="RH", unit="%"),
+ PySDM_products.Time(name="t"),
+ ]
+
+ self.particulator = builder.build(attributes, products)
+
+ def save(self, output):
+ cell_id = 0
+ output["r"].append(
+ self.particulator.products["radius_m1"].get(unit=const.si.m)[cell_id]
+ )
+
+ output["z"].append(self.particulator.products["z"].get()[cell_id])
+ output["S"].append(self.particulator.products["RH"].get()[cell_id] / 100 - 1)
+ output["t"].append(self.particulator.products["t"].get())
+
+ def run(self):
+ output = {
+ "r": [],
+ "S": [],
+ "z": [],
+ "t": [],
+ }
+
+ self.save(output)
+ while self.particulator.environment["z"][0] > 0 and output["r"][-1] > 1e-6:
+ self.particulator.run(1)
+ self.save(output)
+
+ return output
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/__init__.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py
new file mode 100644
index 0000000000..75f57bb05d
--- /dev/null
+++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py
@@ -0,0 +1,231 @@
+import pytest
+import os
+import numpy as np
+from scipy.optimize import fsolve
+
+from PySDM import Formulae
+from PySDM.physics import si
+from PySDM_examples.Loftus_and_Wordsworth_2021 import Settings, Simulation
+from PySDM_examples.Loftus_and_Wordsworth_2021.planet import EarthLike
+
+
+class GroundTruthLoader:
+ def __init__(self, groundtruth_dir_path, n_samples=2, random_seed=2137):
+ self.dir_path = groundtruth_dir_path
+ self.RHs = None
+ self.r0grid = None
+ self.m_frac_evap = None
+ self.n_samples = n_samples # Number of random samples to test
+ np.random.seed(random_seed) # reproducible random samples during debugging
+
+ def __enter__(self):
+ try:
+ self.RHs = np.load(os.path.join(self.dir_path, "RHs.npy"))
+ self.r0grid = np.load(os.path.join(self.dir_path, "r0grid.npy"))
+ self.m_frac_evap = np.load(os.path.join(self.dir_path, "m_frac_evap.npy"))
+ return self
+ except FileNotFoundError as e:
+ print(f"Error loading ground truth files: {e}")
+ raise
+ except Exception as e:
+ print(f"An unexpected error occurred while loading ground truth data: {e}")
+ raise
+
+ def __exit__(self, exc_type, exc_val, exc_tb):
+ pass
+
+
+class TestNPYComparison:
+ @staticmethod
+ def _mix(dry_prop, vap_prop, ratio):
+ return (dry_prop + ratio * vap_prop) / (1 + ratio)
+
+ def _calculate_cloud_properties(
+ self, planet: EarthLike, surface_RH: float, formulae_instance: Formulae
+ ):
+ const = formulae_instance.constants
+
+ planet.RH_zref = surface_RH
+
+ pvs_stp = formulae_instance.saturation_vapour_pressure.pvs_water(planet.T_STP)
+
+ initial_water_vapour_mixing_ratio = const.eps / (
+ planet.p_STP / planet.RH_zref / pvs_stp - 1
+ )
+
+ R_air_mix = self._mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio)
+ cp_mix = self._mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio)
+
+ def solve_Tcloud(T_candidate):
+ pv_ad = (
+ initial_water_vapour_mixing_ratio
+ / (initial_water_vapour_mixing_ratio + const.eps)
+ * planet.p_STP
+ * (T_candidate / planet.T_STP) ** (cp_mix / R_air_mix)
+ )
+ pvs_tc = formulae_instance.saturation_vapour_pressure.pvs_water(T_candidate)
+ return pv_ad - pvs_tc
+
+ Tcloud_solutions = fsolve(solve_Tcloud, [150.0, 300.0])
+ Tcloud = np.max(Tcloud_solutions)
+
+ Zcloud = (planet.T_STP - Tcloud) * cp_mix / planet.g_std
+
+ th_std = formulae_instance.trivia.th_std(planet.p_STP, planet.T_STP)
+
+ pcloud = formulae_instance.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio(
+ planet.p_STP, th_std, initial_water_vapour_mixing_ratio, Zcloud
+ )
+ return initial_water_vapour_mixing_ratio, Tcloud, Zcloud, pcloud
+
+ def test_figure_2_replication_accuracy(self):
+ current_dir = os.path.dirname(os.path.abspath(__file__))
+ groundtruth_dir = os.path.abspath(os.path.join(current_dir, "ground_truth"))
+ if not os.path.isdir(groundtruth_dir):
+ pytest.fail(f"Groundtruth directory not found at {groundtruth_dir}")
+
+ formulae = Formulae(
+ ventilation="PruppacherAndRasmussen1979",
+ saturation_vapour_pressure="AugustRocheMagnus",
+ diffusion_coordinate="WaterMassLogarithm",
+ )
+
+ with GroundTruthLoader(groundtruth_dir) as gt:
+ if gt.RHs is None or gt.r0grid is None or gt.m_frac_evap is None:
+ pytest.fail("Ground truth data (.npy files) not loaded properly.")
+
+ n_rh_values = len(gt.RHs)
+ n_radius_values = gt.r0grid.shape[1]
+
+ total_points = n_rh_values * n_radius_values
+ n_samples = min(gt.n_samples, total_points)
+
+ if n_samples == 0:
+ pytest.skip("No data points available to sample.")
+
+ all_indices = np.array(
+ [(i, j) for i in range(n_rh_values) for j in range(n_radius_values)]
+ )
+
+ sampled_indices_flat = np.random.choice(len(all_indices), n_samples, replace=False)
+ sampled_ij_pairs = all_indices[sampled_indices_flat]
+
+ for i_rh, j_r in sampled_ij_pairs:
+ current_planet_state = EarthLike()
+
+ current_rh = gt.RHs[i_rh]
+ current_r_m = gt.r0grid[0, j_r]
+ expected_m_frac_evap = gt.m_frac_evap[i_rh, j_r]
+
+ try:
+ iwvmr, Tcloud, Zcloud, pcloud = self._calculate_cloud_properties(
+ current_planet_state, current_rh, formulae
+ )
+ simulated_m_frac_evap_point = self.calc_simulated_m_frac_evap_point(
+ current_planet_state,
+ formulae,
+ i_rh,
+ j_r,
+ current_rh,
+ current_r_m,
+ expected_m_frac_evap,
+ iwvmr,
+ Tcloud,
+ Zcloud,
+ pcloud,
+ )
+ except Exception as e:
+ print(
+ f"Warning: Error in _calculate_cloud_properties for RH={current_rh} (sample idx {i_rh},{j_r}): {e}."
+ )
+
+ error_context = (
+ f"Sample (RH_idx={i_rh}, R_idx={j_r}), "
+ f"RH={current_rh:.4f}, R_m={current_r_m:.3e}. "
+ f"Expected: {expected_m_frac_evap}, Got: {simulated_m_frac_evap_point}"
+ )
+
+ if np.isnan(expected_m_frac_evap):
+ assert np.isnan(
+ simulated_m_frac_evap_point
+ ), f"NaN Mismatch. {error_context} (Expected NaN, got non-NaN)"
+ else:
+ assert not np.isnan(
+ simulated_m_frac_evap_point
+ ), f"NaN Mismatch. {error_context} (Expected non-NaN, got NaN)"
+ np.testing.assert_allclose(
+ simulated_m_frac_evap_point,
+ expected_m_frac_evap,
+ rtol=1e-1, # Relative tolerance
+ atol=1e-1, # Absolute tolerance
+ err_msg=f"Value Mismatch. {error_context}",
+ )
+
+ def calc_simulated_m_frac_evap_point(
+ self,
+ current_planet_state,
+ formulae,
+ i_rh,
+ j_r,
+ current_rh,
+ current_r_m,
+ expected_m_frac_evap,
+ iwvmr,
+ Tcloud,
+ Zcloud,
+ pcloud,
+ ):
+
+ simulated_m_frac_evap_point = np.nan
+
+ if np.isnan(current_r_m) or current_r_m <= 0:
+ print(f"Warning: Invalid radius current_r_m={current_r_m} for sample idx {i_rh},{j_r}.")
+ else:
+ settings = Settings(
+ planet=current_planet_state,
+ r_wet=current_r_m,
+ mass_of_dry_air=1e5 * si.kg,
+ initial_water_vapour_mixing_ratio=iwvmr,
+ pcloud=pcloud,
+ Zcloud=Zcloud,
+ Tcloud=Tcloud,
+ formulae=formulae,
+ )
+ simulation = Simulation(settings)
+
+ try:
+ output = simulation.run()
+ if (
+ output
+ and "r" in output
+ and len(output["r"]) > 0
+ and "z" in output
+ and len(output["z"]) > 0
+ ):
+ final_radius_um = output["r"][-1]
+ if np.isnan(final_radius_um) or final_radius_um < 0:
+ final_radius_m = final_radius_um * 1e-6
+ if final_radius_m < 0: # Non-physical radius
+ simulated_m_frac_evap_point = 1.0 # 1.0 means fully evaporated
+ else:
+ simulated_m_frac_evap_point = np.nan
+ else:
+ final_radius_m = final_radius_um * 1e-6
+ if current_r_m == 0:
+ frac_evap = 1.0
+ else:
+ frac_evap = 1.0 - (final_radius_m / current_r_m) ** 3
+ frac_evap = np.clip(frac_evap, 0.0, 1.0)
+ simulated_m_frac_evap_point = frac_evap
+ else:
+ simulated_m_frac_evap_point = np.nan
+ except Exception as e:
+ print(
+ f"Warning: Simulation run failed for RH={current_rh}, r={current_r_m} (sample idx {i_rh},{j_r}): {e}."
+ )
+ if np.isclose(expected_m_frac_evap, 1.0, atol=1e-6):
+ simulated_m_frac_evap_point = 1.0
+ else:
+ simulated_m_frac_evap_point = np.nan
+
+ return simulated_m_frac_evap_point
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHgrid.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHgrid.npy
new file mode 100644
index 0000000000..e02761b427
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHgrid.npy differ
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHs.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHs.npy
new file mode 100644
index 0000000000..73132c6203
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHs.npy differ
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py
new file mode 100644
index 0000000000..ff8d780ac4
--- /dev/null
+++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py
@@ -0,0 +1,97 @@
+################################################################
+# make LoWo21 Figure 2
+# r_min, fraction raindrop mass evaporated as functions of RH
+################################################################
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+
+# load results
+root_dir = os.path.dirname(os.path.abspath(__file__))
+RHs = np.load(os.path.join(root_dir, "RHs.npy"))
+r0grid = np.load(os.path.join(root_dir, "r0grid.npy"))
+RHgrid = np.load(os.path.join(root_dir, "RHgrid.npy"))
+m_frac_evap = np.load(os.path.join(root_dir, "m_frac_evap.npy"))
+r_mins = np.load(os.path.join(root_dir, "r_mins.npy"))
+
+i_RH75 = 29 # index for RH=0.75
+
+# make figure 2
+# only put colorbar on lower panel, still line up x-axes
+f, axs = plt.subplots(
+ 2,
+ 2,
+ sharex="col",
+ figsize=(6, 7),
+ gridspec_kw={"height_ratios": [1, 3], "width_ratios": [20, 1]},
+)
+plt.subplots_adjust(hspace=0.05)
+axs[0, 0].set_xscale("log")
+axs[1, 0].set_xlabel(r"$r_0$ [mm]")
+axs[1, 0].set_ylabel("surface RH [ ]")
+axs[0, 0].tick_params(right=True, which="both")
+axs[1, 0].tick_params(right=True, which="both")
+axs[1, 0].tick_params(top=True, which="both")
+axs[0, 0].tick_params(top=True, which="both")
+axs[0, 0].set_ylabel("fraction mass \n evaporated [ ]")
+axs[0, 0].set_xlim(r0grid[0, 0] * 1e3, r0grid[0, -1] * 1e3)
+axs[0, 0].set_ylim(-0.04, 1.04)
+levels_smooth = np.linspace(0, 1, 250)
+cmesh = axs[1, 0].contourf(
+ r0grid * 1e3,
+ RHgrid,
+ m_frac_evap,
+ cmap=plt.cm.binary,
+ vmin=0,
+ vmax=1,
+ levels=levels_smooth,
+)
+
+
+cb = f.colorbar(cmesh, cax=axs[1, 1])
+cb.solids.set_edgecolor("face")
+axs[0, 1].axis("off")
+cb.solids.set_edgecolor("face")
+cb.solids.set_linewidth(1e-5)
+cb.set_label("fraction mass evaporated [ ]")
+cb.set_ticks([0, 0.1, 0.25, 0.5, 0.75, 1])
+axs[1, 0].axhline(0.75, lw=0.5, ls="--", c="plum")
+axs[1, 0].plot(r_mins * 1e3, RHs, lw=3, c="darkviolet", zorder=10)
+c_10 = axs[1, 0].contour(
+ r0grid * 1e3,
+ RHgrid,
+ m_frac_evap,
+ colors="indigo",
+ linewidths=1,
+ linestyles="--",
+ levels=[0.1],
+)
+cb.add_lines(c_10)
+axs[1, 0].clabel(c_10, c_10.levels, fmt={0.1: "10% mass evaporated"}, fontsize="smaller")
+axs[1, 0].fill_betweenx(RHs, 1e-2, (r_mins - 1e-6) * 1e3, edgecolor="k", facecolor="w", hatch="//")
+axs[1, 0].annotate(text="TOTAL \\n EVAPORATION", xy=(0.04, 0.55), c="k", backgroundcolor="w")
+axs[1, 0].annotate(
+ text=r"$r_\mathrm{min}$",
+ xy=(0.35, 0.27),
+ c="darkviolet",
+ backgroundcolor="w",
+ size=8,
+)
+
+axs[0, 0].scatter(r_mins[i_RH75] * 1e3, 0.99, color="darkviolet", zorder=10)
+axs[0, 0].axvline(r_mins[i_RH75] * 1e3, lw=0.5, c="darkviolet", ls="--")
+axs[0, 0].plot(r0grid[0, :] * 1e3, m_frac_evap[i_RH75, :], lw=2.05, c="k")
+axs[0, 0].axhline(1, c="w", lw=3)
+axs[0, 0].plot([1e-2, r_mins[i_RH75] * 1e3], [1, 1], c="k", lw=2.05, ls="--")
+axs[0, 0].annotate(text="surface RH = 0.75", xy=(0.8, 0.85), c="plum", size=8)
+axs[0, 0].annotate(text=r"$r_\mathrm{min}$", xy=(0.13, 0.05), c="darkviolet", size=8)
+
+figs_path = os.path.join(dir, "figs")
+os.mkdir(figs_path) if not os.path.exists(figs_path) else None
+plt.savefig(
+ os.path.join(figs_path, "fig02.pdf"),
+ transparent=True,
+ bbox_inches="tight",
+ pad_inches=0.5,
+)
+plt.close()
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/m_frac_evap.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/m_frac_evap.npy
new file mode 100644
index 0000000000..657765da4b
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/m_frac_evap.npy differ
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r0grid.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r0grid.npy
new file mode 100644
index 0000000000..5f4e75e31d
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r0grid.npy differ
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r_mins.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r_mins.npy
new file mode 100644
index 0000000000..dba178b0df
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r_mins.npy differ
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py
new file mode 100644
index 0000000000..d5b6425e39
--- /dev/null
+++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py
@@ -0,0 +1,238 @@
+from contextlib import contextmanager
+from PySDM import Formulae
+from PySDM.physics import si
+import pytest
+import numpy as np
+from scipy.optimize import fsolve
+
+from PySDM_examples.Loftus_and_Wordsworth_2021.planet import (
+ EarthLike,
+ Earth,
+ EarlyMars,
+ Jupiter,
+ Saturn,
+ K2_18B,
+)
+from PySDM_examples.Loftus_and_Wordsworth_2021.simulation import Simulation
+from PySDM_examples.Loftus_and_Wordsworth_2021.parcel import AlienParcel
+from PySDM_examples.Loftus_and_Wordsworth_2021 import Settings
+
+
+class TestLoftusWordsworth2021:
+
+ @contextmanager
+ def _get_test_resources(self):
+ formulae = Formulae(
+ ventilation="PruppacherAndRasmussen1979",
+ saturation_vapour_pressure="AugustRocheMagnus",
+ diffusion_coordinate="WaterMassLogarithm",
+ )
+ earth_like = EarthLike()
+ try:
+ yield formulae, earth_like
+ finally:
+ pass
+
+ def test_planet_classes(self):
+ """Test planet class instantiation and basic properties."""
+ planets = [EarthLike(), Earth(), EarlyMars(), Jupiter(), Saturn(), K2_18B()]
+
+ for planet in planets:
+ assert planet.g_std > 0
+ assert planet.T_STP > 0
+ assert planet.p_STP > 0
+ assert planet.RH_zref >= 0
+ assert planet.RH_zref <= 1
+
+ # atmospheric composition sums to 1 or less
+ total_conc = (
+ planet.dry_molar_conc_H2
+ + planet.dry_molar_conc_He
+ + planet.dry_molar_conc_N2
+ + planet.dry_molar_conc_O2
+ + planet.dry_molar_conc_CO2
+ )
+ assert (
+ total_conc <= 1.01
+ ), f"Total molar concentration {total_conc} exceeds 1.01 for {planet.__class__.__name__}"
+
+ def test_water_vapour_mixing_ratio_calculation(self):
+ """Test water vapour mixing ratio calculation."""
+ with self._get_test_resources() as (formulae, earth_like):
+ const = formulae.constants
+ planet = earth_like
+
+ pvs = formulae.saturation_vapour_pressure.pvs_water(planet.T_STP)
+ initial_water_vapour_mixing_ratio = const.eps / (
+ planet.p_STP / planet.RH_zref / pvs - 1
+ )
+
+ assert initial_water_vapour_mixing_ratio > 0
+ assert initial_water_vapour_mixing_ratio < 0.1 # Should be less than 10%
+
+ def test_alien_parcel_initialization(self):
+ """Test AlienParcel class initialization."""
+ parcel = AlienParcel(
+ dt=1.0 * si.second,
+ mass_of_dry_air=1e5 * si.kg,
+ pcloud=90000 * si.pascal,
+ initial_water_vapour_mixing_ratio=0.01,
+ Tcloud=280 * si.kelvin,
+ w=0,
+ zcloud=1000 * si.m,
+ )
+ assert parcel is not None
+
+ @pytest.mark.parametrize(
+ "r_wet_val, mass_of_dry_air_val, iwvmr_val, pcloud_val, Zcloud_val, Tcloud_val",
+ [
+ (1e-4, 1e5, 0.01, 90000, 1000, 280),
+ (1e-5, 1e4, 0.005, 80000, 500, 270),
+ (2e-4, 2e5, 0.02, 95000, 1500, 290),
+ ],
+ )
+ def test_simulation_class(
+ self,
+ r_wet_val,
+ mass_of_dry_air_val,
+ iwvmr_val,
+ pcloud_val,
+ Zcloud_val,
+ Tcloud_val,
+ ):
+ """Test Simulation class initialization and basic functionality with parametrized settings."""
+ with self._get_test_resources() as (formulae, earth_like):
+ planet = earth_like
+
+ settings = Settings(
+ planet=planet,
+ r_wet=r_wet_val * si.m,
+ mass_of_dry_air=mass_of_dry_air_val * si.kg,
+ initial_water_vapour_mixing_ratio=iwvmr_val,
+ pcloud=pcloud_val * si.pascal,
+ Zcloud=Zcloud_val * si.m,
+ Tcloud=Tcloud_val * si.kelvin,
+ formulae=formulae,
+ )
+
+ simulation = Simulation(settings)
+
+ assert hasattr(simulation, "particulator")
+ assert hasattr(simulation, "run")
+ assert hasattr(simulation, "save")
+
+ products = simulation.particulator.products
+ required_products = ["radius_m1", "z", "RH", "t"]
+ for product in required_products:
+ assert product in products
+
+ @pytest.mark.parametrize(
+ "r_wet_val, mass_of_dry_air_val, iwvmr_val, pcloud_val, Zcloud_val, Tcloud_val",
+ [
+ (1e-5, 1e5, 0.01, 90000, 100, 280),
+ (1e-5, 1e4, 0.005, 80000, 500, 270),
+ (2e-4, 2e5, 0.02, 95000, 1500, 290),
+ ],
+ )
+ def test_simulation_run_basic(
+ self,
+ r_wet_val,
+ mass_of_dry_air_val,
+ iwvmr_val,
+ pcloud_val,
+ Zcloud_val,
+ Tcloud_val,
+ ):
+ """Test basic simulation run functionality."""
+ with self._get_test_resources() as (formulae, earth_like):
+ planet = earth_like
+
+ settings = Settings(
+ planet=planet,
+ r_wet=r_wet_val * si.m,
+ mass_of_dry_air=mass_of_dry_air_val * si.kg,
+ initial_water_vapour_mixing_ratio=iwvmr_val,
+ pcloud=pcloud_val * si.pascal,
+ Zcloud=Zcloud_val * si.m,
+ Tcloud=Tcloud_val * si.kelvin,
+ formulae=formulae,
+ )
+
+ simulation = Simulation(settings)
+ output = simulation.run()
+
+ assert "r" in output
+ assert "S" in output
+ assert "z" in output
+ assert "t" in output
+
+ assert output["r"] is not None
+ assert output["S"] is not None
+ assert output["z"] is not None
+ assert output["t"] is not None
+
+ assert len(output["r"]) > 0, "Output array 'r' is empty"
+ assert len(output["S"]) > 0, "Output array 'S' is empty"
+ assert len(output["z"]) > 0, "Output array 'z' is empty"
+ assert len(output["t"]) > 0, "Output array 't' is empty"
+
+ lengths = [len(output[key]) for key in output.keys()]
+ assert all(
+ l == lengths[0] for l in lengths
+ ), "Not all output arrays have the same length"
+
+ def test_saturation_at_cloud_base(self):
+ formulae = Formulae(
+ ventilation="PruppacherAndRasmussen1979",
+ saturation_vapour_pressure="AugustRocheMagnus",
+ diffusion_coordinate="WaterMassLogarithm",
+ )
+
+ new_Earth = EarthLike()
+ new_Earth.T_STP
+
+ RH_array = np.linspace(0.25, 0.99, 50)
+ const = formulae.constants
+
+ def mix(dry, vap, ratio):
+ return (dry + ratio * vap) / (1 + ratio)
+
+ for i, RH in enumerate(RH_array[::-1]):
+ new_Earth.RH_zref = RH
+
+ pvs = formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP)
+ initial_water_vapour_mixing_ratio = const.eps / (
+ new_Earth.p_STP / new_Earth.RH_zref / pvs - 1
+ )
+
+ Rair = mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio)
+ c_p = mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio)
+
+ def f(x):
+ return initial_water_vapour_mixing_ratio / (
+ initial_water_vapour_mixing_ratio + const.eps
+ ) * new_Earth.p_STP * (x / new_Earth.T_STP) ** (
+ c_p / Rair
+ ) - formulae.saturation_vapour_pressure.pvs_water(
+ x
+ )
+
+ tdews = fsolve(f, [150, 300])
+ Tcloud = np.max(tdews)
+ Zcloud = (new_Earth.T_STP - Tcloud) * c_p / new_Earth.g_std
+ thstd = formulae.trivia.th_std(new_Earth.p_STP, new_Earth.T_STP)
+
+ pcloud = formulae.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio(
+ new_Earth.p_STP, thstd, initial_water_vapour_mixing_ratio, Zcloud
+ )
+
+ np.testing.assert_approx_equal(
+ actual=pcloud
+ * (
+ initial_water_vapour_mixing_ratio
+ / (initial_water_vapour_mixing_ratio + const.eps)
+ )
+ / formulae.saturation_vapour_pressure.pvs_water(Tcloud),
+ desired=1,
+ significant=4,
+ )
diff --git a/tests/examples_tests/conftest.py b/tests/examples_tests/conftest.py
index 74cafe3e12..581f8abdbe 100644
--- a/tests/examples_tests/conftest.py
+++ b/tests/examples_tests/conftest.py
@@ -83,6 +83,7 @@ def findfiles(path, regex):
"utils",
"Zaba_et_al_2025",
"Gonfiantini_1986",
+ "Loftus_and_Wordsworth_2021",
],
}