Fix wbar and add tests for ebar and wbar

pttools/bubble/thermo.py

@@ -123,11 +123,13 @@ def ubarf2(v: np.ndarray, w: np.ndarray, xi: np.ndarray, v_wall: float, ek_bva:
     return ek_bva / w[-1]
 
 
-def wbar(w: np.ndarray, xi: np.ndarray, v_wall: float, wn: float):
-    logger.warning("wbar is not properly tested and may return false results.")
+def wbar(w: np.ndarray, xi: np.ndarray, v_wall: float, wn: float) -> float:
+    # https://stackoverflow.com/a/8768734
     w_reverse = w[::-1]
-    i_max = len(w_reverse) - np.argmax(w_reverse != w[-1]) - 1
-    ret = 1/v_wall**3 * np.trapezoid(w[:i_max], xi[:i_max]**3)
+    i_max = w.size - np.argmax(w_reverse != w[-1]) - 1
+    if i_max == 0:
+        i_max = -1
+    ret = 1/(xi[i_max]**3) * np.trapezoid(w[:i_max+1], xi[:i_max+1]**3)
     if not (ret is None or np.isnan(ret)) and ret <= wn:
         logger.warning(f"Should have wbar > wn. Got: wbar={wn}, wn={wn}")
     return ret

tests/bubble/test_thermo.py

@@ -4,6 +4,7 @@
 import numpy as np
 
 from pttools.bubble.bubble import Bubble
+from pttools.bubble import thermo
 from pttools.models.model import Model
 from pttools.models.bag import BagModel
 from pttools.models.const_cs import ConstCSModel
@@ -30,6 +31,19 @@ def setUpClass(cls) -> None:
             for v_wall, alpha_n in zip(cls.V_WALLS, cls.ALPHA_NS)
         ]
 
+    def test_ebar(self):
+        assert_allclose(
+            [thermo.ebar(model=bubble.model, wn=bubble.wn) for bubble in self.bubbles],
+            [bubble.en() for bubble in self.bubbles]
+        )
+
+    def test_wbar(self):
+        """If there is no bubble, then wbar=wn"""
+        assert_allclose(
+            [thermo.wbar(w=np.ones_like(bubble.w)*bubble.wn, xi=bubble.xi, v_wall=bubble.v_wall, wn=bubble.wn) for bubble in self.bubbles],
+            [bubble.wn for bubble in self.bubbles]
+        )
+
     def test_kappa(self):
         assert_allclose([bubble.kappa for bubble in self.bubbles], self.KAPPA_REF, rtol=1.5e-2)