diff --git a/apps/HeatEquationKernelWithMaterial.py b/apps/HeatEquationKernelWithMaterial.py
index 9daeff693c63de02f89bdb35b7e5b7e3daf3396e..af283966d048492bdb03b9a6ad786524d1ff7e3e 100644
--- a/apps/HeatEquationKernelWithMaterial.py
+++ b/apps/HeatEquationKernelWithMaterial.py
@@ -4,8 +4,9 @@ from importlib.resources import files
 from pystencilssfg import SourceFileGenerator
 from sfg_walberla import Sweep
 from pymatlib.core.yaml_parser import create_alloy_from_yaml
-from pymatlib.core.assignment_converter import assignment_converter
+from pymatlib.core.property_array_extractor import PropertyArrayExtractor
 from pymatlib.core.codegen.interpolation_array_container import InterpolationArrayContainer
+from pymatlib.core.assignment_converter import assignment_converter
 
 with SourceFileGenerator() as sfg:
     data_type = "float64"  # if ctx.double_accuracy else "float32"
@@ -22,8 +23,9 @@ with SourceFileGenerator() as sfg:
     heat_pde_discretized = heat_pde_discretized.args[1] + heat_pde_discretized.args[0].simplify()
 
     yaml_path = files('pymatlib.data.alloys.SS304L').joinpath('SS304L.yaml')
-    mat = create_alloy_from_yaml(yaml_path, u.center())
-    arr_container = InterpolationArrayContainer.from_material("SS304L", mat)
+    mat, temperature_array = create_alloy_from_yaml(yaml_path, u.center())
+    array_extractor = PropertyArrayExtractor(mat, temperature_array, u.center)
+    arr_container = InterpolationArrayContainer("SS304L", temperature_array, array_extractor.energy_density_array)
     sfg.generate(arr_container)
 
     # Convert assignments to pystencils format
diff --git a/src/pymatlib/core/alloy.py b/src/pymatlib/core/alloy.py
index 190252a6dd182c2ca94a9f2f86088e7ef3b37200..f7fc563b8470cd61092dae4792634b099a0fe4a1 100644
--- a/src/pymatlib/core/alloy.py
+++ b/src/pymatlib/core/alloy.py
@@ -60,15 +60,11 @@ class Alloy:
     temperature_boil: float = field(init=False)
 
     # Optional properties with default values
-    base_temperature: float = None
-    base_density: float = None
     density: PropertyTypes = None
     dynamic_viscosity: PropertyTypes = None
     energy_density: PropertyTypes = None
     energy_density_solidus: float = None
     energy_density_liquidus: float = None
-    energy_density_temperature_array: np.ndarray = field(default_factory=lambda: np.array([]))
-    energy_density_array: np.ndarray = field(default_factory=lambda: np.array([]))
     heat_capacity: PropertyTypes = None
     heat_conductivity: PropertyTypes = None
     kinematic_viscosity: PropertyTypes = None
@@ -76,8 +72,6 @@ class Alloy:
     latent_heat_of_vaporization: PropertyTypes = None
     specific_enthalpy: PropertyTypes = None
     surface_tension: PropertyTypes = None
-    temperature: PropertyTypes = None
-    temperature_array: np.ndarray = field(default_factory=lambda: np.array([]))
     thermal_diffusivity: PropertyTypes = None
     thermal_expansion_coefficient: PropertyTypes = None
 
diff --git a/src/pymatlib/core/codegen/interpolation_array_container.py b/src/pymatlib/core/codegen/interpolation_array_container.py
index 240ea1c9e93aad6811e5dca351829abb0285f696..a3d738317f243b7aa5e42f23bf65daa1ca9276a5 100644
--- a/src/pymatlib/core/codegen/interpolation_array_container.py
+++ b/src/pymatlib/core/codegen/interpolation_array_container.py
@@ -7,27 +7,27 @@ from pymatlib.core.interpolators import prepare_interpolation_arrays
 
 
 class InterpolationArrayContainer(CustomGenerator):
-    """Container for energy-temperature interpolation arrays and methods.
+    """Container for x-y interpolation arrays and methods.
 
-    This class stores temperature and energy density arrays and generates C++ code
-    for efficient bilateral conversion between these properties. It supports both
+    This class stores x and y arrays and generates C++ code
+    for efficient conversion to compute y for a given x. It supports both
     binary search interpolation (O(log n)) and double lookup interpolation (O(1))
     with automatic method selection based on data characteristics.
 
     Attributes:
         name (str): Name for the generated C++ class.
-        T_array (np.ndarray): Array of temperature values (must be monotonically increasing).
-        E_array (np.ndarray): Array of energy density values corresponding to T_array.
+        x_array (np.ndarray): Array of x values (must be monotonically increasing).
+        y_array (np.ndarray): Array of y values corresponding to x_array.
         method (str): Interpolation method selected ("binary_search" or "double_lookup").
-        T_bs (np.ndarray): Temperature array prepared for binary search.
-        E_bs (np.ndarray): Energy array prepared for binary search.
+        x_bs (np.ndarray): x array prepared for binary search.
+        y_bs (np.ndarray): y array prepared for binary search.
         has_double_lookup (bool): Whether double lookup interpolation is available.
 
     If has_double_lookup is True, the following attributes are also available:
-        T_eq (np.ndarray): Equidistant temperature array for double lookup.
-        E_neq (np.ndarray): Non-equidistant energy array for double lookup.
-        E_eq (np.ndarray): Equidistant energy array for double lookup.
-        inv_delta_E_eq (float): Inverse of the energy step size for double lookup.
+        x_eq (np.ndarray): Equidistant x array for double lookup.
+        y_neq (np.ndarray): Non-equidistant y array for double lookup.
+        y_eq (np.ndarray): Equidistant y array for double lookup.
+        inv_delta_y_eq (float): Inverse of the y step size for double lookup.
         idx_map (np.ndarray): Index mapping array for double lookup.
 
     Examples:
@@ -44,14 +44,14 @@ class InterpolationArrayContainer(CustomGenerator):
         >>>     container = InterpolationArrayContainer("MyMaterial", T, E)
         >>>     sfg.generate(container)
     """
-    def __init__(self, name: str, temperature_array: np.ndarray, energy_density_array: np.ndarray):
+    def __init__(self, name: str, x_array: np.ndarray, y_array: np.ndarray):
         """Initialize the interpolation container.
         Args:
             name (str): Name for the generated C++ class.
-            temperature_array (np.ndarray): Array of temperature values (K).
+            x_array (np.ndarray): Array of x values.
                 Must be monotonically increasing.
-            energy_density_array (np.ndarray): Array of energy density values (J/m³)
-                corresponding to temperature_array.
+            y_array (np.ndarray): Array of y values
+                corresponding to x_array.
         Raises:
             ValueError: If arrays are empty, have different lengths, or are not monotonic.
             TypeError: If name is not a string or arrays are not numpy arrays.
@@ -65,28 +65,28 @@ class InterpolationArrayContainer(CustomGenerator):
         if not re.match(r'^[a-zA-Z_][a-zA-Z0-9_]*$', name):
             raise ValueError(f"'{name}' is not a valid C++ class name")
 
-        if not isinstance(temperature_array, np.ndarray) or not isinstance(energy_density_array, np.ndarray):
+        if not isinstance(x_array, np.ndarray) or not isinstance(y_array, np.ndarray):
             raise TypeError("Temperature and energy arrays must be numpy arrays")
 
         self.name = name
-        self.T_array = temperature_array
-        self.E_array = energy_density_array
+        self.x_array = x_array
+        self.y_array = y_array
 
         # Prepare arrays and determine best method
         try:
-            self.data = prepare_interpolation_arrays(T_array=self.T_array, E_array=self.E_array, verbose=False)
+            self.data = prepare_interpolation_arrays(x_array=self.x_array, y_array=self.y_array, verbose=False)
             self.method = self.data["method"]
 
             # Store arrays for binary search (always available)
-            self.T_bs = self.data["T_bs"]
-            self.E_bs = self.data["E_bs"]
+            self.x_bs = self.data["x_bs"]
+            self.y_bs = self.data["y_bs"]
 
             # Store arrays for double lookup if available
             if self.method == "double_lookup":
-                self.T_eq = self.data["T_eq"]
-                self.E_neq = self.data["E_neq"]
-                self.E_eq = self.data["E_eq"]
-                self.inv_delta_E_eq = self.data["inv_delta_E_eq"]
+                self.x_eq = self.data["x_eq"]
+                self.y_neq = self.data["y_neq"]
+                self.y_eq = self.data["y_eq"]
+                self.inv_delta_y_eq = self.data["inv_delta_y_eq"]
                 self.idx_map = self.data["idx_map"]
                 self.has_double_lookup = True
             else:
@@ -94,18 +94,6 @@ class InterpolationArrayContainer(CustomGenerator):
         except Exception as e:
             raise ValueError(f"Failed to prepare interpolation arrays: {e}") from e
 
-    @classmethod
-    def from_material(cls, name: str, material):
-        """Create an interpolation container from a material object.
-        Args:
-            name (str): Name for the generated C++ class.
-            material: Material object with temperature and energy properties.
-                Must have energy_density_temperature_array and energy_density_array attributes.
-        Returns:
-            InterpolationArrayContainer: Container with arrays for interpolation.
-        """
-        return cls(name, material.energy_density_temperature_array, material.energy_density_array)
-
     def _generate_binary_search(self, sfg: SfgComposer):
         """Generate code for binary search interpolation.
         Args:
@@ -113,19 +101,19 @@ class InterpolationArrayContainer(CustomGenerator):
         Returns:
             list: List of public members for the C++ class.
         """
-        T_bs_arr_values = ", ".join(str(v) for v in self.T_bs)
-        E_bs_arr_values = ", ".join(str(v) for v in self.E_bs)
+        x_bs_arr_values = ", ".join(str(v) for v in self.x_bs)
+        y_bs_arr_values = ", ".join(str(v) for v in self.y_bs)
 
-        E_target = sfg.var("E_target", "double")
+        y_target = sfg.var("y_target", "double")
 
         return [
             # Binary search arrays
-            f"static constexpr std::array< double, {self.T_bs.shape[0]} > T_bs {{ {T_bs_arr_values} }}; \n"
-            f"static constexpr std::array< double, {self.E_bs.shape[0]} > E_bs {{ {E_bs_arr_values} }}; \n",
+            f"static constexpr std::array< double, {self.x_bs.shape[0]} > x_bs {{ {x_bs_arr_values} }}; \n"
+            f"static constexpr std::array< double, {self.y_bs.shape[0]} > y_bs {{ {y_bs_arr_values} }}; \n",
 
             # Binary search method
             sfg.method("interpolateBS", returns=PsCustomType("[[nodiscard]] double"), inline=True, const=True)(
-                sfg.expr("return interpolate_binary_search_cpp({}, *this);", E_target)
+                sfg.expr("return interpolate_binary_search_cpp({}, *this);", y_target)
             )
         ]
 
@@ -139,24 +127,24 @@ class InterpolationArrayContainer(CustomGenerator):
         if not self.has_double_lookup:
             return []
 
-        T_eq_arr_values = ", ".join(str(v) for v in self.T_eq)
-        E_neq_arr_values = ", ".join(str(v) for v in self.E_neq)
-        E_eq_arr_values = ", ".join(str(v) for v in self.E_eq)
+        x_eq_arr_values = ", ".join(str(v) for v in self.x_eq)
+        y_neq_arr_values = ", ".join(str(v) for v in self.y_neq)
+        y_eq_arr_values = ", ".join(str(v) for v in self.y_eq)
         idx_mapping_arr_values = ", ".join(str(v) for v in self.idx_map)
 
-        E_target = sfg.var("E_target", "double")
+        y_target = sfg.var("y_target", "double")
 
         return [
             # Double lookup arrays
-            f"static constexpr std::array< double, {self.T_eq.shape[0]} > T_eq {{ {T_eq_arr_values} }}; \n"
-            f"static constexpr std::array< double, {self.E_neq.shape[0]} > E_neq {{ {E_neq_arr_values} }}; \n"
-            f"static constexpr std::array< double, {self.E_eq.shape[0]} > E_eq {{ {E_eq_arr_values} }}; \n"
-            f"static constexpr double inv_delta_E_eq = {self.inv_delta_E_eq}; \n"
+            f"static constexpr std::array< double, {self.x_eq.shape[0]} > x_eq {{ {x_eq_arr_values} }}; \n"
+            f"static constexpr std::array< double, {self.y_neq.shape[0]} > y_neq {{ {y_neq_arr_values} }}; \n"
+            f"static constexpr std::array< double, {self.y_eq.shape[0]} > y_eq {{ {y_eq_arr_values} }}; \n"
+            f"static constexpr double inv_delta_y_eq = {self.inv_delta_y_eq}; \n"
             f"static constexpr std::array< int, {self.idx_map.shape[0]} > idx_map {{ {idx_mapping_arr_values} }}; \n",
 
             # Double lookup method
             sfg.method("interpolateDL", returns=PsCustomType("[[nodiscard]] double"), inline=True, const=True)(
-                sfg.expr("return interpolate_double_lookup_cpp({}, *this);", E_target)
+                sfg.expr("return interpolate_double_lookup_cpp({}, *this);", y_target)
             )
         ]
 
@@ -178,17 +166,17 @@ class InterpolationArrayContainer(CustomGenerator):
             public_members.extend(self._generate_double_lookup(sfg))
 
         # Add interpolate method that uses recommended approach
-        E_target = sfg.var("E_target", "double")
+        y_target = sfg.var("y_target", "double")
         if self.has_double_lookup:
             public_members.append(
                 sfg.method("interpolate", returns=PsCustomType("[[nodiscard]] double"), inline=True, const=True)(
-                    sfg.expr("return interpolate_double_lookup_cpp({}, *this);", E_target)
+                    sfg.expr("return interpolate_double_lookup_cpp({}, *this);", y_target)
                 )
             )
         else:
             public_members.append(
                 sfg.method("interpolate", returns=PsCustomType("[[nodiscard]] double"), inline=True, const=True)(
-                    sfg.expr("return interpolate_binary_search_cpp({}, *this);", E_target)
+                    sfg.expr("return interpolate_binary_search_cpp({}, *this);", y_target)
                 )
             )
 
diff --git a/src/pymatlib/core/cpp/include/pymatlib_interpolators/interpolate_binary_search_cpp.h b/src/pymatlib/core/cpp/include/pymatlib_interpolators/interpolate_binary_search_cpp.h
index 62dc9db7597a5308d91882fd9174ed773a7215e7..49fb0d6474cb11b5ff72b27fd2fd078d9860ff5a 100644
--- a/src/pymatlib/core/cpp/include/pymatlib_interpolators/interpolate_binary_search_cpp.h
+++ b/src/pymatlib/core/cpp/include/pymatlib_interpolators/interpolate_binary_search_cpp.h
@@ -5,15 +5,15 @@
 
 template<typename ArrayContainer>
 double interpolate_binary_search_cpp(
-    double E_target,
+    double y_target,
     const ArrayContainer& arrs) {
 
     static constexpr double EPSILON = 1e-6;
-    const size_t n = arrs.T_bs.size();
+    const size_t n = arrs.x_bs.size();
 
     // Quick boundary checks
-    if (E_target <= arrs.E_bs[0]) return arrs.T_bs[0];
-    if (E_target >= arrs.E_bs.back()) return arrs.T_bs.back();
+    if (y_target <= arrs.y_bs[0]) return arrs.x_bs[0];
+    if (y_target >= arrs.y_bs.back()) return arrs.x_bs.back();
 
     // Binary search
     size_t left = 0;
@@ -22,11 +22,11 @@ double interpolate_binary_search_cpp(
     while (right - left > 1) {
         const size_t mid = left + (right - left) / 2;
 
-        if (std::abs(arrs.E_bs[mid] - E_target) < EPSILON) {
-            return arrs.T_bs[mid];
+        if (std::abs(arrs.y_bs[mid] - y_target) < EPSILON) {
+            return arrs.x_bs[mid];
         }
 
-        if (arrs.E_bs[mid] > E_target) {
+        if (arrs.y_bs[mid] > y_target) {
             right = mid;
         } else {
             left = mid;
@@ -34,11 +34,11 @@ double interpolate_binary_search_cpp(
     }
 
     // Linear interpolation
-    const double x0 = arrs.E_bs[right];
-    const double x1 = arrs.E_bs[left];
-    const double y0 = arrs.T_bs[right];
-    const double y1 = arrs.T_bs[left];
+    const double x1 = arrs.y_bs[right];
+    const double x2 = arrs.y_bs[left];
+    const double y1 = arrs.x_bs[right];
+    const double y2 = arrs.x_bs[left];
 
-    const double slope = (y1 - y0) / (x1 - x0);
-    return y0 + slope * (E_target - x0);
+    const double slope = (y2 - y1) / (x2 - x1);
+    return y1 + slope * (y_target - x1);
 }
diff --git a/src/pymatlib/core/cpp/include/pymatlib_interpolators/interpolate_double_lookup_cpp.h b/src/pymatlib/core/cpp/include/pymatlib_interpolators/interpolate_double_lookup_cpp.h
index d60134ed6b8e53129857ac9036587efe8342e258..15010fc6969ab482ca00aa922811c58a5a492463 100644
--- a/src/pymatlib/core/cpp/include/pymatlib_interpolators/interpolate_double_lookup_cpp.h
+++ b/src/pymatlib/core/cpp/include/pymatlib_interpolators/interpolate_double_lookup_cpp.h
@@ -5,33 +5,33 @@
 
 template<typename ArrayContainer>
 double interpolate_double_lookup_cpp(
-    double E_target,
+    double y_target,
     const ArrayContainer& arrs) {
 
     // Handle boundary cases
-    if (E_target <= arrs.E_neq[0]) return arrs.T_eq[0];
-    if (E_target >= arrs.E_neq.back()) return arrs.T_eq.back();
+    if (y_target <= arrs.y_neq[0]) return arrs.x_eq[0];
+    if (y_target >= arrs.y_neq.back()) return arrs.x_eq.back();
 
     // Calculate index in equidistant grid
-    int idx_E_eq = static_cast<int>((E_target - arrs.E_eq[0]) * arrs.inv_delta_E_eq);
-    idx_E_eq = std::min(idx_E_eq, static_cast<int>(arrs.idx_map.size() - 1));
+    int idx_y_eq = static_cast<int>((y_target - arrs.y_eq[0]) * arrs.inv_delta_y_eq);
+    idx_y_eq = std::min(idx_y_eq, static_cast<int>(arrs.idx_map.size() - 1));
 
     // Get index from mapping
-    int idx_E_neq = arrs.idx_map[idx_E_eq];
+    int idx_y_neq = arrs.idx_map[idx_y_eq];
 
     // Make sure we don't go out of bounds
-    idx_E_neq = std::min(idx_E_neq, static_cast<int>(arrs.E_neq.size() - 2));
+    idx_y_neq = std::min(idx_y_neq, static_cast<int>(arrs.y_neq.size() - 2));
 
     // Adjust index if needed
-    idx_E_neq += arrs.E_neq[idx_E_neq + 1] < E_target;
+    idx_y_neq += arrs.y_neq[idx_y_neq + 1] < y_target;
 
     // Get interpolation points
-    const double E1 = arrs.E_neq[idx_E_neq];
-    const double E2 = arrs.E_neq[idx_E_neq + 1];
-    const double T1 = arrs.T_eq[idx_E_neq];
-    const double T2 = arrs.T_eq[idx_E_neq + 1];
+    const double y1 = arrs.y_neq[idx_y_neq];
+    const double y2 = arrs.y_neq[idx_y_neq + 1];
+    const double x1 = arrs.x_eq[idx_y_neq];
+    const double x2 = arrs.x_eq[idx_y_neq + 1];
 
     // Linear interpolation
-    const double slope = (T2 - T1) / (E2 - E1);
-    return T1 + slope * (E_target - E1);
+    const double slope = (x2 - x1) / (y2 - y1);
+    return x1 + slope * (y_target - y1);
 }
diff --git a/src/pymatlib/core/data_handler.py b/src/pymatlib/core/data_handler.py
index dd2b189078de8dbb4f5ea2a4e589942e449a8302..4e4c791eb75454bb02607f825d9fe54d0bac01e1 100644
--- a/src/pymatlib/core/data_handler.py
+++ b/src/pymatlib/core/data_handler.py
@@ -141,7 +141,7 @@ def read_data_from_file(file_config: Union[str, Dict], header: bool = True) -> T
         temp_col = file_config['temp_col']
         prop_col = file_config['prop_col']
 
-    print(f"Reading data from file: {file_path}")
+    # print(f"Reading data from file: {file_path}")
 
     if file_path.endswith('.xlsx'):
         df = pd.read_excel(file_path, header=0 if header else None)
diff --git a/src/pymatlib/core/interpolators.py b/src/pymatlib/core/interpolators.py
index 0286ff899f89d8f5329104c7cc61cab04cbd5622..64e07a4b3594c36621302b60c8e56964f03695d3 100644
--- a/src/pymatlib/core/interpolators.py
+++ b/src/pymatlib/core/interpolators.py
@@ -274,93 +274,92 @@ def interpolate_binary_search(
     return float(y0 + (y1 - y0) * (h_in - x0) / (x1 - x0))
 
 
-def E_eq_from_E_neq(E_neq: np.ndarray) -> Tuple[np.ndarray, float]:
-    delta_min: float = np.min(np.diff(E_neq))
-    # delta_min: float = np.min(np.abs(np.diff(E_neq)))
+def y_eq_from_y_neq(y_neq: np.ndarray) -> Tuple[np.ndarray, float]:
+    delta_min: float = np.min(np.diff(y_neq))
     if delta_min < 1.:
-        raise ValueError(f"Energy density array points are very closely spaced, delta = {delta_min}")
-    delta_E_eq = max(np.floor(delta_min * 0.95), 1.)
-    E_eq = np.arange(E_neq[0], E_neq[-1] + delta_E_eq, delta_E_eq, dtype=np.float64)
-    inv_delta_E_eq: float = float(1. / (E_eq[1] - E_eq[0]))
-    return E_eq, inv_delta_E_eq
-
-
-def create_idx_mapping(E_neq: np.ndarray, E_eq: np.ndarray) -> np.ndarray:
-    """idx_map = np.zeros(len(E_eq), dtype=int)
-    for i, e in enumerate(E_eq):
-        idx = int(np.searchsorted(E_neq, e) - 1)
-        idx = max(0, min(idx, len(E_neq) - 2))  # Bound check
+        raise ValueError(f"Array points are very closely spaced, delta = {delta_min}")
+    delta_y_eq = max(np.floor(delta_min * 0.95), 1.)
+    y_eq = np.arange(y_neq[0], y_neq[-1] + delta_y_eq, delta_y_eq, dtype=np.float64)
+    inv_delta_y_eq: float = float(1. / (y_eq[1] - y_eq[0]))
+    return y_eq, inv_delta_y_eq
+
+
+def create_idx_mapping(y_neq: np.ndarray, y_eq: np.ndarray) -> np.ndarray:
+    """idx_map = np.zeros(len(y_eq), dtype=int)
+    for i, e in enumerate(y_eq):
+        idx = int(np.searchsorted(y_neq, e) - 1)
+        idx = max(0, min(idx, len(y_neq) - 2))  # Bound check
         idx_map[i] = idx"""
-    # idx_map = np.searchsorted(E_neq, E_eq) - 1
-    idx_map = np.searchsorted(E_neq, E_eq, side='right') - 1
-    idx_map = np.clip(idx_map, 0, len(E_neq) - 2)
+    # idx_map = np.searchsorted(y_neq, y_eq) - 1
+    idx_map = np.searchsorted(y_neq, y_eq, side='right') - 1
+    idx_map = np.clip(idx_map, 0, len(y_neq) - 2)
     return idx_map.astype(np.int32)
 
 
-def prepare_interpolation_arrays(T_array: np.ndarray, E_array: np.ndarray, verbose=False) -> dict:
+def prepare_interpolation_arrays(x_array: np.ndarray, y_array: np.ndarray, verbose=False) -> dict:
     """
     Validates input arrays and prepares data for interpolation.
 
     Args:
-        T_array: Array of temperature values
-        E_array: Array of energy density values corresponding to temperatures
+        x_array: Array of x values
+        y_array: Array of y values corresponding to x
         verbose: If True, prints diagnostic information during processing
 
     Returns:
         A dictionary with arrays and metadata for the appropriate interpolation method:
-        - Common keys: 'T_bs', 'E_bs', 'method', 'is_equidistant', 'increment'
-        - Additional keys for double_lookup method: 'T_eq', 'E_neq', 'E_eq',
-          'inv_delta_E_eq', 'idx_map'
+        - Common keys: 'x_bs', 'y_bs', 'method', 'is_equidistant', 'increment'
+        - Additional keys for double_lookup method: 'x_eq', 'y_neq', 'y_eq',
+          'inv_delta_y_eq', 'idx_map'
 
     Raises:
         ValueError: If arrays don't meet requirements for interpolation
     """
     # Convert to numpy arrays if not already
-    T_array = np.asarray(T_array)
-    E_array = np.asarray(E_array)
+    x_array = np.asarray(x_array)
+    y_array = np.asarray(y_array)
 
     # Basic validation
-    if len(T_array) != len(E_array):
-        raise ValueError(f"Energy density array (length {len(E_array)}) and temperature array (length {len(T_array)}) must have the same length!")
+    if len(x_array) != len(y_array):
+        raise ValueError(f"y array (length {len(y_array)}) and x array (length {len(x_array)}) must have the same length!")
 
-    if not E_array.size or not T_array.size:
-        raise ValueError("Energy density and temperature arrays must not be empty!")
+    if not y_array.size or not x_array.size:
+        raise ValueError("x and y arrays must not be empty!")
 
-    if len(E_array) < 2:
-        raise ValueError(f"Arrays must contain at least 2 elements, got {len(E_array)}!")
+    if len(y_array) < 2:
+        raise ValueError(f"Arrays must contain at least 2 elements, got {len(y_array)}!")
 
     # Check if arrays are monotonic
-    T_increasing = np.all(np.diff(T_array) > 0)
-    T_decreasing = np.all(np.diff(T_array) < 0)
-    E_increasing = np.all(np.diff(E_array) > 0)
-    E_decreasing = np.all(np.diff(E_array) < 0)
+    x_increasing = np.all(np.diff(x_array) > 0)
+    x_decreasing = np.all(np.diff(x_array) < 0)
+    y_increasing = np.all(np.diff(y_array) > 0)
+    y_decreasing = np.all(np.diff(y_array) < 0)
 
-    if not (T_increasing or T_decreasing):
+    if not (x_increasing or x_decreasing):
         raise ValueError("Temperature array must be monotonically increasing or decreasing")
 
-    if not (E_increasing or E_decreasing):
+    if not (y_increasing or y_decreasing):
         raise ValueError("Energy density array must be monotonically increasing or decreasing")
 
     # Check if energy increases with temperature (both increasing or both decreasing)
-    if (T_increasing and E_decreasing) or (T_decreasing and E_increasing):
+    if (x_increasing and y_decreasing) or (x_decreasing and y_increasing):
         raise ValueError("Energy density must increase with temperature")
 
     # Create copies to avoid modifying original arrays
-    T_bs = np.copy(T_array)
-    E_bs = np.copy(E_array)
+    x_bs = np.copy(x_array)
+    y_bs = np.copy(y_array)
 
     # Flip arrays if temperature is in descending order
-    if T_decreasing:
+    if x_decreasing:
         if verbose:
             print("Temperature array is descending, flipping arrays for processing")
-        T_bs = np.flip(T_bs)
-        E_bs = np.flip(E_bs)
+        x_bs = np.flip(x_bs)
+        y_bs = np.flip(y_bs)
 
     # Check for strictly increasing values
     has_warnings = False
     try:
-        check_strictly_increasing(T_bs, "Temperature array")
-        check_strictly_increasing(E_bs, "Energy density array")
+        check_strictly_increasing(x_bs, "x array")
+        check_strictly_increasing(y_bs, "y array")
     except ValueError as e:
         has_warnings = True
         if verbose:
@@ -368,34 +367,34 @@ def prepare_interpolation_arrays(T_array: np.ndarray, E_array: np.ndarray, verbo
             print("Continuing with interpolation, but results may be less accurate")
 
     # Use the existing check_equidistant function to determine if suitable for double lookup
-    T_incr = check_equidistant(T_bs)
-    is_equidistant = T_incr != 0.0
+    x_incr = check_equidistant(x_bs)
+    is_equidistant = x_incr != 0.0
 
     # Initialize result with common fields
     result = {
-        "T_bs": T_bs,
-        "E_bs": E_bs,
+        "x_bs": x_bs,
+        "y_bs": y_bs,
         "method": "binary_search",
         "is_equidistant": is_equidistant,
-        "increment": T_incr if is_equidistant else 0.0,
+        "increment": x_incr if is_equidistant else 0.0,
         "has_warnings": has_warnings
     }
 
     # If temperature is equidistant, prepare for double lookup
     if is_equidistant:
         if verbose:
-            print(f"Temperature array is equidistant with increment {T_incr}, using double lookup")
+            print(f"x array is equidistant with increment {x_incr}, using double lookup")
 
         try:
             # Create equidistant energy array and mapping
-            E_eq, inv_delta_E_eq = E_eq_from_E_neq(E_bs)
-            idx_mapping = create_idx_mapping(E_bs, E_eq)
+            y_eq, inv_delta_y_eq = y_eq_from_y_neq(y_bs)
+            idx_mapping = create_idx_mapping(y_bs, y_eq)
 
             result.update({
-                "T_eq": T_bs,
-                "E_neq": E_bs,
-                "E_eq": E_eq,
-                "inv_delta_E_eq": inv_delta_E_eq,
+                "x_eq": x_bs,
+                "y_neq": y_bs,
+                "y_eq": y_eq,
+                "inv_delta_y_eq": inv_delta_y_eq,
                 "idx_map": idx_mapping,
                 "method": "double_lookup"
             })
@@ -404,7 +403,7 @@ def prepare_interpolation_arrays(T_array: np.ndarray, E_array: np.ndarray, verbo
                 print(f"Warning: Failed to create double lookup tables: {e}")
                 print("Falling back to binary search method")
     elif verbose:
-        print("Temperature array is not equidistant, using binary search")
+        print("x array is not equidistant, using binary search")
 
     return result
 
diff --git a/src/pymatlib/core/models.py b/src/pymatlib/core/models.py
index d45512b86e6a005dbe3a285cbe6ef1f0b66dabc1..6f3e60503f58daf78916d695694d52b28a4ebaf6 100644
--- a/src/pymatlib/core/models.py
+++ b/src/pymatlib/core/models.py
@@ -52,14 +52,17 @@ def _prepare_material_expressions(*properties):
     """Prepare expressions and collect assignments from material properties."""
     sub_assignments = []
     expressions = []
-
     for prop in properties:
         if isinstance(prop, MaterialProperty):
             sub_assignments.extend(prop.assignments)
             expressions.append(prop.expr)
         else:
             expressions.append(sympy_wrapper(prop))
-    return expressions, sub_assignments
+    # If there's only one expression, return it directly instead of in a list
+    if len(expressions) == 1:
+        return expressions[0], sub_assignments
+    else:
+        return expressions, sub_assignments
 
 
 def density_by_thermal_expansion(
@@ -88,10 +91,10 @@ def density_by_thermal_expansion(
     if isinstance(temperature, ArrayTypes):
         raise TypeError(f"Incompatible input type for temperature. Expected float or sp.Expr, got {type(temperature)}")
     try:
-        tec_expr = thermal_expansion_coefficient.expr if isinstance(thermal_expansion_coefficient, MaterialProperty) else sympy_wrapper(thermal_expansion_coefficient)
-        sub_assignments = thermal_expansion_coefficient.assignments if isinstance(thermal_expansion_coefficient, MaterialProperty) else []
-        density = density_base * (1 + tec_expr * (temperature - temperature_base)) ** (-3)
-        return MaterialProperty(density, sub_assignments)
+        tec_expr, sub_assignments \
+            = _prepare_material_expressions(thermal_expansion_coefficient)
+        density_expr = density_base * (1 + tec_expr * (temperature - temperature_base)) ** (-3)
+        return MaterialProperty(density_expr, sub_assignments)
     except ZeroDivisionError:
         raise ValueError("Division by zero encountered in density calculation")
 
@@ -123,42 +126,32 @@ def thermal_diffusivity_by_heat_conductivity(
         raise ValueError("Division by zero encountered in thermal diffusivity calculation")
 
 
-def energy_density_standard(
+def specific_enthalpy_sensible(
         temperature: Union[float, sp.Symbol],
-        density: Union[float, MaterialProperty],
-        heat_capacity: Union[float, MaterialProperty],
-        latent_heat: Union[float, MaterialProperty]) \
+        heat_capacity: Union[float, MaterialProperty]) \
         -> MaterialProperty:
 
-    (density_expr, heat_capacity_expr, latent_heat_expr), sub_assignments \
-        = _prepare_material_expressions(density, heat_capacity, latent_heat)
-
-    density_expr = density.expr if isinstance(density, MaterialProperty) else sympy_wrapper(density)
-    heat_capacity_expr = heat_capacity.expr if isinstance(heat_capacity, MaterialProperty) else sympy_wrapper(heat_capacity)
-    latent_heat_expr = latent_heat.expr if isinstance(latent_heat, MaterialProperty) else sympy_wrapper(latent_heat)
+    heat_capacity_expr, sub_assignments \
+        = _prepare_material_expressions(heat_capacity)
 
-    energy_density_expr = density_expr * (temperature * heat_capacity_expr + latent_heat_expr)
-    return MaterialProperty(energy_density_expr, sub_assignments)
+    specific_enthalpy_expr = temperature * heat_capacity_expr
+    return MaterialProperty(specific_enthalpy_expr, sub_assignments)
 
 
-# For backward compatibility
-energy_density = energy_density_standard
-
-
-def energy_density_enthalpy_based(
-        density: Union[float, MaterialProperty],
-        specific_enthalpy: Union[float, MaterialProperty],
+def specific_enthalpy_with_latent_heat(
+        temperature: Union[float, sp.Symbol],
+        heat_capacity: Union[float, MaterialProperty],
         latent_heat: Union[float, MaterialProperty]) \
         -> MaterialProperty:
 
-    (density_expr, specific_enthalpy_expr, latent_heat_expr), sub_assignments \
-        = _prepare_material_expressions(density, specific_enthalpy, latent_heat)
+    (heat_capacity_expr, latent_heat_expr), sub_assignments \
+        = _prepare_material_expressions(heat_capacity, latent_heat)
 
-    energy_density_expr = density_expr * (specific_enthalpy_expr + latent_heat_expr)
-    return MaterialProperty(energy_density_expr, sub_assignments)
+    specific_enthalpy_expr = temperature * heat_capacity_expr + latent_heat_expr
+    return MaterialProperty(specific_enthalpy_expr, sub_assignments)
 
 
-def energy_density_total_enthalpy(
+def energy_density(
         density: Union[float, MaterialProperty],
         specific_enthalpy: Union[float, MaterialProperty]) \
         -> MaterialProperty:
diff --git a/src/pymatlib/core/property_array_extractor.py b/src/pymatlib/core/property_array_extractor.py
new file mode 100644
index 0000000000000000000000000000000000000000..a105518ecb83b8e48d4d2ca2c66a49052e5afa66
--- /dev/null
+++ b/src/pymatlib/core/property_array_extractor.py
@@ -0,0 +1,62 @@
+import numpy as np
+import sympy as sp
+from dataclasses import dataclass, field
+from pymatlib.core.alloy import Alloy
+from pymatlib.core.typedefs import MaterialProperty
+
+
+@dataclass
+class PropertyArrayExtractor:
+    """
+    Extracts arrays of property values from an Alloy at specified temperatures.
+    Attributes:
+        alloy (Alloy): The alloy object containing material properties.
+        temperature_array (np.ndarray): Array of temperature values to evaluate properties at.
+        symbol (sp.Symbol): Symbol to use for property evaluation (e.g., u.center()).
+        specific_enthalpy_array (np.ndarray): Extracted specific enthalpy values.
+        energy_density_array (np.ndarray): Extracted energy density values.
+    """
+    alloy: Alloy
+    temperature_array: np.ndarray
+    symbol: sp.Symbol
+
+    # Arrays will be populated during extraction
+    specific_enthalpy_array: np.ndarray = field(default_factory=lambda: np.array([]))
+    energy_density_array: np.ndarray = field(default_factory=lambda: np.array([]))
+
+    def __post_init__(self):
+        """Initialize arrays after instance creation."""
+        # Extract property arrays after initialization if temperature array is provided.
+        if len(self.temperature_array) >= 2:
+            self.extract_all_arrays()
+
+    def extract_property_array(self, property_name: str) -> np.ndarray:
+        """
+        Extract array of property values at each temperature point using the provided symbol.
+        Args:
+            property_name (str): Name of the property to extract.
+        Returns:
+            np.ndarray: Array of property values corresponding to temperature_array.
+        Raises:
+            ValueError: If property doesn't exist or isn't a MaterialProperty.
+        """
+        # Get the property from the alloy
+        property_obj = getattr(self.alloy, property_name, None)
+        # Check if property exists
+        if property_obj is None:
+            raise ValueError(f"Property '{property_name}' not found in alloy")
+        # Check if it's a MaterialProperty
+        if isinstance(property_obj, MaterialProperty):
+            # Use the symbolic temperature variable from the MaterialProperty
+            return property_obj.evalf(self.symbol, self.temperature_array)
+        else:
+            raise ValueError(f"Property '{property_name}' is not a MaterialProperty")
+
+    def extract_all_arrays(self) -> None:
+        """Extract arrays for all supported properties."""
+        # Extract specific enthalpy array if available
+        if hasattr(self.alloy, 'specific_enthalpy') and self.alloy.specific_enthalpy is not None:
+            self.specific_enthalpy_array = self.extract_property_array('specific_enthalpy')
+        # Extract energy density array if available
+        if hasattr(self.alloy, 'energy_density') and self.alloy.energy_density is not None:
+            self.energy_density_array = self.extract_property_array('energy_density')
diff --git a/src/pymatlib/core/yaml_parser.py b/src/pymatlib/core/yaml_parser.py
index c00efd123322f3e68e8f41f4e58172ae827922c0..dccaa12a742b7da99f6c0ecab52d86fddf520e37 100644
--- a/src/pymatlib/core/yaml_parser.py
+++ b/src/pymatlib/core/yaml_parser.py
@@ -8,7 +8,8 @@ from pymatlib.core.interpolators import interpolate_property
 from pymatlib.core.data_handler import read_data_from_file
 from pymatlib.core.models import (density_by_thermal_expansion,
                                   thermal_diffusivity_by_heat_conductivity,
-                                  energy_density_standard, energy_density_enthalpy_based, energy_density_total_enthalpy)
+                                  specific_enthalpy_sensible, specific_enthalpy_with_latent_heat,
+                                  energy_density)
 from pymatlib.core.typedefs import MaterialProperty
 from ruamel.yaml import YAML, constructor, scanner
 from difflib import get_close_matches
@@ -20,7 +21,6 @@ class PropertyType(Enum):
     FILE = auto()
     KEY_VAL = auto()
     COMPUTE = auto()
-    TUPLE_STRING = auto()
     INVALID = auto()
 
 
@@ -35,16 +35,12 @@ class MaterialConfigParser:
     ABSOLUTE_ZERO = 0.0  # Kelvin
 
     # Define valid properties as class-level constants
-    VALID_PROPERTIES = {
+    VALID_YAML_PROPERTIES = {
         'base_temperature',
         'base_density',
         'density',
         'dynamic_viscosity',
         'energy_density',
-        'energy_density_solidus',
-        'energy_density_liquidus',
-        'energy_density_temperature_array',
-        'energy_density_array',
         'heat_capacity',
         'heat_conductivity',
         'kinematic_viscosity',
@@ -52,8 +48,6 @@ class MaterialConfigParser:
         'latent_heat_of_vaporization',
         'specific_enthalpy',
         'surface_tension',
-        'temperature',
-        'temperature_array',
         'thermal_diffusivity',
         'thermal_expansion_coefficient',
     }
@@ -76,6 +70,7 @@ class MaterialConfigParser:
         self.base_dir = Path(yaml_path).parent
         self.config: Dict[str, Any] = self._load_yaml()
         self._validate_config()
+        self._process_temperature_range(self.config['temperature_range'])
 
     def _load_yaml(self) -> Dict[str, Any]:
         """Load and parse YAML file.
@@ -115,14 +110,34 @@ class MaterialConfigParser:
         if not isinstance(self.config, dict):
             # raise ValueError("Root YAML element must be a mapping")
             raise ValueError("The YAML file must start with a dictionary/object structure with key-value pairs, not a list or scalar value")
-        if 'properties' not in self.config:
-            raise ValueError("Missing 'properties' section in configuration")
+        self._validate_required_fields()
         properties = self.config.get('properties', {})
         if not isinstance(properties, dict):
-            # raise ValueError("'properties' must be a mapping")
             raise ValueError("The 'properties' section in your YAML file must be a dictionary with key-value pairs")
         self._validate_property_names(properties)
-        self._validate_required_fields()
+
+    def _validate_required_fields(self) -> None:
+        """
+        Validate required configuration fields.
+        Raises:
+            ValueError: If any required field is missing.
+        """
+        required_fields = {'name', 'composition', 'solidus_temperature', 'liquidus_temperature', 'temperature_range', 'properties'}
+        missing_fields = required_fields - set(self.config.keys())
+        if missing_fields:
+            raise ValueError(f"Missing required fields: {', '.join(missing_fields)}")
+        # Check for extra fields
+        extra_fields = set(self.config.keys()) - required_fields
+        if extra_fields:
+            suggestions = {
+                field: get_close_matches(field, required_fields, n=1, cutoff=0.6)
+                for field in extra_fields
+            }
+            error_msg = "Extra fields found in configuration:\n"
+            for field, matches in suggestions.items():
+                suggestion = f" (did you mean '{matches[0]}'?)" if matches else ""
+                error_msg += f"  - '{field}'{suggestion}\n"
+            raise ValueError(error_msg)
 
     def _validate_property_names(self, properties: Dict[str, Any]) -> None:
         """
@@ -132,10 +147,10 @@ class MaterialConfigParser:
         Raises:
             ValueError: If any property name is not in VALID_PROPERTIES.
         """
-        invalid_props = set(properties.keys()) - self.VALID_PROPERTIES
+        invalid_props = set(properties.keys()) - self.VALID_YAML_PROPERTIES
         if invalid_props:
             suggestions = {
-                prop: get_close_matches(prop, self.VALID_PROPERTIES, n=1, cutoff=0.6)
+                prop: get_close_matches(prop, self.VALID_YAML_PROPERTIES, n=1, cutoff=0.6)
                 for prop in invalid_props
             }
             error_msg = "Invalid properties found:\n"
@@ -144,16 +159,66 @@ class MaterialConfigParser:
                 error_msg += f"  - '{prop}'{suggestion}\n"
             raise ValueError(error_msg)
 
-    def _validate_required_fields(self) -> None:
+    ##################################################
+    # Temperature Array Processing
+    ##################################################
+
+    def _process_temperature_range(self, array_def: List[Union[float, int]]) -> np.ndarray:
         """
-        Validate required configuration fields.
+        Process temperature array definition with format [start, end, points/delta].
+        Args:
+            array_def (List[Union[float, int]]): A list defining the temperature array in the format
+                                                 [start, end, points/delta].
+        Returns:
+            np.ndarray: An array of temperature values.
         Raises:
-            ValueError: If any required field is missing.
+            ValueError: If the array definition is invalid, improperly formatted,
+                        or contains physically impossible temperatures.
+        Examples:
+            >>> self._process_temperature_range([300, 3000, 5])
+            array([ 300., 975., 1650., 2325., 3000.])
+            >>> self._process_temperature_range([3000, 300, -1350.])
+            array([3000., 1650., 300.])
         """
-        required_fields = {'name', 'composition', 'solidus_temperature', 'liquidus_temperature'}
-        missing_fields = required_fields - set(self.config.keys())
-        if missing_fields:
-            raise ValueError(f"Missing required fields: {', '.join(missing_fields)}")
+        if not (isinstance(array_def, list) and len(array_def) == 3):
+            raise ValueError("Temperature array must be defined as [start, end, points/delta]")
+        try:
+            start, end, step = float(array_def[0]), float(array_def[1]), array_def[2]
+            if start <= self.ABSOLUTE_ZERO or end <= self.ABSOLUTE_ZERO:
+                raise ValueError(f"Temperature must be above absolute zero ({self.ABSOLUTE_ZERO}K)")
+            if abs(float(step)) < self.EPSILON:
+                raise ValueError("Delta or number of points cannot be zero.")
+            # Check if step represents delta (float) or points (int)
+            if isinstance(step, float):
+                temperature_array = self._process_float_step(start, end, step)
+            else:
+                temperature_array = self._process_int_step(start, end, int(step))
+            self.temperature_array = temperature_array  # Store the temperature array
+            # print(self.temperature_array)
+            return temperature_array
+        except ValueError as e:
+            raise ValueError(f"Invalid temperature array definition \n -> {e}")
+
+    @staticmethod
+    def _process_float_step(start: float, end: float, delta: float) -> np.ndarray:
+        """Process temperature array with float step (delta)."""
+        if start < end and delta <= 0:
+            raise ValueError("Delta must be positive for increasing range")
+        if start > end and delta >= 0:
+            raise ValueError("Delta must be negative for decreasing range")
+        max_delta = abs(end - start)
+        if abs(delta) > max_delta:
+            raise ValueError(f"Absolute value of delta ({abs(delta)}) is too large for the range. It should be <= {max_delta}")
+        return np.arange(start, end + delta/2, delta)
+
+    def _process_int_step(self, start: float, end: float, points: int) -> np.ndarray:
+        """Process temperature array with integer step (number of points)."""
+        print(f"Processing TA as int: start={start}, end={end}, points={points}")
+        if points <= 0:
+            raise ValueError(f"Number of points must be positive, got {points}!")
+        if points < self.MIN_POINTS:
+            raise ValueError(f"Number of points must be at least {self.MIN_POINTS}, got {points}!")
+        return np.linspace(start, end, points)
 
     @staticmethod
     def _validate_property_value(prop: str, value: Union[float, np.ndarray]) -> None:
@@ -168,12 +233,11 @@ class MaterialConfigParser:
         """
         # Define property constraints with SI units
         BASE_PROPERTIES = {'base_temperature', 'base_density'}  # always single float values
-        POSITIVE_PROPERTIES = {'density', 'heat_capacity', 'heat_conductivity', 'temperature',
+        POSITIVE_PROPERTIES = {'density', 'heat_capacity', 'heat_conductivity',
                                'dynamic_viscosity', 'kinematic_viscosity', 'thermal_diffusivity',
                                'surface_tension'}
-        NON_NEGATIVE_PROPERTIES = {'latent_heat_of_fusion', 'latent_heat_of_vaporization', 'energy_density',
-                                   'energy_density_solidus', 'energy_density_liquidus', 'energy_density_array'}
-        ARRAY_PROPERTIES = {'temperature_array', 'energy_density_temperature_array', 'energy_density_array'}
+        NON_NEGATIVE_PROPERTIES = {'latent_heat_of_fusion', 'latent_heat_of_vaporization', 'energy_density',}
+                                   # 'energy_density_solidus', 'energy_density_liquidus'}
 
         # Property-specific range constraints
         PROPERTY_RANGES = {
@@ -182,10 +246,6 @@ class MaterialConfigParser:
             'density': (100, 22000),  # kg/m³
             'dynamic_viscosity': (1e-4, 1e5),  # Pa·s
             'energy_density': (0, 1e8),  # J/m³
-            'energy_density_solidus': (0, 1e8),  # J/m³
-            'energy_density_liquidus': (0, 1e8),  # J/m³
-            'energy_density_temperature_array': (0, 5000),  # K
-            'energy_density_array': (0, 1e8),  # J/m³
             'heat_capacity': (100, 10000),  # J/(kg·K)
             'heat_conductivity': (1, 600),  # W/(m·K)
             'kinematic_viscosity': (1e-8, 1e-3),  # m²/s
@@ -193,8 +253,6 @@ class MaterialConfigParser:
             'latent_heat_of_vaporization': (50000, 12000000),  # J/kg
             'specific_enthalpy': (0, 15000000),  # J/kg
             'surface_tension': (0.1, 3.0),  # N/m
-            'temperature': (0, 5000),  # K
-            'temperature_array': (0, 5000),  # K
             'thermal_diffusivity': (1e-8, 1e-3),  # m²/s
             'thermal_expansion_coefficient': (-3e-5, 3e-5),  # 1/K
         }
@@ -213,7 +271,7 @@ class MaterialConfigParser:
                         bad_values = value[value <= 0]
                         raise ValueError(f"All '{prop}' values must be positive. Found {len(bad_indices)} invalid values "
                                          f"at indices {bad_indices}: {bad_values}.")
-                if prop in NON_NEGATIVE_PROPERTIES or prop in ARRAY_PROPERTIES:
+                if prop in NON_NEGATIVE_PROPERTIES:
                     if (value < 0).any():
                         bad_indices = np.where(value < 0)[0]
                         bad_values = value[value < 0]
@@ -253,6 +311,37 @@ class MaterialConfigParser:
         except Exception as e:
             raise ValueError(f"Failed to validate property value \n -> {e}")
 
+    def _validate_temperature_range(self, prop: str, temp_array: np.ndarray) -> None:
+        """
+        Validate that temperature arrays from properties fall within the global temperature range.
+        Args:
+            prop (str): The name of the property being validated
+            temp_array (np.ndarray): The temperature array to validate
+        Raises:
+            ValueError: If temperature values fall outside the global temperature range
+        """
+        # Extract global temperature range from config
+        # temp_range = self.config.get('temperature_range', [])
+        temp_range = self.temperature_array
+        min_temp, max_temp = np.min(temp_range), np.max(temp_range)
+        # print(min_temp, max_temp)
+        # Check for temperatures outside the global range
+        if ((temp_array < min_temp) | (temp_array > max_temp)).any():
+            out_of_range = np.where((temp_array < min_temp) | (temp_array > max_temp))[0]
+            out_values = temp_array[out_of_range]
+            # Create a more detailed error message
+            if len(out_of_range) > 5:
+                # Show just a few examples if there are many out-of-range values
+                sample_indices = out_of_range[:5]
+                sample_values = temp_array[sample_indices]
+                raise ValueError(f"Property '{prop}' contains temperature values outside global range [{min_temp}, {max_temp}] "
+                                 f"\n -> Found {len(out_of_range)} out-of-range values, first 5 at indices {sample_indices}: {sample_values}"
+                                 f"\n -> Min value: {temp_array.min()}, Max value: {temp_array.max()}")
+            else:
+                # Show all out-of-range values if there aren't too many
+                raise ValueError(f"Property '{prop}' contains temperature values outside global range [{min_temp}, {max_temp}] "
+                                 f"\n -> Found {len(out_of_range)} out-of-range values at indices {out_of_range}: {out_values}")
+
     ##################################################
     # Alloy Creation
     ##################################################
@@ -425,8 +514,6 @@ class MaterialConfigParser:
                 return PropertyType.KEY_VAL
             elif self._is_compute_property(config):
                 return PropertyType.COMPUTE
-            elif prop_name == 'energy_density_temperature_array' and isinstance(config, str) and config.startswith('(') and config.endswith(')'):
-                return PropertyType.TUPLE_STRING
             else:
                 return PropertyType.INVALID
         except Exception as e:
@@ -447,7 +534,6 @@ class MaterialConfigParser:
             PropertyType.FILE: [],
             PropertyType.KEY_VAL: [],
             PropertyType.COMPUTE: [],
-            PropertyType.TUPLE_STRING: []
         }
         for prop_name, config in properties.items():
             try:
@@ -486,9 +572,6 @@ class MaterialConfigParser:
                         self._process_key_val_property(alloy, prop_name, config, T)
                     elif prop_type == PropertyType.COMPUTE:
                         self._process_computed_property(alloy, prop_name, T)
-                    elif prop_type == PropertyType.TUPLE_STRING:
-                        # Handle tuple string properties if needed
-                        pass
         except Exception as e:
             raise ValueError(f"Failed to process properties \n -> {e}")
 
@@ -538,6 +621,8 @@ class MaterialConfigParser:
                 # For string configuration, construct the full path
                 file_path = str(yaml_dir / file_config)
                 temp_array, prop_array = read_data_from_file(file_path)
+            # Validate the temperature array
+            self._validate_temperature_range(prop_name, temp_array)
             # Validate the property array
             self._validate_property_value(prop_name, prop_array)
             self._process_property_data(alloy, prop_name, T, temp_array, prop_array)
@@ -563,6 +648,8 @@ class MaterialConfigParser:
             val_array = np.array(prop_config['val'], dtype=float)
             if len(key_array) != len(val_array):
                 raise ValueError(f"Length mismatch in {prop_name}: key and val arrays must have same length")
+            # Validate the temperature array
+            self._validate_temperature_range(prop_name, key_array)
             # Validate the value array
             self._validate_property_value(prop_name, val_array)
             self._process_property_data(alloy, prop_name, T, key_array, val_array)
@@ -665,7 +752,7 @@ class MaterialConfigParser:
         try:
             if isinstance(T, sp.Symbol):
                 self._process_symbolic_temperature(alloy, prop_name, T, temp_array, prop_array)
-            elif isinstance(T, (float, int)):
+            elif isinstance(T, float):
                 self._process_constant_temperature(alloy, prop_name, T, temp_array, prop_array)
             else:
                 raise ValueError(f"Unexpected type for T: {type(T)}")
@@ -683,13 +770,10 @@ class MaterialConfigParser:
             temp_array (np.ndarray): Array of temperature values.
             prop_array (np.ndarray): Array of property values.
         """
-        # If T is symbolic, store the full temperature array if not already set then interpolate
-        if getattr(alloy, 'temperature_array', None) is None or len(alloy.temperature_array) == 0:
-            alloy.temperature_array = temp_array
         material_property = interpolate_property(T, temp_array, prop_array)
         setattr(alloy, prop_name, material_property)
         if prop_name == 'energy_density':
-            self._process_energy_density(alloy, material_property, T, temp_array, prop_array)
+            self._process_energy_density(alloy, material_property, T)
 
     @staticmethod
     def _process_constant_temperature(alloy: Alloy, prop_name: str, T: Union[float, int], temp_array: np.ndarray, prop_array: np.ndarray) -> None:
@@ -709,18 +793,16 @@ class MaterialConfigParser:
         setattr(alloy, prop_name, material_property)
 
     @staticmethod
-    def _process_energy_density(alloy: Alloy, material_property: Any, T: sp.Symbol, temp_array: np.ndarray, prop_array: np.ndarray) -> None:
+    def _process_energy_density(alloy: Alloy, material_property: Any, T: sp.Symbol) -> None:
         """
         Process additional properties for energy density.
         Args:
             alloy (Alloy): The alloy object to update.
             material_property (Any): The interpolated material property.
             T (sp.Symbol): Symbolic temperature.
-            temp_array (np.ndarray): Array of temperature values.
-            prop_array (np.ndarray): Array of property values.
+            # temp_array (np.ndarray): Array of temperature values.
+            # prop_array (np.ndarray): Array of property values.
         """
-        alloy.energy_density_temperature_array = temp_array
-        alloy.energy_density_array = prop_array
         alloy.energy_density_solidus = material_property.evalf(T, alloy.temperature_solidus)
         alloy.energy_density_liquidus = material_property.evalf(T, alloy.temperature_liquidus)
 
@@ -761,7 +843,7 @@ class MaterialConfigParser:
         setattr(alloy, prop_name, material_property)
         # Handle special case for energy_density
         if prop_name == 'energy_density' and isinstance(T, sp.Symbol):
-            self._handle_energy_density(alloy, material_property, T, method_dependencies)
+            self._handle_energy_density(alloy, material_property, T)
 
     @staticmethod
     def _get_computation_methods(alloy: Alloy, T: Union[float, sp.Symbol]):
@@ -789,21 +871,21 @@ class MaterialConfigParser:
                     alloy.heat_capacity
                 )
             },
-            'energy_density': {
-                'default': lambda: energy_density_standard(
+            'specific_enthalpy': {
+                'default': lambda: specific_enthalpy_sensible(
+                    T,
+                    alloy.heat_capacity
+                ),
+                'latent_heat_based': lambda: specific_enthalpy_with_latent_heat(
                     T,
-                    alloy.density,
                     alloy.heat_capacity,
                     alloy.latent_heat_of_fusion
-                ),
-                'enthalpy_based': lambda: energy_density_enthalpy_based(
+                )
+            },
+            'energy_density': {
+                'default': lambda: energy_density(
                     alloy.density,
                     alloy.specific_enthalpy,
-                    alloy.latent_heat_of_fusion
-                ),
-                'total_enthalpy': lambda: energy_density_total_enthalpy(
-                    alloy.density,
-                    alloy.specific_enthalpy
                 ),
             },
         }
@@ -823,10 +905,12 @@ class MaterialConfigParser:
             'thermal_diffusivity': {
                 'default': ['heat_conductivity', 'density', 'heat_capacity'],
             },
+            'specific_enthalpy': {
+                'default': ['heat_capacity'],
+                'latent_heat_based': ['heat_capacity', 'latent_heat_of_fusion'],
+            },
             'energy_density': {
-                'default': ['density', 'heat_capacity', 'latent_heat_of_fusion'],
-                'enthalpy_based': ['density', 'specific_enthalpy', 'latent_heat_of_fusion'],
-                'total_enthalpy': ['density', 'specific_enthalpy'],
+                'default': ['density', 'specific_enthalpy'],
             },
         }
 
@@ -854,7 +938,7 @@ class MaterialConfigParser:
         if missing_deps:
             raise ValueError(f"Cannot compute {prop_name}. Missing dependencies: {missing_deps}")
 
-    def _handle_energy_density(self, alloy: Alloy, material_property: MaterialProperty, T: sp.Symbol, dependencies: List[str]):
+    def _handle_energy_density(self, alloy: Alloy, material_property: MaterialProperty, T: sp.Symbol):
         """
         Handle the special case of energy density computation.
         This method computes additional properties related to energy density when T is symbolic.
@@ -863,92 +947,32 @@ class MaterialConfigParser:
             alloy (Alloy): The alloy object to process.
             material_property (MaterialProperty): The computed energy density property.
             T (sp.Symbol): The symbolic temperature variable.
-            dependencies (List[str]): List of dependencies for energy density computation.
+            # dependencies (List[str]): List of dependencies for energy density computation.
         Raises:
             ValueError: If T is not symbolic or if energy_density_temperature_array is not defined in the config.
         """
         # Ensure T is symbolic
         if not isinstance(T, sp.Symbol):
             raise ValueError("_handle_energy_density should only be called with symbolic T")
-        # Check dependencies
-        deps_to_check = [getattr(alloy, dep) for dep in dependencies if hasattr(alloy, dep)]
-        if any(isinstance(dep, MaterialProperty) for dep in deps_to_check):
-            if 'energy_density_temperature_array' not in self.config['properties']:
-                raise ValueError(f"energy_density_temperature_array must be defined when energy_density is computed with symbolic T")
-            # Process energy_density_temperature_array
-            edta = self.config['properties']['energy_density_temperature_array']
-            alloy.energy_density_temperature_array = self._process_edta(edta)
-        if len(alloy.energy_density_temperature_array) >= 2:
-            alloy.energy_density_array = material_property.evalf(T, alloy.energy_density_temperature_array)
+        if len(self.temperature_array) >= 2:
+            # alloy.energy_density_array = material_property.evalf(T, self.temperature_array)
             alloy.energy_density_solidus = material_property.evalf(T, alloy.temperature_solidus)
             alloy.energy_density_liquidus = material_property.evalf(T, alloy.temperature_liquidus)
 
-    ##################################################
-    # Energy Density Temperature Array Processing
-    ##################################################
-
-    def _process_edta(self, array_def: str) -> np.ndarray:
-        """
-        Process temperature array definition with format (start, end, points/delta).
-        Args:
-            array_def (str): A string defining the temperature array in the format
-                             "(start, end, points/delta)".
-        Returns:
-            np.ndarray: An array of temperature values.
-        Raises:
-            ValueError: If the array definition is invalid, improperly formatted,
-                        or contains physically impossible temperatures.
-        Examples:
-            >>> self._process_edta("(300, 3000, 5)")
-            array([ 300., 975., 1650., 2325., 3000.])
-            >>> self._process_edta("(3000, 300, -1350.)")
-            array([3000., 1650., 300.])
-        """
-        if not (isinstance(array_def, str) and array_def.startswith('(') and array_def.endswith(')')):
-            raise ValueError("Temperature array must be defined as (start, end, points/delta)")
-        try:
-            # Parse the tuple string
-            values = [v.strip() for v in array_def.strip('()').split(',')]
-            if len(values) != 3:
-                raise ValueError("'energy_density_temperature_array' must be a tuple of three comma-separated values representing (start, end, points/step)")
-            start, end, step = float(values[0]), float(values[1]), values[2]
-            if start <= self.ABSOLUTE_ZERO or end <= self.ABSOLUTE_ZERO:
-                raise ValueError(f"Temperature must be above absolute zero ({self.ABSOLUTE_ZERO}K)")
-            if abs(float(step)) < self.EPSILON:
-                raise ValueError("Delta or number of points cannot be zero.")
-            # Check if step represents delta (float) or points (int)
-            if '.' in step or 'e' in step.lower():
-                return self._process_float_step(start, end, float(step))
-            else:
-                return self._process_int_step(start, end, int(step))
-        except ValueError as e:
-            raise ValueError(f"Invalid temperature array definition \n -> {e}")
-
-    @staticmethod
-    def _process_float_step(start: float, end: float, delta: float) -> np.ndarray:
-        """Process temperature array with float step (delta)."""
-        if start < end and delta <= 0:
-            raise ValueError("Delta must be positive for increasing range")
-        if start > end and delta >= 0:
-            raise ValueError("Delta must be negative for decreasing range")
-        max_delta = abs(end - start)
-        if abs(delta) > max_delta:
-            raise ValueError(f"Absolute value of delta ({abs(delta)}) is too large for the range. It should be <= {max_delta}")
-        return np.arange(start, end + delta/2, delta)
-
-    def _process_int_step(self, start: float, end: float, points: int) -> np.ndarray:
-        """Process temperature array with integer step (number of points)."""
-        if points <= 0:
-            raise ValueError(f"Number of points must be positive, got {points}!")
-        if points < self.MIN_POINTS:
-            raise ValueError(f"Number of points must be at least {self.MIN_POINTS}, got {points}!")
-        return np.linspace(start, end, points)
-
 ##################################################
 # External Function
 ##################################################
 
 def create_alloy_from_yaml(yaml_path: Union[str, Path], T: Union[float, sp.Symbol]) -> Alloy:
-    """Create alloy instance from YAML configuration file"""
+    """
+    Create alloy instance from YAML configuration file
+    Args:
+        yaml_path: Path to the YAML configuration file
+        T: Temperature value or symbol for property evaluation
+    Returns:
+        Tuple containing the alloy instance and the temperature array
+    """
     parser = MaterialConfigParser(yaml_path)
-    return parser.create_alloy(T)
+    alloy = parser.create_alloy(T)
+    temp_array = parser.temperature_array
+    return alloy, temp_array
diff --git a/src/pymatlib/data/alloys/SS304L/SS304L.py b/src/pymatlib/data/alloys/SS304L/SS304L.py
deleted file mode 100644
index 7483217f8a39de102dbf4b81ec1b4ff67e3bfd79..0000000000000000000000000000000000000000
--- a/src/pymatlib/data/alloys/SS304L/SS304L.py
+++ /dev/null
@@ -1,181 +0,0 @@
-import time
-import numpy as np
-import sympy as sp
-from pathlib import Path
-from typing import Union
-from pymatlib.core.alloy import Alloy
-from pymatlib.data.element_data import Fe, Cr, Ni, Mo, Mn
-from pymatlib.core.models import thermal_diffusivity_by_heat_conductivity, density_by_thermal_expansion, energy_density, energy_density_total_enthalpy
-from pymatlib.core.data_handler import read_data_from_txt, celsius_to_kelvin, thousand_times, read_data_from_excel, read_data_from_file, plot_arrays
-from pymatlib.core.interpolators import interpolate_property, prepare_interpolation_arrays, interpolate_binary_search, interpolate_double_lookup
-
-
-def create_SS304L(T: Union[float, sp.Symbol]) -> Alloy:
-    """
-    Creates an Alloy instance for SS304L stainless steel with specific properties.
-    Args:
-        T (Union[float, sp.Symbol]): Temperature as a symbolic variable or numeric value.
-    Returns:
-        Alloy: Initialized SS304L alloy with physical properties.
-    Notes:
-        - **Material Properties**:
-            - **Density**: 8.0 g/cm³ (8000 kg/m³) at room temperature
-            - **Composition**:
-                - Iron (Fe): 67.5 wt% (Balance)
-                - Chromium (Cr): 17.0 wt%
-                - Nickel (Ni): 12.0 wt%
-                - Molybdenum (Mo): 2.5 wt%
-                - Manganese (Mn): 1.0 wt%
-            - **Phase Transitions**:
-                - Solidus: 1658.0 K (1385°C)
-                - Liquidus: 1723.0 K (1450°C)
-            - **Thermal Expansion**: 16.3 × 10^-6 /K at room temperature
-        - **Data Units**: All input data should be in SI units:
-            - **Temperature**: Kelvin (K)
-            - **Density**: kg/m³
-            - **Heat Capacity**: J/(kg·K)
-            - **Heat Conductivity**: W/(m·K)
-        - **Temperature Range**: Valid from room temperature (273.15K) to 2000K
-        - **Property Variations**: Properties are temperature-dependent and implemented as piecewise functions
-        - **Data Sources**: Property values based on experimental data and literature
-        - **Input Data Files**: Required files in same directory:
-            - density_temperature.txt
-            - heat_capacity_temperature.txt
-            - heat_conductivity_temperature.txt
-    Example:
-        >>> T = sp.Symbol('T')
-        >>> ss316l = create_SS304L(T)
-        >>> density_at_1000K = ss316l.density.evalf(T, 1000.0)
-    """
-    # Define the alloy with specific elemental composition and phase transition temperatures
-    SS304L = Alloy(
-        elements=[Fe, Cr, Ni, Mo, Mn],
-        composition=[0.675, 0.17, 0.12, 0.025, 0.01],  # Fe: 67.5%, Cr: 17%, Ni: 12%, Mo: 2.5%, Mn: 1%
-        temperature_solidus=1605.,  # Solidus temperature in Kelvin (test at 1653.15 K = 1380 C)
-        temperature_liquidus=1735.,  # Liquidus temperature in Kelvin (test at 1723.15 K = 1450 C)
-        thermal_expansion_coefficient=16.3e-6  # in 1/K
-    )
-    # density_data_file_path = "/local/ca00xebo/repos/pymatlib/src/pymatlib/data/alloys/SS304L/density_temperature.txt"
-    # Determine the base directory
-    base_dir = Path(__file__).parent  # Directory of the current file
-
-    # Paths to data files using relative paths
-    # density_data_file_path = str(base_dir / 'density_temperature_edited.txt')
-    density_data_file_path = str(base_dir / '304L_Erstarrungsdaten_edited.xlsx')
-    # heat_capacity_data_file_path = str(base_dir / 'heat_capacity_temperature_edited.txt')
-    heat_capacity_data_file_path = str(base_dir / '304L_Erstarrungsdaten_edited.xlsx')
-    heat_conductivity_data_file_path = str(base_dir / '..' / 'SS304L' / '304L_Erstarrungsdaten_edited.xlsx')
-    latent_heat_of_fusion_data_file_path = str(base_dir / '..' / 'SS304L' / '304L_Erstarrungsdaten_edited.xlsx')
-
-    # Read temperature and material property data from the files
-    # density_temp_array, density_array = read_data_from_txt(density_data_file_path)
-    density_temp_array, density_array = read_data_from_excel(density_data_file_path, temp_col="T (K)", prop_col="Density (kg/(m)^3)")
-    # heat_capacity_temp_array, heat_capacity_array = read_data_from_txt(heat_capacity_data_file_path)
-    heat_capacity_temp_array, heat_capacity_array = read_data_from_excel(heat_capacity_data_file_path, temp_col="T (K)", prop_col="Specific heat regular weighted")
-    heat_conductivity_temp_array, heat_conductivity_array = read_data_from_excel(heat_conductivity_data_file_path, temp_col="T (K)", prop_col="Thermal conductivity (W/(m*K))-TOTAL-10000.0(K/s)")
-    latent_heat_of_fusion_temp_array, latent_heat_of_fusion_array = read_data_from_excel(latent_heat_of_fusion_data_file_path, temp_col="T (K)", prop_col="Latent heat (J/Kg)")
-    _, sp_enthalpy_array = read_data_from_excel(density_data_file_path, temp_col="T (K)", prop_col="Enthalpy (J/g)-TOTAL-10000.0(K/s)")
-    _, sp_enthalpy_weighted = read_data_from_excel(density_data_file_path, temp_col="T (K)", prop_col="Enthalpy (J/kg) weighted")
-    u1_temp_array, u1_array = read_data_from_excel(density_data_file_path, temp_col="T (K)", prop_col="u1 (rho*h)")
-    _, u2 = read_data_from_excel(density_data_file_path, temp_col="T (K)", prop_col="u2 (rho*h + rho*L)")
-    _, u3 = read_data_from_excel(density_data_file_path, temp_col="T (K)", prop_col="u3 (rho*h_weighted + rho*L)")
-
-    # Ensure the data was loaded correctly
-    if any(arr.size == 0 for arr in [density_temp_array, density_array, heat_capacity_temp_array, heat_capacity_array,
-                                     heat_conductivity_temp_array, heat_conductivity_array]):
-        raise ValueError("Failed to load temperature or material data from the file.")
-
-    SS304L.heat_conductivity = interpolate_property(T, heat_conductivity_temp_array, heat_conductivity_array)
-    SS304L.density = interpolate_property(T, density_temp_array, density_array)
-    SS304L.heat_capacity = interpolate_property(T, heat_capacity_temp_array, heat_capacity_array)
-    SS304L.thermal_diffusivity = thermal_diffusivity_by_heat_conductivity(SS304L.heat_conductivity, SS304L.density, SS304L.heat_capacity)
-    # SS304L.latent_heat_of_fusion = interpolate_property(T, SS304L.solidification_interval(), np.array([171401.0, 0.0]))
-    SS304L.latent_heat_of_fusion = interpolate_property(T, latent_heat_of_fusion_temp_array, latent_heat_of_fusion_array)
-    SS304L.specific_enthalpy = interpolate_property(T, density_temp_array, sp_enthalpy_array)
-    SS304L.energy_density = energy_density_total_enthalpy(SS304L.density, SS304L.specific_enthalpy)
-    SS304L.energy_density_solidus = SS304L.energy_density.evalf(T, SS304L.temperature_solidus)
-    SS304L.energy_density_liquidus = SS304L.energy_density.evalf(T, SS304L.temperature_liquidus)
-
-    # Populate temperature_array and energy_density_array
-    SS304L.temperature_array = density_temp_array
-    SS304L.energy_density_temperature_array = u1_temp_array
-
-    SS304L.energy_density_array = np.array([
-        SS304L.energy_density.evalf(T, temp) for temp in density_temp_array
-    ])
-
-    args = (SS304L.temperature_array,
-            SS304L.energy_density_liquidus,
-            SS304L.energy_density_array)
-
-    E = SS304L.energy_density.evalf(T, SS304L.temperature_liquidus)
-    result = prepare_interpolation_arrays(
-        SS304L.temperature_array,
-        SS304L.energy_density_array
-    )
-
-    if result["method"] == "double_lookup":
-        start_time3 = time.perf_counter()
-        T_interpolate3 = interpolate_double_lookup(
-            E,
-            result["T_eq"],
-            result["E_neq"],
-            result["E_eq"],
-            result["inv_delta_E_eq"],
-            result["idx_map"]
-        )
-        execution_time3 = time.perf_counter() - start_time3
-        print(f"Interpolated temperature: {T_interpolate3}")
-        print(f"Execution time for double lookup: {execution_time3:.8f} seconds")
-    elif result["method"] == "binary_search":  # Fixed syntax here
-        start_time3 = time.perf_counter()
-        T_interpolate3 = interpolate_binary_search(  # Changed function name to match method
-            E,
-            result["T_bs"],
-            result["E_bs"]
-        )
-        execution_time3 = time.perf_counter() - start_time3
-        print(f"Interpolated temperature: {T_interpolate3}")
-        print(f"Execution time for binary search: {execution_time3:.8f} seconds")
-    else:
-        print("Unknown interpolation method")
-
-    T_star = interpolate_binary_search(*args)
-    print(f"T_star: {T_star}")
-
-    return SS304L
-
-
-if __name__ == '__main__':
-    Temp = sp.Symbol('T')
-    alloy = create_SS304L(Temp)
-
-    # Print the composition of each element in the alloy
-    for i in range(len(alloy.composition)):
-        print(f"Element {alloy.elements[i]}: {alloy.composition[i]}")
-
-    print("\nTesting SS304L with symbolic temperature:")
-    for field in vars(alloy):
-        print(f"{field} = {alloy.__getattribute__(field)}")
-
-    # Test interpolate_property
-    print("\nTesting interpolate_property:")
-    test_temp = 2003.15  # Example temperature value
-
-    # Interpolate density, heat capacity, and heat conductivity
-    density_result = alloy.density.evalf(Temp, test_temp)
-    print(f"Interpolated density at T={test_temp} using evalf: {density_result}")
-
-    heat_capacity_result = alloy.heat_capacity.evalf(Temp, test_temp)
-    print(f"Interpolated heat capacity at T={test_temp} using evalf: {heat_capacity_result}")
-
-    heat_conductivity_result = alloy.heat_conductivity.evalf(Temp, test_temp)
-    print(f"Interpolated heat conductivity at T={test_temp} using evalf: {heat_conductivity_result}")
-
-    # Test thermal diffusivity calculation
-    heat_conductivity = heat_conductivity_result  # Example value for testing
-    density = density_result  # Example value for testing
-    heat_capacity = heat_capacity_result  # Example value for testing
-    thermal_diffusivity = thermal_diffusivity_by_heat_conductivity(
-        heat_conductivity, density, heat_capacity)
-    print(f"Calculated thermal diffusivity at T={test_temp}: {thermal_diffusivity}")
diff --git a/src/pymatlib/data/alloys/SS304L/SS304L.yaml b/src/pymatlib/data/alloys/SS304L/SS304L.yaml
index e17cfc717e14a8b786d631f8d104f07bf6d44fc3..19bb2124f4d67cf18427080b518748acc6848a78 100644
--- a/src/pymatlib/data/alloys/SS304L/SS304L.yaml
+++ b/src/pymatlib/data/alloys/SS304L/SS304L.yaml
@@ -128,6 +128,9 @@ composition:
 solidus_temperature: 1605.
 liquidus_temperature: 1735.
 
+# temperature_range: [300, 3000, 541]  # 541 (int) -> points
+#OR
+temperature_range: [3000, 300, -5.0]  # 5.0 (float) -> increment
 
 properties:
   # energy_density:
@@ -136,10 +139,9 @@ properties:
   #OR
   # energy_density: compute
   #OR
-  energy_density:
-    compute: total_enthalpy  # Explicit model selection
+  energy_density: compute
   # User can specify either:
-  energy_density_temperature_array: (300, 3000, 541)  # int for number of points
+  # energy_density_temperature_array: (300, 3000, 541)  # int for number of points
   # OR
   #energy_density_temperature_array: (300, 3000, 5.0)   # float for delta (increment)
   # save the energy_density and temperature as arrays always.
@@ -147,14 +149,20 @@ properties:
 
   #density: 7950.
   #OR
-  density: compute  # computed by thermal expansion coefficient, should be acceptable even if TEC is defined later in the file
+  #density: compute  # computed by thermal expansion coefficient, should be acceptable even if TEC is defined later in the file
   base_temperature: 2273.
   base_density: 6.591878918e3
   #OR
+  density:
+    file: ./304L_Erstarrungsdaten_edited.xlsx
+    temp_col: T (K)
+    prop_col: Density (kg/(m)^3)
 
 
-
-
+  heat_conductivity:
+    file: ./304L_Erstarrungsdaten_edited.xlsx
+    temp_col: T (K)
+    prop_col: Thermal conductivity (W/(m*K))-TOTAL-10000.0(K/s)
 
 
   heat_capacity:
@@ -173,27 +181,27 @@ properties:
     file: ./304L_Erstarrungsdaten_edited.xlsx
     temp_col: T (K)
     prop_col: Latent heat (J/Kg)
-  #OR
-  # latent_heat_of_fusion:
+    #OR
+    # latent_heat_of_fusion:
     # key: [solidus_temperature, liquidus_temperature]
     # val: [171401, 0]
 
 
-  heat_conductivity:
-    val: [25, 30, 33, 35]  # corresponding heat_conductivity values
-    key: [25, 30, 33, 35]
+    # heat_conductivity:
+    # key: [1200, 1800, 2200, 2400]  # temperature values
+    # val: [25, 30, 33, 35]  # corresponding heat_conductivity values
 
 
-  # heat_capacity:
+    # heat_capacity:
     # key: [solidus_temperature, liquidus_temperature]
     # val: [600, 800]
-  #OR
-  # heat_capacity:
+    #OR
+    # heat_capacity:
     # key: (1000, 200)  # generates equidistant values starting 1000 with an increment of 200 until the length of val
     # [1000, 1200, 1400, 1600, 1800, 2000, 2200, 2400, 2600, 2800, 3000, 3200]
     # val: [580, 590, 600, 600, 600, 610, 620, 630, 660, 700, 750, 750]
-  #OR
-  # heat_capacity:
+    #OR
+    # heat_capacity:
     # key: [1000, 1200, 1400, 1600, ...]
     # val: [580, 590, 600, 600, 600, 610, 620, 630, 660, 700, 750, 780, 789, 799, 800, 800, 800, ...]
   #OR
diff --git a/src/pymatlib/data/alloys/SS304L/SS304L_comprehensive.yaml b/src/pymatlib/data/alloys/SS304L/SS304L_comprehensive.yaml
index 7d2278a33c7158e6e4304c4035571d21481ecc5d..9e579c59cda81d0951b5d6b320cdbe418c8c5f8b 100644
--- a/src/pymatlib/data/alloys/SS304L/SS304L_comprehensive.yaml
+++ b/src/pymatlib/data/alloys/SS304L/SS304L_comprehensive.yaml
@@ -10,12 +10,16 @@ composition:
 solidus_temperature: 1605.
 liquidus_temperature: 1735.
 
+# temperature_range: [300, 3000, 541]  # 541 (int) -> points
+#OR
+temperature_range: [300, 3000, 5.0]  # 5.0 (float) -> increment
+
 properties:
   # 1. Constant float property
   density: compute  # kg/m³
 
   # 2.1 File-based property (simple format)
-  heat_conductivity: ./heat_conductivity_temperature.txt
+  # heat_conductivity: ./heat_conductivity_temperature.txt
 
   # 2.2 File-based property (advanced format)
   specific_enthalpy:
@@ -23,6 +27,11 @@ properties:
     temp_col: T (K)
     prop_col: Enthalpy (J/kg)
 
+  heat_conductivity:
+    file: ./304L_Erstarrungsdaten_edited.xlsx
+    temp_col: T (K)
+    prop_col: Thermal conductivity (W/(m*K))-TOTAL-10000.0(K/s)
+
   # 3.1 Key-val pair with explicit temperature list
   heat_capacity:
     key: [1000, 1200, 1400, 1600, 1800, 2000, 2200, 2400, 2600, 2800]
@@ -47,14 +56,7 @@ properties:
   thermal_diffusivity: compute  # k/(rho*cp)
 
   # 4.2 Computed property (advanced format with model selection)
-  energy_density:
-    compute: total_enthalpy
-
-  # Required for energy_density computation (with number of points)
-  # energy_density_temperature_array: (300, 3000, 541)
-
-  # Alternative specification (with temperature increment)
-  energy_density_temperature_array: (300, 3000, 5.0)
+  energy_density: compute
 
   # Properties required for density computation by thermal expansion
   base_temperature: 2273.
diff --git a/tests/cpp/TestArrayContainer.py b/tests/cpp/TestArrayContainer.py
index fb9881d5dd891392be19366da0f62b0b79831d18..73bb6107291fc2f4057fbdb96ddee1ad7b74ac3d 100644
--- a/tests/cpp/TestArrayContainer.py
+++ b/tests/cpp/TestArrayContainer.py
@@ -4,6 +4,7 @@ from importlib.resources import files
 from pystencilssfg import SourceFileGenerator
 from pymatlib.core.yaml_parser import create_alloy_from_yaml
 from pymatlib.core.codegen.interpolation_array_container import InterpolationArrayContainer
+from pymatlib.core.property_array_extractor import PropertyArrayExtractor
 
 
 with SourceFileGenerator() as sfg:
@@ -22,6 +23,10 @@ with SourceFileGenerator() as sfg:
     sfg.generate(custom_container)
 
     yaml_path = files('pymatlib.data.alloys.SS304L').joinpath('SS304L.yaml')
-    mat = create_alloy_from_yaml(yaml_path, u.center())
-    arr_container = InterpolationArrayContainer.from_material("SS304L", mat)
+    mat, temp_array = create_alloy_from_yaml(yaml_path, u.center())
+    array_extractor = PropertyArrayExtractor(mat, temp_array, u.center)
+    arr_container = InterpolationArrayContainer("SS304L", temp_array, array_extractor.energy_density_array)
     sfg.generate(arr_container)
+    print(f"extractor.temperature_array: {array_extractor.temperature_array}")
+    print(f"extractor.specific_enthalpy_array: {array_extractor.specific_enthalpy_array}")
+    print(f"extractor.energy_density_array: {array_extractor.energy_density_array}")
diff --git a/tests/cpp/test_interpolation.cpp b/tests/cpp/test_interpolation.cpp
index 7b3e10ab8a3f45b5a02677242f4def5231a24748..714636097f6a69db6dfd5744825cb6fd7c461000 100644
--- a/tests/cpp/test_interpolation.cpp
+++ b/tests/cpp/test_interpolation.cpp
@@ -201,9 +201,9 @@ void test_performance() {
     // Generate random values
     std::random_device rd;
     std::mt19937 gen(rd());
-    constexpr double E_min = SS304L::E_neq.front() * 0.8;
-    constexpr double E_max = SS304L::E_neq.back() * 1.2;
-    std::uniform_real_distribution<double> dist(E_min, E_max);
+    constexpr double y_min = SS304L::y_neq.front() * 0.8;
+    constexpr double y_max = SS304L::y_neq.back() * 1.2;
+    std::uniform_real_distribution<double> dist(y_min, y_max);
 
     // Fill random energies
     for(auto& E : random_energies) {
diff --git a/tests/python/test_alloy.py b/tests/python/test_alloy.py
index cfbb44cfe278a952fa659fe4d5ee29b765025e56..354fa3d12e4f3f7e6dd8b0e159034e31a3be3211 100644
--- a/tests/python/test_alloy.py
+++ b/tests/python/test_alloy.py
@@ -3,7 +3,6 @@ import numpy as np
 import sympy as sp
 from pymatlib.core.alloy import Alloy, AlloyCompositionError, AlloyTemperatureError
 from pymatlib.data.element_data import Ti, Al, V, Fe, Cr, Mn, Ni
-from src.pymatlib.data.alloys.SS304L.SS304L import create_SS304L
 from src.pymatlib.core.typedefs import MaterialProperty
 
 def test_alloy_creation():
@@ -21,68 +20,6 @@ def test_alloy_creation():
     with pytest.raises(ValueError):
         Alloy(elements=[Fe, Cr, Mn, Ni], composition=[0.7, 0.2, 0.05, 0.05], temperature_solidus=1900., temperature_liquidus=1800.)
 
-def test_create_SS316L():
-    """Test the creation and properties of SS304L alloy."""
-    # Test with float temperature input
-    alloy = create_SS304L(1400.0)
-
-    # Check if properties are set and have correct types
-    assert hasattr(alloy, 'density')
-    assert hasattr(alloy, 'heat_capacity')
-    assert hasattr(alloy, 'heat_conductivity')
-    assert hasattr(alloy, 'thermal_diffusivity')
-
-    # Check if properties have correct values and types
-    assert alloy.density is not None
-    assert alloy.heat_capacity is not None
-    assert alloy.heat_conductivity is not None
-    assert alloy.thermal_diffusivity is not None
-
-    # Test with symbolic temperature input
-    T = sp.Symbol('T')
-    alloy_symbolic = create_SS304L(T)
-
-    # Check symbolic properties
-    assert alloy_symbolic.density is not None
-    assert alloy_symbolic.heat_capacity is not None
-    assert alloy_symbolic.heat_conductivity is not None
-    assert alloy_symbolic.thermal_diffusivity is not None
-
-    # Test property values
-    assert isinstance(float(alloy.density.expr), float)
-    assert isinstance(float(alloy.heat_capacity.expr), float)
-    assert isinstance(float(alloy.heat_conductivity.expr), float)
-    assert isinstance(float(alloy.thermal_diffusivity.expr), float)
-
-def test_create_SS316L2():
-    alloy = create_SS304L(1400.0)
-
-    # Check if density exists and has the correct type
-    assert hasattr(alloy, 'density')
-    assert alloy.density is not None
-    assert type(alloy.density).__name__ == "MaterialProperty", f"Expected MaterialProperty, got {type(alloy.density)}"
-
-def test_alloy_single_element():
-    # Test creating an alloy with a single element
-    alloy = Alloy(elements=[Fe], composition=[1.0], temperature_solidus=1800, temperature_liquidus=1900)
-    assert len(alloy.elements) == 1
-    assert alloy.composition == [1.0]
-
-def test_alloy_property_modification():
-    # Test accessing and modifying individual alloy properties
-    alloy = create_SS304L(1400.0)
-    # assert isinstance(alloy.density, MaterialProperty)
-    # Set the density to a MaterialProperty instance with a constant expression
-    alloy.density = MaterialProperty(expr=sp.Float(8000.0))
-    assert alloy.density.expr == sp.Float(8000.0)
-
-'''def test_create_SS316L_invalid_temperature():
-    # Test creating SS304L alloy with invalid temperature inputs
-    with pytest.raises(ValueError):
-        create_SS316L(-100.0)  # Negative temperature
-    with pytest.raises(ValueError):
-        create_SS316L(3500.0)  # Temperature outside valid range'''
-
 def test_valid_alloy():
     Ti64 = Alloy([Ti, Al, V], [0.90, 0.06, 0.04], 1878, 1928)
     assert np.isclose(Ti64.atomic_number, 0.90 * Ti.atomic_number + 0.06 * Al.atomic_number + 0.04 * V.atomic_number)
diff --git a/tests/python/test_yaml_config.py b/tests/python/test_yaml_config.py
index f121e0c10456047b460f66c5e8de679d910082fb..edcab721e244fdddeb740a3d7d5fc461448848ee 100644
--- a/tests/python/test_yaml_config.py
+++ b/tests/python/test_yaml_config.py
@@ -17,7 +17,7 @@ current_file = Path(__file__)
 yaml_path = current_file.parent.parent.parent / "src" / "pymatlib" / "data" / "alloys" / "SS304L" / "SS304L.yaml"
 
 # Create alloy from YAML
-ss316l = create_alloy_from_yaml(yaml_path, T)
+ss316l, temp = create_alloy_from_yaml(yaml_path, T)
 
 # Test various properties
 print(f"Elements: {ss316l.elements}")
@@ -47,4 +47,3 @@ print(f"Latent heat: {ss316l.latent_heat_of_fusion.evalf(T, test_temp)}")
 # Test array generation for energy density
 if hasattr(ss316l, 'energy_density_array'):
     print(f"\nEnergy Density Array Shape: {ss316l.energy_density_array.shape}")
-    print(f"Energy Density Temperature Array Shape: {ss316l.energy_density_temperature_array.shape}")