Check x bounds before calculating i0 (similar to what is done for Eeq)

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## [Unreleased](https://github.com/NREL/BATMODS-lite/)
 
 ### New Features
+- Check bounds of intercalation fraction before calculating exchange current density ([#19](https://github.com/NREL/BATMODS-lite/pull/19))
 - Update `materials` module with new LFP properties, from ICI ([#18](https://github.com/NREL/BATMODS-lite/pull/18))
 - Use `CubicSpline` interpolations for LFP properties ([#15](https://github.com/NREL/BATMODS-lite/pull/15))
 - Add degradation parameters for the electrolyte (applies to P2D only) ([#14](https://github.com/NREL/BATMODS-lite/pull/14))

src/bmlite/materials/_graphite.py

@@ -1,3 +1,5 @@
+from numbers import Real
+
 import numpy as np
 
 
@@ -93,6 +95,15 @@ def get_i0(self, x: float | np.ndarray, C_Li: float | np.ndarray,
 
         c = Constants()
 
+        # Avoid floating point errors
+        if isinstance(x, Real):
+            if (x < 0 and self.alpha_c < 1) or (x > 1 and self.alpha_a < 1):
+                raise ValueError('x is out of [0, 1] during i0 calculation')
+        elif isinstance(x, np.ndarray):
+            if ((any(x.flatten() < 0) and self.alpha_c < 1)
+               or (any(x.flatten() > 1) and self.alpha_a < 1)):
+                raise ValueError('x is out of [0, 1] during i0 calculation')
+
         i0 = 2.5 * 0.27 * np.exp(-30e6 / c.R * (1 / T - 1 / 303.15)) \
             * C_Li**self.alpha_a * (self.Li_max * x)**self.alpha_c \
             * (self.Li_max - self.Li_max * x)**self.alpha_a

src/bmlite/materials/_lfp.py

@@ -1,3 +1,5 @@
+from numbers import Real
+
 import numpy as np
 
 
@@ -101,6 +103,15 @@ def get_i0(self, x: float | np.ndarray, C_Li: float | np.ndarray,
 
         c = Constants()
 
+        # Avoid floating point errors
+        if isinstance(x, Real):
+            if (x < 0 and self.alpha_c < 1) or (x > 1 and self.alpha_a < 1):
+                raise ValueError('x is out of [0, 1] during i0 calculation')
+        elif isinstance(x, np.ndarray):
+            if ((any(x.flatten() < 0) and self.alpha_c < 1)
+               or (any(x.flatten() > 1) and self.alpha_a < 1)):
+                raise ValueError('x is out of [0, 1] during i0 calculation')
+
         i0 = 0.27 * np.exp(-30e6 / c.R * (1 / T - 1 / 303.15)) \
            * C_Li**self.alpha_a * (self.Li_max * x)**self.alpha_c \
            * (self.Li_max - self.Li_max * x)**self.alpha_a