From 28434d7d5641c03e94bc82befd59e1b96a1baca4 Mon Sep 17 00:00:00 2001
From: Mazhar331 <84920914+Mazhar331@users.noreply.github.com>
Date: Fri, 17 Jul 2026 15:06:25 -0400
Subject: [PATCH 1/5] Fix drying latent heat factor (#24)

---
 PharmaPy/Drying_Model.py                | 40 +++++++++++-
 tests/test_drying_energy_rate_basis.py  |  7 +--
 tests/test_drying_latent_heat_factor.py | 83 +++++++++++++++++++++++++
 3 files changed, 125 insertions(+), 5 deletions(-)
 create mode 100644 tests/test_drying_latent_heat_factor.py

diff --git a/PharmaPy/Drying_Model.py b/PharmaPy/Drying_Model.py
index c5ce47a..610de87 100644
--- a/PharmaPy/Drying_Model.py
+++ b/PharmaPy/Drying_Model.py
@@ -424,6 +424,44 @@ def material_balance(self, time, satur, temp_gas, temp_sol, y_gas, x_liq,
 
     def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
                        u_gas, rho_gas, dry_rate, inputs, return_terms=False):
+        """Evaluate gas and condensed-phase energy balances.
+
+        Parameters
+        ----------
+        time : float
+            Current integration time [s].
+        temp_gas : ndarray
+            Gas-phase temperature by spatial node [K].
+        temp_sol : ndarray
+            Condensed-phase temperature by spatial node [K].
+        satur : ndarray
+            Cake saturation by spatial node [-].
+        y_gas : ndarray
+            Gas-phase mass fractions by node and species [-].
+        x_liq : ndarray
+            Liquid mass fractions by node and volatile species [-].
+        u_gas : ndarray
+            Gas velocity by spatial node [m/s].
+        rho_gas : ndarray
+            Gas density by spatial node [kg/m**3].
+        dry_rate : ndarray
+            Component drying rates on a mass basis [kg/m**3/s].
+        inputs : dict
+            Inlet condition dictionary; ``temp`` is the gas inlet
+            temperature [K].
+        return_terms : bool, optional
+            Return stored diagnostic terms instead of temperature derivatives [-].
+
+        Returns
+        -------
+        list of ndarray
+            ``dTg_dt`` and ``dTcond_dt`` by spatial node [K/s].
+
+        Notes
+        -----
+        ``latent_heat`` is requested on a mass basis [J/kg], so multiplying by
+        ``dry_rate`` gives the full latent power [J/m**3/s].
+        """
 
         # temp_ref = 298
         mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)
@@ -476,7 +514,7 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
         # heat_loss_cond = self.h_T_loss * self.a_V * (temp_sol - self.T_ambient)
         heat_loss_cond = self.h_T_loss * self.cake_height * (2*np.pi*1.5/2/100) *(temp_sol - self.T_ambient)
         heat_loss_cond = 0
-        drying_terms = (dry_rate[:, self.idx_volatiles] * latent_heat * 2).sum(axis=1)
+        drying_terms = (dry_rate[:, self.idx_volatiles] * latent_heat).sum(axis=1)
         denom_cond = self.rho_sol * (1 - self.porosity) * self.cp_sol + \
             self.porosity * satur * cpl_mix * dens_liq
 
diff --git a/tests/test_drying_energy_rate_basis.py b/tests/test_drying_energy_rate_basis.py
index 0245c72..37040c4 100644
--- a/tests/test_drying_energy_rate_basis.py
+++ b/tests/test_drying_energy_rate_basis.py
@@ -116,12 +116,11 @@ def test_unit_model_uses_mass_drying_rate_in_energy_terms(monkeypatch):
         2475.0 / 700.0,
     ])  # [K/s]
 
-    # The mass-rate latent powers are [256000, 338000, 128000] J/m**3/s. Issue #24
-    # separately tracks the existing factor of two applied to those powers.
+    # The mass-rate latent powers are [256000, 338000, 128000] J/m**3/s.
     expected_condensed_temperature_rate = np.array([
-        -512000.0 / 900000.0,
-        -676000.0 / 900000.0,
         -256000.0 / 900000.0,
+        -338000.0 / 900000.0,
+        -128000.0 / 900000.0,
     ])  # [K/s]
 
     np.testing.assert_allclose(
diff --git a/tests/test_drying_latent_heat_factor.py b/tests/test_drying_latent_heat_factor.py
new file mode 100644
index 0000000..6b28ade
--- /dev/null
+++ b/tests/test_drying_latent_heat_factor.py
@@ -0,0 +1,83 @@
+"""Regression coverage for Drying condensed-phase latent heat."""
+
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+
+pytest.importorskip("assimulo")
+import PharmaPy.Drying_Model as drying_module
+from PharmaPy.Drying_Model import Drying
+
+pytestmark = pytest.mark.assimulo
+
+
+def test_condensed_energy_uses_single_latent_heat_factor(monkeypatch):
+    """Mass drying rates multiplied by latent heat already give full power."""
+    dryer = Drying(number_nodes=3, supercrit_names=["nitrogen"])
+    dryer.idx_volatiles = np.array([0, 2])  # component indices [-]
+    dryer.idx_supercrit = np.array([1])  # component indices [-]
+    dryer.porosity = 0.5  # [-]
+    dryer.rho_sol = 1000.0  # [kg/m**3]
+    dryer.cp_sol = 1000.0  # [J/kg/K]
+    dryer.rho_liq = np.full(3, 800.0)  # [kg/m**3]
+    dryer.h_T_j = 0.0  # [W/m**2/K]
+    dryer.a_V = 1.0  # [m**2/m**3]
+    dryer.h_T_loss = 0.0  # [W/m**2/K]
+    dryer.cake_height = 1.0  # [m]
+    dryer.T_ambient = 298.15  # [K]
+    dryer.dz = np.ones(3)  # [m]
+
+    dryer.Vapor_1 = SimpleNamespace(
+        mw=np.array([28.0, 28.0, 28.0]),  # [g/mol]
+        getCp=lambda temp, mass_frac, basis: np.full(3, 1000.0),  # [J/kg/K]
+        getHeatVaporization=lambda temp, basis: np.array([2.0e6, 1.0e6]),  # [J/kg]
+    )
+    dryer.Liquid_1 = SimpleNamespace(
+        getCp=lambda temp, mass_frac, basis: np.full(3, 2000.0),  # [J/kg/K]
+    )
+
+    monkeypatch.setattr(
+        drying_module,
+        "high_resolution_fvm",
+        lambda values, boundary_cond: np.zeros(values.size + 1),  # [K]
+    )
+
+    temp_gas = np.full(3, 295.0)  # [K], matches energy_balance wet-bulb default
+    temp_sol = np.array([299.0, 304.0, 309.0])  # [K]
+    satur = np.full(3, 0.5)  # [-]
+    y_gas = np.array([
+        [0.10, 0.80, 0.10],
+        [0.20, 0.70, 0.10],
+        [0.30, 0.60, 0.10],
+    ])  # gas mass fractions [-]
+    x_liq = np.array([
+        [0.0, 1.0],
+        [0.0, 1.0],
+        [0.0, 1.0],
+    ])  # liquid volatile mass fractions [-]
+    dry_rate = np.array([
+        [0.036, 0.0, 0.184],
+        [0.054, 0.0, 0.230],
+        [0.018, 0.0, 0.092],
+    ])  # [kg/m**3/s]
+
+    _, dTcond_dt = dryer.energy_balance(
+        time=0.0,
+        temp_gas=temp_gas,
+        temp_sol=temp_sol,
+        satur=satur,
+        y_gas=y_gas,
+        x_liq=x_liq,
+        u_gas=np.zeros(3),  # [m/s]
+        rho_gas=np.ones(3),  # [kg/m**3]
+        dry_rate=dry_rate,
+        inputs={"temp": 295.0},  # [K]
+    )
+
+    expected_condensed_temperature_rate = np.array([
+        -256000.0 / 900000.0,
+        -338000.0 / 900000.0,
+        -128000.0 / 900000.0,
+    ])  # [K/s]
+    np.testing.assert_allclose(dTcond_dt, expected_condensed_temperature_rate)

From 790cc52c03fe4f2f05595bf38f1f27c10b32b319 Mon Sep 17 00:00:00 2001
From: Mazhar331 <84920914+Mazhar331@users.noreply.github.com>
Date: Fri, 17 Jul 2026 15:29:33 -0400
Subject: [PATCH 2/5] Document drying energy diagnostic returns

---
 PharmaPy/Drying_Model.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/PharmaPy/Drying_Model.py b/PharmaPy/Drying_Model.py
index 610de87..4ab8755 100644
--- a/PharmaPy/Drying_Model.py
+++ b/PharmaPy/Drying_Model.py
@@ -455,7 +455,11 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
         Returns
         -------
         list of ndarray
-            ``dTg_dt`` and ``dTcond_dt`` by spatial node [K/s].
+            ``dTg_dt`` and ``dTcond_dt`` by spatial node [K/s] when
+            ``return_terms`` is False.
+        tuple of ndarray
+            When ``return_terms`` is True, returns the diagnostic terms
+            ``(convec_term, drying, heat_cond, heat_loss_emp)`` instead.
 
         Notes
         -----

From 733f149d993821126738d3871b31e4c7c142a9d9 Mon Sep 17 00:00:00 2001
From: Mazhar331 <84920914+Mazhar331@users.noreply.github.com>
Date: Tue, 21 Jul 2026 13:03:52 -0400
Subject: [PATCH 3/5] Clarify drying energy diagnostics

---
 PharmaPy/Drying_Model.py               | 77 +++++++++++++-------------
 tests/test_drying_energy_rate_basis.py |  4 +-
 2 files changed, 40 insertions(+), 41 deletions(-)

diff --git a/PharmaPy/Drying_Model.py b/PharmaPy/Drying_Model.py
index 4ab8755..541c778 100644
--- a/PharmaPy/Drying_Model.py
+++ b/PharmaPy/Drying_Model.py
@@ -69,9 +69,9 @@ def __init__(self, number_nodes, supercrit_names, diam_unit=0.01,
         self.T_ambient = 298
 
         # Transfer coefficients
-        self.k_y = 1e-3  # mol/s/m**2 (Seader, Separation process)
-        # self.h_T_j = 30  # W/m**2/K
-        self.h_T_j = 10  # W/m**2/K
+        self.k_y = 1e-3  # [mol/s/m**2] (Seader, Separation process)
+        # self.h_T_j = 30  # [W/m**2/K]
+        self.h_T_j = 10  # [W/m**2/K]
         self.h_T_loss = 30
 
         self._Phases = None
@@ -291,13 +291,16 @@ def unit_model(self, time, states, sw=None):
         sat_red = (satur - self.s_inf) / (1 - self.s_inf)
         sat_red = np.maximum(0, sat_red)
         k_ra = (1 - sat_red)**2 * (1 - sat_red**1.4)
-        # vel_gas = self.k_perm * k_ra * self.dPg_dz / visc_gas
-        vel_gas = self.dPg_dz/ \
-        (self.CakePhase.alpha * visc_gas * self.rho_sol * (1 - self.porosity))/ np.mean(satur)
+        vel_gas = (
+            self.dPg_dz /
+            (self.CakePhase.alpha * visc_gas * self.rho_sol *
+             (1 - self.porosity)) /
+            np.mean(satur)
+        )  # [m/s]
         
         # ---------- Drying rate term
         mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)
-        rho_gas = self.pres_gas / gas_ct / temp_gas * mw_avg_gas / 1000# kg/m**3
+        rho_gas = self.pres_gas / gas_ct / temp_gas * mw_avg_gas / 1000  # [kg/m**3]
         rho_liq_ = self.Liquid_1.rho_liq[self.idx_volatiles]
         self.rho_liq =  1 / np.sum((x_liq/ rho_liq_), axis=1)
         # Dry correction
@@ -333,9 +336,6 @@ def unit_model(self, time, states, sw=None):
                                           satur, y_gas, x_liq, vel_gas,
                                           rho_gas, self.dry_rate, inputs)
 
-        # print(satur.min())
-        # print(inputs.values())
-
         model_eqns = np.column_stack(material_eqns + energy_eqns)
 
         self.derivatives = model_eqns.ravel()
@@ -374,8 +374,10 @@ def material_balance(self, time, satur, temp_gas, temp_sol, y_gas, x_liq,
 
         Returns
         -------
-        list of ndarray
+        list of ndarray or int
             ``dsat_dt`` [1/s], ``dygas_dt`` [1/s], and ``dxliq_dt`` [1/s].
+            When ``return_terms`` is True, returns the legacy
+            ``masstrans_comp`` diagnostic flag [-].
         """
         
         satur[satur < eps] = eps
@@ -414,7 +416,7 @@ def material_balance(self, time, satur, temp_gas, temp_sol, y_gas, x_liq,
 
         dygas_dt = convection + transfer_gas + total_mass_correction  # [1/s]
 
-        if return_terms:    # TODO: check term by term in material balance down this line
+        if return_terms:
             self.masstrans_comp = 1
 
             return self.masstrans_comp
@@ -460,6 +462,9 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
         tuple of ndarray
             When ``return_terms`` is True, returns the diagnostic terms
             ``(convec_term, drying, heat_cond, heat_loss_emp)`` instead.
+            ``convec_term`` is the raw ``u_gas * dTg_dz`` diagnostic
+            [kg*K/m**3/s]; the remaining entries are gas-temperature-rate
+            contributions [K/s].
 
         Notes
         -----
@@ -475,27 +480,28 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
         # ----- Gas phase equations
         cpg_mix = self.Vapor_1.getCp(temp=temp_gas, mass_frac=y_gas,
                                      basis='mass')
-        cvg_mix = cpg_mix - gas_ct / mw_avg_gas * 1000  # J/kg K
+        cvg_mix = cpg_mix - gas_ct / mw_avg_gas * 1000  # [J/kg/K]
 
-        epsilon_gas = self.porosity * (1 - satur)
-        denom_gas = cvg_mix * epsilon_gas * rho_gas
+        epsilon_gas = self.porosity * (1 - satur)  # [-]
+        denom_gas = cvg_mix * epsilon_gas * rho_gas  # [J/m**3/K]
 
         latent_heat = self.Vapor_1.getHeatVaporization(temp_sol,
-                                                       basis='mass')
+                                                       basis='mass')  # [J/kg]
 
         xliq_extended = np.column_stack((x_liq, np.zeros((x_liq.shape[0],
                                                           len(self.idx_supercrit)))))
 
         cpl_mix = self.Liquid_1.getCp(temp=temp_sol, mass_frac=xliq_extended,
                                       basis='mass')
-        temp_wb = 22+273
-        sensible_heat = cpg_mix * (temp_gas - temp_wb) * dry_rate.sum(axis=1)
+        temp_wb = 22 + 273  # [K]
+        sensible_heat = (
+            cpg_mix * (temp_gas - temp_wb) * dry_rate.sum(axis=1)
+        )  # [J/m**3/s]
         
-        heat_transf = self.h_T_j * self.a_V * (temp_gas - temp_sol)
-        drying_terms = rho_gas / self.rho_liq * cpg_mix
-        # heat_loss = self.h_T_loss * self.a_V * (temp_gas - self.T_ambient)
-        heat_loss = self.h_T_loss * self.cake_height * (2*np.pi*1.5/2/100) *(temp_gas - self.T_ambient)
-        heat_loss = 0  # This line is for assumption of no heat loss
+        heat_transf = (
+            self.h_T_j * self.a_V * (temp_gas - temp_sol)
+        )  # [J/m**3/s]
+        heat_loss = np.zeros_like(temp_gas)  # [J/m**3/s]
         fluxes_Tg = high_resolution_fvm(temp_gas,
                                         boundary_cond=temp_gas_inputs)
         
@@ -511,24 +517,17 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
         # dTg_dt = -u_gas * dTg_dz + (-heat_loss) / denom_gas
         # dTg_dt =  (conv_term - heat_loss) / denom_gas
 
-        # print(dTg_dt[0])
-
         # ----- Condensed phases equations
-        dens_liq = self.rho_liq
-        # heat_loss_cond = self.h_T_loss * self.a_V * (temp_sol - self.T_ambient)
-        heat_loss_cond = self.h_T_loss * self.cake_height * (2*np.pi*1.5/2/100) *(temp_sol - self.T_ambient)
-        heat_loss_cond = 0
-        drying_terms = (dry_rate[:, self.idx_volatiles] * latent_heat).sum(axis=1)
+        dens_liq = self.rho_liq  # [kg/m**3]
+        heat_loss_cond = np.zeros_like(temp_sol)  # [J/m**3/s]
+        drying_terms = (
+            dry_rate[:, self.idx_volatiles] * latent_heat
+        ).sum(axis=1)  # [J/m**3/s]
         denom_cond = self.rho_sol * (1 - self.porosity) * self.cp_sol + \
             self.porosity * satur * cpl_mix * dens_liq
 
         dTcond_dt = (-drying_terms + heat_transf - heat_loss_cond) / denom_cond
         
-        ## ---- Trial for lumping both cond/gas phase into one
-        # dTtotal_dt = (conv_term + sensible_heat - drying_terms - heat_loss_cond)/ (denom_gas + denom_cond)
-        
-        # dTg_dt, dTcond_dt = dTtotal_dt, dTtotal_dt
-        
         if return_terms:
             self.convec_term = u_gas * dTg_dz
             self.drying = drying_terms/ denom_gas
@@ -644,9 +643,9 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True,
 
         # Mass transfer
         moments = self.Solid_1.getMoments(mom_num=[0, 1, 2, 3, 4])
-        sauter_diam = moments[1] / moments[0]  # m
+        sauter_diam = moments[1] / moments[0]  # [m]
 
-        # self.a_V = 6 / sauter_diam  # m**2/m**3
+        # self.a_V = 6 / sauter_diam  # [m**2/m**3]
         self.a_V = moments[2] * (1 - porosity) / moments[3]
         # Gas pressure
         # deltaP_media = deltaP*self.resist_medium / \
@@ -668,10 +667,10 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True,
         csd = self.Solid_1.distrib# * 1e6
         mom_zero = self.Solid_1.moments[0]
 
-        # rholiq_mass = np.mean(rho_liq[0] * self.Liquid_1.mw_av[0])  # kg/m**3
+        # rholiq_mass = np.mean(rho_liq[0] * self.Liquid_1.mw_av[0])  # [kg/m**3]
         self.s_inf = get_sat_inf(x_csd, csd, deltaP, porosity,
                                  self.cake_height, mom_zero,
-                                 (np.mean(surf_tens), rho_liq[0]))#rholiq_mass))
+                                 (np.mean(surf_tens), rho_liq[0]))
 
         # ---------- Solve model
         model = self.unit_model
diff --git a/tests/test_drying_energy_rate_basis.py b/tests/test_drying_energy_rate_basis.py
index 37040c4..5214013 100644
--- a/tests/test_drying_energy_rate_basis.py
+++ b/tests/test_drying_energy_rate_basis.py
@@ -109,14 +109,14 @@ def test_unit_model_uses_mass_drying_rate_in_energy_terms(monkeypatch):
     derivatives = derivatives.reshape(3, -1)
 
     # Sensible power [J/m**3/s] divided by gas capacity [J/m**3/K].
-    # power = cp_gas [J/kg/K] * (T_gas - 295 K) * sum(m_dot_i) [kg/m**3/s].
+    # power = cp_gas [J/kg/K] * (T_gas - 295 [K]) * sum(m_dot_i) [kg/m**3/s].
     expected_gas_temperature_rate = np.array([
         1100.0 / 450.0,
         3408.0 / 1100.0,
         2475.0 / 700.0,
     ])  # [K/s]
 
-    # The mass-rate latent powers are [256000, 338000, 128000] J/m**3/s.
+    # The mass-rate latent powers are [256000, 338000, 128000] [J/m**3/s].
     expected_condensed_temperature_rate = np.array([
         -256000.0 / 900000.0,
         -338000.0 / 900000.0,

From 55a530be6826390738e7b2ffd2f292645e52a14b Mon Sep 17 00:00:00 2001
From: Mazhar331 <84920914+Mazhar331@users.noreply.github.com>
Date: Tue, 21 Jul 2026 16:52:31 -0400
Subject: [PATCH 4/5] Annotate drying energy units

---
 PharmaPy/Drying_Model.py | 130 +++++++++++++++++++++------------------
 1 file changed, 70 insertions(+), 60 deletions(-)

diff --git a/PharmaPy/Drying_Model.py b/PharmaPy/Drying_Model.py
index 541c778..935b42d 100644
--- a/PharmaPy/Drying_Model.py
+++ b/PharmaPy/Drying_Model.py
@@ -139,6 +139,19 @@ def Inlet(self, inlet):
         self._Inlet = inlet
 
     def nomenclature(self):
+        """Create drying-model state metadata.
+
+        Returns
+        -------
+        None
+            The method populates inlet, outlet, state, and final-state metadata
+            on the instance.
+
+        Notes
+        -----
+        Saturation and mass fractions are dimensionless [-]. Gas and condensed
+        temperatures are tracked in [K].
+        """
         self.names_states_in = ['temp', 'mass_frac']
         self.names_states_out = self.names_states_in
 
@@ -152,17 +165,17 @@ def nomenclature(self):
         self.states_di = {
             'saturation': {# 'index': index_z,
                            'dim': 1,
-                           'units': '', 'type': 'diff'},
+                           'units': '[-]', 'type': 'diff'},
             'y_gas': {'index': self.name_species,
-                      'dim': len(self.name_species), 'units': '',
+                      'dim': len(self.name_species), 'units': '[-]',
                       'type': 'diff'},
-            'x_liq': {'index': name_liq, 'dim': len(name_liq), 'units': '',
+            'x_liq': {'index': name_liq, 'dim': len(name_liq), 'units': '[-]',
                       'type': 'diff'},
             'temp_gas': {# 'index': index_z,
-                         'dim': 1, 'units': 'K',
+                         'dim': 1, 'units': '[K]',
                          'type': 'diff'},
             'temp_cond': {# 'index': index_z,
-                         'dim': 1, 'units': 'K',
+                         'dim': 1, 'units': '[K]',
                          'type': 'diff'},
             }
 
@@ -273,49 +286,49 @@ def unit_model(self, time, states, sw=None):
             are [K/s].
         """
 
-        num_comp = self.Liquid_1.num_species
-        states_reord = states.reshape(-1, 3 + num_comp + self.num_volatiles)
+        num_comp = self.Liquid_1.num_species  # [-]
+        states_reord = states.reshape(
+            -1, 3 + num_comp + self.num_volatiles)  # [-] and [K]
 
-        satur = states_reord[:, 0]
-        y_gas = states_reord[:, 1:1 + num_comp]
+        satur = states_reord[:, 0]  # [-]
+        y_gas = states_reord[:, 1:1 + num_comp]  # [-]
         x_liq = states_reord[:, 1 + num_comp:
-                             1 + num_comp + self.num_volatiles]
-        x_liq[:, -2] = 0
-        temp_gas = states_reord[:, -2]
-        temp_sol = states_reord[:, -1]
+                             1 + num_comp + self.num_volatiles]  # [-]
+        x_liq[:, -2] = 0  # [-]
+        temp_gas = states_reord[:, -2]  # [K]
+        temp_sol = states_reord[:, -1]  # [K]
 
         # ---------- Darcy's equation
         visc_gas = self.Vapor_1.getViscosity(temp=temp_gas,
-                                             mass_frac=y_gas)
+                                             mass_frac=y_gas)  # [Pa*s]
 
-        sat_red = (satur - self.s_inf) / (1 - self.s_inf)
-        sat_red = np.maximum(0, sat_red)
-        k_ra = (1 - sat_red)**2 * (1 - sat_red**1.4)
+        sat_red = (satur - self.s_inf) / (1 - self.s_inf)  # [-]
+        sat_red = np.maximum(0, sat_red)  # [-]
+        k_ra = (1 - sat_red)**2 * (1 - sat_red**1.4)  # [-]
         vel_gas = (
             self.dPg_dz /
             (self.CakePhase.alpha * visc_gas * self.rho_sol *
              (1 - self.porosity)) /
             np.mean(satur)
         )  # [m/s]
-        
+
         # ---------- Drying rate term
-        mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)
+        mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)  # [g/mol]
         rho_gas = self.pres_gas / gas_ct / temp_gas * mw_avg_gas / 1000  # [kg/m**3]
-        rho_liq_ = self.Liquid_1.rho_liq[self.idx_volatiles]
-        self.rho_liq =  1 / np.sum((x_liq/ rho_liq_), axis=1)
+        rho_liq_ = self.Liquid_1.rho_liq[self.idx_volatiles]  # [kg/m**3]
+        self.rho_liq =  1 / np.sum((x_liq/ rho_liq_), axis=1)  # [kg/m**3]
         # Dry correction
         if self.mass_eta:
-            rho_liq = self.rho_liq
-            
+            rho_liq = self.rho_liq  # [kg/m**3]
+
             sat_eta = (self.porosity * satur * rho_liq)/ \
-                ((1 - self.porosity) * self.rho_sol + self.porosity * satur *rho_liq)
-            # sat_eta = satur * rho_liq / (satur*rho_liq + (1 - satur)*rho_gas)
-            w_eta = x_liq
+                ((1 - self.porosity) * self.rho_sol + self.porosity * satur *rho_liq)  # [-]
+            w_eta = x_liq  # [-]
         else:
-            sat_eta = satur
-            w_eta = x_liq
+            sat_eta = satur  # [-]
+            w_eta = x_liq  # [-]
 
-        limiter_factor = self.eta_fun(sat_eta, w_eta)
+        limiter_factor = self.eta_fun(sat_eta, w_eta)  # [-]
 
         # Drying rate from get_drying_rate is molar [mol/m**3/s].
         self.dry_rate = self.get_drying_rate(x_liq, temp_sol, y_gas,
@@ -323,22 +336,22 @@ def unit_model(self, time, states, sw=None):
 
         # Balances below consume mass-basis drying rates [kg/m**3/s].
         self.dry_rate = self._drying_rate_mass_basis(self.dry_rate)
-        self.dry_rate *= limiter_factor[..., np.newaxis]
+        self.dry_rate *= limiter_factor[..., np.newaxis]  # [kg/m**3/s]
 
         # ---------- Model equations
-        inputs = self.get_inputs(time)['Inlet']
+        inputs = self.get_inputs(time)['Inlet']  # inlet states: [-] and [K]
 
         material_eqns = self.material_balance(
             time, satur, temp_gas, temp_sol, y_gas, x_liq,
-            vel_gas, rho_gas, self.dry_rate, inputs)
+            vel_gas, rho_gas, self.dry_rate, inputs)  # [1/s]
 
         energy_eqns = self.energy_balance(time, temp_gas, temp_sol,
                                           satur, y_gas, x_liq, vel_gas,
-                                          rho_gas, self.dry_rate, inputs)
+                                          rho_gas, self.dry_rate, inputs)  # [K/s]
 
-        model_eqns = np.column_stack(material_eqns + energy_eqns)
+        model_eqns = np.column_stack(material_eqns + energy_eqns)  # [1/s, K/s]
 
-        self.derivatives = model_eqns.ravel()
+        self.derivatives = model_eqns.ravel()  # [1/s, K/s]
 
         return model_eqns.ravel()
 
@@ -469,17 +482,20 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
         Notes
         -----
         ``latent_heat`` is requested on a mass basis [J/kg], so multiplying by
-        ``dry_rate`` gives the full latent power [J/m**3/s].
+        ``dry_rate`` gives the full latent power [J/m**3/s]. The existing
+        gas-convection discretization is preserved in this branch; ``dTg_dz``
+        [kg*K/m**4] and ``conv_term`` [J*kg/m**6/s] are annotated as
+        implemented so that the remaining dimensional debt is explicit without
+        expanding the #24 behavior change.
         """
 
-        # temp_ref = 298
-        mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)
+        mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)  # [g/mol]
         # ----- Reading inputs
-        temp_gas_inputs = inputs['temp']
+        temp_gas_inputs = inputs['temp']  # [K]
 
         # ----- Gas phase equations
         cpg_mix = self.Vapor_1.getCp(temp=temp_gas, mass_frac=y_gas,
-                                     basis='mass')
+                                     basis='mass')  # [J/kg/K]
         cvg_mix = cpg_mix - gas_ct / mw_avg_gas * 1000  # [J/kg/K]
 
         epsilon_gas = self.porosity * (1 - satur)  # [-]
@@ -489,10 +505,10 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
                                                        basis='mass')  # [J/kg]
 
         xliq_extended = np.column_stack((x_liq, np.zeros((x_liq.shape[0],
-                                                          len(self.idx_supercrit)))))
+                                                          len(self.idx_supercrit)))))  # [-]
 
         cpl_mix = self.Liquid_1.getCp(temp=temp_sol, mass_frac=xliq_extended,
-                                      basis='mass')
+                                      basis='mass')  # [J/kg/K]
         temp_wb = 22 + 273  # [K]
         sensible_heat = (
             cpg_mix * (temp_gas - temp_wb) * dry_rate.sum(axis=1)
@@ -503,19 +519,13 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
         )  # [J/m**3/s]
         heat_loss = np.zeros_like(temp_gas)  # [J/m**3/s]
         fluxes_Tg = high_resolution_fvm(temp_gas,
-                                        boundary_cond=temp_gas_inputs)
-        
-        # fluxes_Tg = upwind_fvm(temp_gas, boundary_cond=temp_gas_inputs)
-        # sensible_heat = cpg_mix * np.diff(fluxes_Tg) * dry_rate.sum(axis=1)
-        dTg_dz = np.diff(fluxes_Tg) / self.dz * epsilon_gas * rho_gas
+                                        boundary_cond=temp_gas_inputs)  # [K]
 
-        conv_term = -u_gas * dTg_dz * cpg_mix * rho_gas
+        dTg_dz = np.diff(fluxes_Tg) / self.dz * epsilon_gas * rho_gas  # [kg*K/m**4]
 
-        dTg_dt = (conv_term + sensible_heat - heat_transf - heat_loss) / denom_gas
-        
-        # Empty port
-        # dTg_dt = -u_gas * dTg_dz + (-heat_loss) / denom_gas
-        # dTg_dt =  (conv_term - heat_loss) / denom_gas
+        conv_term = -u_gas * dTg_dz * cpg_mix * rho_gas  # [J*kg/m**6/s]
+
+        dTg_dt = (conv_term + sensible_heat - heat_transf - heat_loss) / denom_gas  # [K/s]
 
         # ----- Condensed phases equations
         dens_liq = self.rho_liq  # [kg/m**3]
@@ -524,15 +534,15 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
             dry_rate[:, self.idx_volatiles] * latent_heat
         ).sum(axis=1)  # [J/m**3/s]
         denom_cond = self.rho_sol * (1 - self.porosity) * self.cp_sol + \
-            self.porosity * satur * cpl_mix * dens_liq
+            self.porosity * satur * cpl_mix * dens_liq  # [J/m**3/K]
+
+        dTcond_dt = (-drying_terms + heat_transf - heat_loss_cond) / denom_cond  # [K/s]
 
-        dTcond_dt = (-drying_terms + heat_transf - heat_loss_cond) / denom_cond
-        
         if return_terms:
-            self.convec_term = u_gas * dTg_dz
-            self.drying = drying_terms/ denom_gas
-            self.heat_cond = heat_transf/ denom_gas
-            self.heat_loss_emp = heat_loss/ denom_gas
+            self.convec_term = u_gas * dTg_dz  # [kg*K/m**3/s]
+            self.drying = drying_terms/ denom_gas  # [K/s]
+            self.heat_cond = heat_transf/ denom_gas  # [K/s]
+            self.heat_loss_emp = heat_loss/ denom_gas  # [K/s]
 
             return self.convec_term, self.drying, self.heat_cond, self.heat_loss_emp
 

From 7c12c1dd8f47c21d7cf4ae0e1a6cbafad384287f Mon Sep 17 00:00:00 2001
From: Mazhar331 <84920914+Mazhar331@users.noreply.github.com>
Date: Tue, 21 Jul 2026 18:12:09 -0400
Subject: [PATCH 5/5] Address drying latent heat audit follow-up

---
 PharmaPy/Drying_Model.py                | 194 ++++++++++++++++++------
 PharmaPy/Plotting.py                    |  40 ++++-
 tests/test_drying_latent_heat_factor.py |  39 +++--
 tests/test_plotting_unit_labels.py      |  35 +++++
 4 files changed, 243 insertions(+), 65 deletions(-)
 create mode 100644 tests/test_plotting_unit_labels.py

diff --git a/PharmaPy/Drying_Model.py b/PharmaPy/Drying_Model.py
index 935b42d..c0d4a2e 100644
--- a/PharmaPy/Drying_Model.py
+++ b/PharmaPy/Drying_Model.py
@@ -24,9 +24,9 @@
 from PharmaPy.Results import DynamicResult
 from PharmaPy.Plotting import plot_distrib
 
-eps = np.finfo(float).eps
+eps = np.finfo(float).eps  # [-]
 
-gas_ct = 8.314
+gas_ct = 8.314  # [J/mol/K]
 
 
 class Drying:
@@ -62,17 +62,17 @@ def __init__(self, number_nodes, supercrit_names, diam_unit=0.01,
         """
         self.supercrit_names = supercrit_names
 
-        self.num_nodes = number_nodes
-        self.station_diameter = diam_unit
-        self.area_cross = self.station_diameter**2 * np.pi/4
-        self.resist_medium = resist_medium
-        self.T_ambient = 298
+        self.num_nodes = number_nodes  # [-]
+        self.station_diameter = diam_unit  # [m]
+        self.area_cross = self.station_diameter**2 * np.pi/4  # [m**2]
+        self.resist_medium = resist_medium  # [1/m]
+        self.T_ambient = 298  # [K]
 
         # Transfer coefficients
         self.k_y = 1e-3  # [mol/s/m**2] (Seader, Separation process)
         # self.h_T_j = 30  # [W/m**2/K]
         self.h_T_j = 10  # [W/m**2/K]
-        self.h_T_loss = 30
+        self.h_T_loss = 30  # [W/m**2/K]
 
         self._Phases = None
         self._Inlet = None
@@ -88,7 +88,7 @@ def __init__(self, number_nodes, supercrit_names, diam_unit=0.01,
         self.state_event_list = state_events
 
         self.outputs = None
-        self.elapsed_time = 0
+        self.elapsed_time = 0  # [s]
 
     @property
     def Phases(self):
@@ -199,17 +199,33 @@ def _eval_state_events(self, time, states, sw):
         return events
 
     def get_y_equilib(self, temp_cond, x_liq, p_gas):
-        mw_liq = self.Liquid_1.mw[self.idx_volatiles]
+        """Calculate equilibrium gas mole fractions above the liquid.
+
+        Parameters
+        ----------
+        temp_cond : ndarray
+            Condensed-phase temperature by spatial node [K].
+        x_liq : ndarray
+            Liquid mass fractions for volatile components [-].
+        p_gas : ndarray or float
+            Gas pressure [Pa].
+
+        Returns
+        -------
+        ndarray
+            Equilibrium gas mole fractions for volatile components [-].
+        """
+        mw_liq = self.Liquid_1.mw[self.idx_volatiles]  # [g/mol]
         x_liq_mole_frac = (x_liq / mw_liq).T / np.dot(1/mw_liq, x_liq.T)
-        x_liq_mole_frac = x_liq_mole_frac.T
+        x_liq_mole_frac = x_liq_mole_frac.T  # [-]
 
-        p_sat = self.Liquid_1.AntoineEquation(temp=temp_cond)
+        p_sat = self.Liquid_1.AntoineEquation(temp=temp_cond)  # [Pa]
 
-        gamma = self.Liquid_1.getActivityCoeff(mole_frac=x_liq_mole_frac)
+        gamma = self.Liquid_1.getActivityCoeff(mole_frac=x_liq_mole_frac)  # [-]
         # p_gas_total = np.sum(x_liq_mole_frac * p_sat[:, self.idx_volatiles],
         #                      axis=1)
-        p_partial = (gamma * x_liq_mole_frac * p_sat[:, self.idx_volatiles]).T
-        y_equil = p_partial  / p_gas
+        p_partial = (gamma * x_liq_mole_frac * p_sat[:, self.idx_volatiles]).T  # [Pa]
+        y_equil = p_partial  / p_gas  # [-]
 
         return y_equil
 
@@ -284,6 +300,14 @@ def unit_model(self, time, states, sw=None):
             Flattened state derivatives ordered like ``states``. Saturation
             and mass-fraction derivatives are [1/s]; temperature derivatives
             are [K/s].
+
+        Notes
+        -----
+        The gas-density path still uses the legacy mass-fraction weighted
+        molecular-weight surrogate [g/mol]; issue #28 owns replacing it with
+        the true mixture molecular weight. The ``x_liq`` supercritical slot
+        reset preserves the existing state layout; issue #42 owns replacing
+        the magic index with explicit metadata.
         """
 
         num_comp = self.Liquid_1.num_species  # [-]
@@ -294,7 +318,7 @@ def unit_model(self, time, states, sw=None):
         y_gas = states_reord[:, 1:1 + num_comp]  # [-]
         x_liq = states_reord[:, 1 + num_comp:
                              1 + num_comp + self.num_volatiles]  # [-]
-        x_liq[:, -2] = 0  # [-]
+        x_liq[:, -2] = 0  # legacy supercritical component slot [-]
         temp_gas = states_reord[:, -2]  # [K]
         temp_sol = states_reord[:, -1]  # [K]
 
@@ -302,9 +326,6 @@ def unit_model(self, time, states, sw=None):
         visc_gas = self.Vapor_1.getViscosity(temp=temp_gas,
                                              mass_frac=y_gas)  # [Pa*s]
 
-        sat_red = (satur - self.s_inf) / (1 - self.s_inf)  # [-]
-        sat_red = np.maximum(0, sat_red)  # [-]
-        k_ra = (1 - sat_red)**2 * (1 - sat_red**1.4)  # [-]
         vel_gas = (
             self.dPg_dz /
             (self.CakePhase.alpha * visc_gas * self.rho_sol *
@@ -313,7 +334,7 @@ def unit_model(self, time, states, sw=None):
         )  # [m/s]
 
         # ---------- Drying rate term
-        mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)  # [g/mol]
+        mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)  # legacy MW surrogate [g/mol]
         rho_gas = self.pres_gas / gas_ct / temp_gas * mw_avg_gas / 1000  # [kg/m**3]
         rho_liq_ = self.Liquid_1.rho_liq[self.idx_volatiles]  # [kg/m**3]
         self.rho_liq =  1 / np.sum((x_liq/ rho_liq_), axis=1)  # [kg/m**3]
@@ -470,14 +491,16 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
         Returns
         -------
         list of ndarray
-            ``dTg_dt`` and ``dTcond_dt`` by spatial node [K/s] when
-            ``return_terms`` is False.
+            ``dTcond_dt`` by spatial node [K/s] and the legacy gas-temperature
+            solver channel ``dTg_dt`` when ``return_terms`` is False. Issue #37
+            owns the gas-convection dimensional correction.
         tuple of ndarray
             When ``return_terms`` is True, returns the diagnostic terms
             ``(convec_term, drying, heat_cond, heat_loss_emp)`` instead.
             ``convec_term`` is the raw ``u_gas * dTg_dz`` diagnostic
-            [kg*K/m**3/s]; the remaining entries are gas-temperature-rate
-            contributions [K/s].
+            [kg*K/m**3/s]. ``drying`` is the condensed-temperature latent
+            contribution [K/s]; ``heat_cond`` and ``heat_loss_emp`` are
+            gas-temperature-rate contributions [K/s].
 
         Notes
         -----
@@ -486,10 +509,11 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
         gas-convection discretization is preserved in this branch; ``dTg_dz``
         [kg*K/m**4] and ``conv_term`` [J*kg/m**6/s] are annotated as
         implemented so that the remaining dimensional debt is explicit without
-        expanding the #24 behavior change.
+        expanding the #24 behavior change. The legacy gas molecular-weight
+        surrogate is documented as issue #28 rather than changed here.
         """
 
-        mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)  # [g/mol]
+        mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw)  # legacy MW surrogate [g/mol]
         # ----- Reading inputs
         temp_gas_inputs = inputs['temp']  # [K]
 
@@ -525,7 +549,9 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
 
         conv_term = -u_gas * dTg_dz * cpg_mix * rho_gas  # [J*kg/m**6/s]
 
-        dTg_dt = (conv_term + sensible_heat - heat_transf - heat_loss) / denom_gas  # [K/s]
+        dTg_dt = (
+            conv_term + sensible_heat - heat_transf - heat_loss
+        ) / denom_gas  # legacy solver output [K/s]; issue #37 owns units.
 
         # ----- Condensed phases equations
         dens_liq = self.rho_liq  # [kg/m**3]
@@ -540,7 +566,7 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
 
         if return_terms:
             self.convec_term = u_gas * dTg_dz  # [kg*K/m**3/s]
-            self.drying = drying_terms/ denom_gas  # [K/s]
+            self.drying = drying_terms / denom_cond  # [K/s]
             self.heat_cond = heat_transf/ denom_gas  # [K/s]
             self.heat_loss_emp = heat_loss/ denom_gas  # [K/s]
 
@@ -552,8 +578,38 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq,
 
     def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True,
                    verbose=True, sundials_opts=None):
+        """Initialize and integrate the drying model.
+
+        Parameters
+        ----------
+        deltaP : float
+            Pressure drop across the drying cake and medium [Pa].
+        runtime : float, optional
+            Duration to simulate from the current elapsed time [s].
+        time_grid : ndarray, optional
+            Output time grid for the CVode simulation [s].
+        any_event : bool, optional
+            Whether any configured state event can terminate the solve [-].
+        verbose : bool, optional
+            If False, suppress solver output [-].
+        sundials_opts : dict, optional
+            Solver option names and values forwarded to CVode.
+
+        Returns
+        -------
+        time : ndarray
+            Simulated time points [s].
+        states : ndarray
+            Flattened state history. Saturation and mass-fraction columns are
+            dimensionless [-]; temperature columns are [K].
+
+        Notes
+        -----
+        Cake permeability is computed as
+        ``1 / (alpha * rho_sol * (1 - porosity))`` [m**2].
+        """
         
-        p_atm=101325
+        p_atm=101325  # [Pa]
         # ---------- Initialization
         # Volatile components
         idx_liquid = np.arange(0, self.Liquid_1.num_species)
@@ -574,15 +630,15 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True,
         # y_gas_init = np.tile(self.Vapor_1.mole_frac, (self.num_nodes,1))
         # x_liq_init = self.CakePhase.Liquid_1.mole_frac[:, idx_volatiles]
 
-        y_gas_init = self.Vapor_1.mass_frac
-        x_liq_init = self.CakePhase.Liquid_1.mass_frac
+        y_gas_init = self.Vapor_1.mass_frac  # [-]
+        x_liq_init = self.CakePhase.Liquid_1.mass_frac  # [-]
 
-        satur_init = self.CakePhase.saturation
+        satur_init = self.CakePhase.saturation  # [-]
 
         # Temperatures
-        temp_cond_init = self.CakePhase.Solid_1.temp
-        temp_gas_init = self.Vapor_1.temp
-        z_cake = self.CakePhase.z_external  # For drying_script_inyoung
+        temp_cond_init = self.CakePhase.Solid_1.temp  # [K]
+        temp_gas_init = self.Vapor_1.temp  # [K]
+        z_cake = self.CakePhase.z_external  # [m], for drying_script_inyoung
         # z_cake = self.CakePhase.z_external # This line for 2MSMPR_Filter.py
 
         if x_liq_init.ndim == 1:
@@ -631,9 +687,9 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True,
         #z_cake = self.cake_height * self.CakePhase.z_external
 
         # Physical properties
-        alpha = self.CakePhase.alpha
-        rho_sol = self.Solid_1.getDensity()
-        porosity = self.CakePhase.porosity
+        alpha = self.CakePhase.alpha  # [m/kg]
+        rho_sol = self.Solid_1.getDensity()  # [kg/m**3]
+        porosity = self.CakePhase.porosity  # [-]
 
         xliq = states_prev[:, num_comp + 1: num_comp + 1 + self.num_volatiles]
         # xliq = states_init[:, num_comp + 1: num_comp + 1 + self.num_volatiles]
@@ -641,36 +697,36 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True,
         xliq_init = np.zeros((self.num_nodes, num_comp))
         xliq_init[:, self.idx_volatiles] = xliq
         rho_liq = self.Liquid_1.getDensity(temp=temp_cond_init,
-                                           mass_frac=xliq_init, basis='mass')
+                                           mass_frac=xliq_init, basis='mass')  # [kg/m**3]
         surf_tens = self.Liquid_1.getSurfTension(temp=temp_cond_init,
-                                                 mass_frac=xliq_init)
+                                                 mass_frac=xliq_init)  # [N/m]
 
-        self.k_perm = 1 / alpha / rho_sol / (1 - porosity)
+        self.k_perm = 1 / alpha / rho_sol / (1 - porosity)  # [m**2]
         # self.rho_liq = rho_liq
         self.rho_sol = rho_sol
         self.porosity = porosity
-        self.cp_sol = self.Solid_1.getCp()
+        self.cp_sol = self.Solid_1.getCp()  # [J/kg/K]
 
         # Mass transfer
         moments = self.Solid_1.getMoments(mom_num=[0, 1, 2, 3, 4])
         sauter_diam = moments[1] / moments[0]  # [m]
 
         # self.a_V = 6 / sauter_diam  # [m**2/m**3]
-        self.a_V = moments[2] * (1 - porosity) / moments[3]
+        self.a_V = moments[2] * (1 - porosity) / moments[3]  # [m**2/m**3]
         # Gas pressure
         # deltaP_media = deltaP*self.resist_medium / \
         #     (alpha*rho_sol*self.cake_height + self.resist_medium)
 
         deltaP_media = deltaP*self.resist_medium / \
             (alpha*rho_sol*(1 - porosity)*self.cake_height +
-              self.resist_medium)
-        deltaP -= deltaP_media
-        self.deltaP = deltaP
-        p_top = p_atm + deltaP
+              self.resist_medium)  # [Pa]
+        deltaP -= deltaP_media  # [Pa]
+        self.deltaP = deltaP  # [Pa]
+        p_top = p_atm + deltaP  # [Pa]
 
-        self.dPg_dz = deltaP / self.cake_height
+        self.dPg_dz = deltaP / self.cake_height  # [Pa/m]
         self.pres_gas = np.linspace(p_top, p_top - deltaP,
-                                    num=self.num_nodes)
+                                    num=self.num_nodes)  # [Pa]
 
         # Irreducible saturation
         x_csd = self.Solid_1.x_distrib #* 1e-6
@@ -680,7 +736,7 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True,
         # rholiq_mass = np.mean(rho_liq[0] * self.Liquid_1.mw_av[0])  # [kg/m**3]
         self.s_inf = get_sat_inf(x_csd, csd, deltaP, porosity,
                                  self.cake_height, mom_zero,
-                                 (np.mean(surf_tens), rho_liq[0]))
+                                 (np.mean(surf_tens), rho_liq[0]))  # [-]
 
         # ---------- Solve model
         model = self.unit_model
@@ -729,6 +785,21 @@ def new_handle(solver, info):
         return time, states
 
     def retrieve_results(self, time, states):
+        """Store drying simulation outputs in a ``DynamicResult``.
+
+        Parameters
+        ----------
+        time : array_like
+            Simulated time points [s].
+        states : ndarray
+            Flattened state history. Saturation and mass-fraction columns are
+            dimensionless [-]; temperature columns are [K].
+
+        Returns
+        -------
+        None
+            The method updates ``result`` and cake-grid metadata in place.
+        """
         time = np.array(time)
         self.timeProf = time
         self.elapsed_time += time[-1]
@@ -749,6 +820,13 @@ def retrieve_results(self, time, states):
         self.CakePhase.z_external = self.z_centers
 
     def flatten_states(self):
+        """Placeholder for future drying-state flattening support.
+
+        Returns
+        -------
+        None
+            This method is not implemented.
+        """
         pass
 
     def plot_profiles(self, times=None, z_pos=None, pick_comp=None, **fig_kw):
@@ -832,6 +910,22 @@ def plot_profiles(self, times=None, z_pos=None, pick_comp=None, **fig_kw):
         return fig, axes
 
     def plot_rates(self, z_pos=0, fig_size=None):
+        """Plot drying rates at one axial position.
+
+        Parameters
+        ----------
+        z_pos : float, optional
+            Axial position used to select the nearest drying node [m].
+        fig_size : tuple, optional
+            Figure size forwarded to matplotlib [-].
+
+        Returns
+        -------
+        fig : matplotlib.figure.Figure
+            Created figure.
+        axis : matplotlib.axes.Axes
+            Axes containing drying-rate profiles [mol/m**3/s].
+        """
         z_idx = np.argmin(abs(z_pos - self.z_centers))
         temp_cond = self.tempLiqProf[:, z_idx]
 
diff --git a/PharmaPy/Plotting.py b/PharmaPy/Plotting.py
index 2e3117c..abb7a26 100644
--- a/PharmaPy/Plotting.py
+++ b/PharmaPy/Plotting.py
@@ -64,6 +64,37 @@ def latexify_name(name, units=False):
     return out
 
 
+def format_unit_label(units):
+    """Return a display label for unit metadata.
+
+    Parameters
+    ----------
+    units : str or None
+        Unit metadata. Bracketed values such as ``[K]`` and ``[-]`` follow the
+        repository documentation convention.
+
+    Returns
+    -------
+    str
+        LaTeX-formatted unit label without metadata brackets. Dimensionless
+        metadata ``[-]`` and empty values return an empty string.
+    """
+    if units is None:
+        return ''
+
+    units = str(units).strip()
+    if units == '':
+        return ''
+
+    if units.startswith('[') and units.endswith(']'):
+        units = units[1:-1].strip()
+
+    if units in ('', '-'):
+        return ''
+
+    return latexify_name(units, units=True)
+
+
 def color_axis(ax, color):
     ax.spines['right'].set_color(color)
     ax.tick_params(axis='y', colors=color, which='both')
@@ -182,8 +213,8 @@ def name_yaxes(ax, states_fstates, names, ylabels, legend):
             ylabel = latexify_name(ylabels[ind])
 
         units = states_fstates[name].get('units', '')
-        if len(units) > 0:
-            unit_name = latexify_name(units, units=True)
+        unit_name = format_unit_label(units)
+        if len(unit_name) > 0:
             ylabel = ylabel + ' (' + unit_name + ')'
 
         axis.set_ylabel(ylabel)
@@ -265,9 +296,8 @@ def plot_function(uo, state_names, axes=None, fig_map=None, ylabels=None,
             ylabel = latexify_name(ylabels[ind])
 
         units = states_and_fstates[name].get('units', '')
-        if len(units) > 0:
-            unit_name = latexify_name(states_and_fstates[name]['units'],
-                                      units=True)
+        unit_name = format_unit_label(units)
+        if len(unit_name) > 0:
             ylabel = ylabel + ' (' + unit_name + ')'
 
         if index_y:
diff --git a/tests/test_drying_latent_heat_factor.py b/tests/test_drying_latent_heat_factor.py
index 6b28ade..50352ff 100644
--- a/tests/test_drying_latent_heat_factor.py
+++ b/tests/test_drying_latent_heat_factor.py
@@ -14,27 +14,27 @@
 
 def test_condensed_energy_uses_single_latent_heat_factor(monkeypatch):
     """Mass drying rates multiplied by latent heat already give full power."""
-    dryer = Drying(number_nodes=3, supercrit_names=["nitrogen"])
+    dryer = Drying(number_nodes=4, supercrit_names=["nitrogen"])
     dryer.idx_volatiles = np.array([0, 2])  # component indices [-]
     dryer.idx_supercrit = np.array([1])  # component indices [-]
     dryer.porosity = 0.5  # [-]
     dryer.rho_sol = 1000.0  # [kg/m**3]
     dryer.cp_sol = 1000.0  # [J/kg/K]
-    dryer.rho_liq = np.full(3, 800.0)  # [kg/m**3]
+    dryer.rho_liq = np.full(4, 800.0)  # [kg/m**3]
     dryer.h_T_j = 0.0  # [W/m**2/K]
     dryer.a_V = 1.0  # [m**2/m**3]
     dryer.h_T_loss = 0.0  # [W/m**2/K]
     dryer.cake_height = 1.0  # [m]
     dryer.T_ambient = 298.15  # [K]
-    dryer.dz = np.ones(3)  # [m]
+    dryer.dz = np.ones(4)  # [m]
 
     dryer.Vapor_1 = SimpleNamespace(
         mw=np.array([28.0, 28.0, 28.0]),  # [g/mol]
-        getCp=lambda temp, mass_frac, basis: np.full(3, 1000.0),  # [J/kg/K]
+        getCp=lambda temp, mass_frac, basis: np.full(4, 1000.0),  # [J/kg/K]
         getHeatVaporization=lambda temp, basis: np.array([2.0e6, 1.0e6]),  # [J/kg]
     )
     dryer.Liquid_1 = SimpleNamespace(
-        getCp=lambda temp, mass_frac, basis: np.full(3, 2000.0),  # [J/kg/K]
+        getCp=lambda temp, mass_frac, basis: np.full(4, 2000.0),  # [J/kg/K]
     )
 
     monkeypatch.setattr(
@@ -43,23 +43,26 @@ def test_condensed_energy_uses_single_latent_heat_factor(monkeypatch):
         lambda values, boundary_cond: np.zeros(values.size + 1),  # [K]
     )
 
-    temp_gas = np.full(3, 295.0)  # [K], matches energy_balance wet-bulb default
-    temp_sol = np.array([299.0, 304.0, 309.0])  # [K]
-    satur = np.full(3, 0.5)  # [-]
+    temp_gas = np.full(4, 295.0)  # [K], matches energy_balance wet-bulb default
+    temp_sol = np.array([299.0, 304.0, 309.0, 314.0])  # [K]
+    satur = np.full(4, 0.5)  # [-]
     y_gas = np.array([
         [0.10, 0.80, 0.10],
         [0.20, 0.70, 0.10],
         [0.30, 0.60, 0.10],
+        [0.40, 0.50, 0.10],
     ])  # gas mass fractions [-]
     x_liq = np.array([
         [0.0, 1.0],
         [0.0, 1.0],
         [0.0, 1.0],
+        [0.0, 1.0],
     ])  # liquid volatile mass fractions [-]
     dry_rate = np.array([
         [0.036, 0.0, 0.184],
         [0.054, 0.0, 0.230],
         [0.018, 0.0, 0.092],
+        [0.027, 0.0, 0.046],
     ])  # [kg/m**3/s]
 
     _, dTcond_dt = dryer.energy_balance(
@@ -69,8 +72,8 @@ def test_condensed_energy_uses_single_latent_heat_factor(monkeypatch):
         satur=satur,
         y_gas=y_gas,
         x_liq=x_liq,
-        u_gas=np.zeros(3),  # [m/s]
-        rho_gas=np.ones(3),  # [kg/m**3]
+        u_gas=np.zeros(4),  # [m/s]
+        rho_gas=np.ones(4),  # [kg/m**3]
         dry_rate=dry_rate,
         inputs={"temp": 295.0},  # [K]
     )
@@ -79,5 +82,21 @@ def test_condensed_energy_uses_single_latent_heat_factor(monkeypatch):
         -256000.0 / 900000.0,
         -338000.0 / 900000.0,
         -128000.0 / 900000.0,
+        -100000.0 / 900000.0,
     ])  # [K/s]
     np.testing.assert_allclose(dTcond_dt, expected_condensed_temperature_rate)
+
+    _, drying_term, _, _ = dryer.energy_balance(
+        time=0.0,
+        temp_gas=temp_gas,
+        temp_sol=temp_sol,
+        satur=satur,
+        y_gas=y_gas,
+        x_liq=x_liq,
+        u_gas=np.zeros(4),  # [m/s]
+        rho_gas=np.ones(4),  # [kg/m**3]
+        dry_rate=dry_rate,
+        inputs={"temp": 295.0},  # [K]
+        return_terms=True,
+    )
+    np.testing.assert_allclose(drying_term, -expected_condensed_temperature_rate)
diff --git a/tests/test_plotting_unit_labels.py b/tests/test_plotting_unit_labels.py
new file mode 100644
index 0000000..03f0b7b
--- /dev/null
+++ b/tests/test_plotting_unit_labels.py
@@ -0,0 +1,35 @@
+"""Regression tests for plotting unit metadata labels."""
+
+import matplotlib.pyplot as plt
+
+from PharmaPy.Plotting import format_unit_label, latexify_name, name_yaxes
+
+
+def test_format_unit_label_accepts_bracketed_metadata():
+    """Bracketed source metadata renders without nested brackets."""
+    temp_units = "[K]"
+    density_units = "[kg/m**3]"
+
+    assert format_unit_label(temp_units) == latexify_name(temp_units[1:-1], units=True)
+    assert format_unit_label(density_units) == latexify_name(
+        density_units[1:-1], units=True)
+
+
+def test_name_yaxes_suppresses_dimensionless_metadata():
+    """Dimensionless ``[-]`` metadata does not add a plot-label suffix."""
+    temp_units = "[K]"
+    fig, axes = plt.subplots(1, 2)
+    states = {
+        "temp": {"units": temp_units},
+        "x_liq": {"units": "[-]"},
+    }
+
+    try:
+        name_yaxes(axes, states, ("temp", "x_liq"), ("T", "x_liq"), True)
+
+        assert axes[0].get_ylabel() == (
+            f"{latexify_name('T')} ({latexify_name(temp_units[1:-1], units=True)})"
+        )
+        assert axes[1].get_ylabel() == latexify_name("x_liq")
+    finally:
+        plt.close(fig)