diff --git a/src/pymatlib/core/cpp/CMakeLists.txt b/src/pymatlib/core/cpp/CMakeLists.txt
index 26261a13834dbf58386a6d01f5ef6ec414494048..30dbbf063ab63c0e65691060a8c01345a329f481 100644
--- a/src/pymatlib/core/cpp/CMakeLists.txt
+++ b/src/pymatlib/core/cpp/CMakeLists.txt
@@ -33,19 +33,3 @@ target_link_libraries(fast_interpolation PRIVATE
         pybind11::module
         Python::NumPy
 )
-
-# Add C++ test executable
-add_executable(fast_interpolation_test_cpp test_interpolation.cpp)
-target_compile_features(fast_interpolation_test_cpp PRIVATE cxx_std_11)
-target_compile_options(fast_interpolation_test_cpp PRIVATE -O3)
-target_link_libraries(fast_interpolation_test_cpp PRIVATE
-        interpolation_template
-)
-
-# Add custom target for running Python tests
-add_custom_target(fast_interpolation_test_py
-        COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_interpolation.py
-        DEPENDS fast_interpolation
-        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-        COMMENT "Running Python tests"
-)
diff --git a/src/pymatlib/data/alloys/SS316L/SS316L.py b/src/pymatlib/data/alloys/SS316L/SS316L.py
index cb4014911f25fa21b4e57fef3e57187cd99bb0ff..b346e0b3024cce7788b1fd9247526d99a4784fcd 100644
--- a/src/pymatlib/data/alloys/SS316L/SS316L.py
+++ b/src/pymatlib/data/alloys/SS316L/SS316L.py
@@ -6,7 +6,7 @@ from typing import Union
 from matplotlib import pyplot as plt
 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
+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
 from pymatlib.core.cpp.fast_interpolation import interpolate_binary_search, interpolate_double_lookup
@@ -76,9 +76,9 @@ def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy:
     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 = read_data_from_excel(density_data_file_path, temp_col="T (K)", prop_col="Enthalpy (J/g)-TOTAL-10000.0(K/s)")
+    _, 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 = read_data_from_excel(density_data_file_path, temp_col="T (K)", prop_col="u1 (rho*h)")
+    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)")
 
@@ -93,12 +93,14 @@ def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy:
     SS316L.thermal_diffusivity = thermal_diffusivity_by_heat_conductivity(SS316L.heat_conductivity, SS316L.density, SS316L.heat_capacity)
     # SS316L.latent_heat_of_fusion = interpolate_property(T, SS316L.solidification_interval(), np.array([171401.0, 0.0]))
     SS316L.latent_heat_of_fusion = interpolate_property(T, latent_heat_of_fusion_temp_array, latent_heat_of_fusion_array)
-    SS316L.energy_density = energy_density(T, SS316L.density, SS316L.heat_capacity, SS316L.latent_heat_of_fusion)
+    SS316L.specific_enthalpy = interpolate_property(T, density_temp_array, sp_enthalpy_array)
+    SS316L.energy_density = energy_density_total_enthalpy(SS316L.density, SS316L.specific_enthalpy)
     SS316L.energy_density_solidus = SS316L.energy_density.evalf(T, SS316L.temperature_solidus)
     SS316L.energy_density_liquidus = SS316L.energy_density.evalf(T, SS316L.temperature_liquidus)
 
     # Populate temperature_array and energy_density_array
     SS316L.temperature_array = density_temp_array
+    SS316L.energy_density_temperature_array = u1_temp_array
 
     SS316L.energy_density_array = np.array([
         SS316L.energy_density.evalf(T, temp) for temp in density_temp_array
@@ -109,17 +111,36 @@ def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy:
             SS316L.energy_density_array)
 
     E = SS316L.energy_density.evalf(T, SS316L.temperature_liquidus)
-    T_eq, E_neq, E_eq, inv_delta_E_eq, idx_map = prepare_interpolation_arrays(
+    result = prepare_interpolation_arrays(
         SS316L.temperature_array,
         SS316L.energy_density_array
     )
-    args3 = (E, T_eq, E_neq, E_eq, inv_delta_E_eq, idx_map)
 
-    start_time3 = time.perf_counter()
-    T_interpolate3 = interpolate_double_lookup(*args3)
-    execution_time3 = time.perf_counter() - start_time3
-    print(f"Interpolated temperature: {T_interpolate3}")
-    print(f"Execution time: {execution_time3:.8f} seconds")
+    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}")
diff --git a/src/pymatlib/core/constants.py b/src/pymatlib/data/constants.py
similarity index 100%
rename from src/pymatlib/core/constants.py
rename to src/pymatlib/data/constants.py
diff --git a/src/pymatlib/data/element_data.py b/src/pymatlib/data/element_data.py
index b52b1433f053b33e150d5e34cc1857ff3aea751d..4f85e2ac183942bfb452a752fe0c4392539608e0 100644
--- a/src/pymatlib/data/element_data.py
+++ b/src/pymatlib/data/element_data.py
@@ -1,4 +1,4 @@
-from pymatlib.core.constants import Constants
+from pymatlib.data.constants import Constants
 from pymatlib.core.elements import ChemicalElement
 
 # NIST: National Institute of Standards and Technology
diff --git a/tests/cpp/CMakeLists.txt b/tests/cpp/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0e359b80bec46fbbb7afc73004169a1de6b6e3e8
--- /dev/null
+++ b/tests/cpp/CMakeLists.txt
@@ -0,0 +1,33 @@
+cmake_minimum_required(VERSION 3.10)
+project(pymatlib_tests)
+
+# Find required packages
+find_package(Python COMPONENTS Interpreter Development NumPy REQUIRED)
+execute_process(
+        COMMAND ${Python_EXECUTABLE} -m pybind11 --cmakedir
+        OUTPUT_VARIABLE PYBIND11_CMAKE_DIR
+        OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+find_package(pybind11 PATHS ${PYBIND11_CMAKE_DIR} NO_DEFAULT_PATH REQUIRED)
+
+# Get the path to the source directory
+set(SRC_DIR "${CMAKE_CURRENT_SOURCE_DIR}/../../src/pymatlib/core/cpp")
+
+# Create interface library for header-only template implementation
+add_library(interpolation_template INTERFACE)
+target_include_directories(interpolation_template INTERFACE
+        ${SRC_DIR}
+        ${SRC_DIR}/include
+)
+
+# Add C++ test executable
+add_executable(fast_interpolation_test_cpp test_interpolation.cpp)
+target_include_directories(fast_interpolation_test_cpp PRIVATE
+        ${SRC_DIR}
+        ${SRC_DIR}/include
+)
+target_compile_features(fast_interpolation_test_cpp PRIVATE cxx_std_11)
+target_compile_options(fast_interpolation_test_cpp PRIVATE -O3)
+target_link_libraries(fast_interpolation_test_cpp PRIVATE
+        interpolation_template
+)
diff --git a/src/pymatlib/core/cpp/TestArrayContainer.py b/tests/cpp/TestArrayContainer.py
similarity index 100%
rename from src/pymatlib/core/cpp/TestArrayContainer.py
rename to tests/cpp/TestArrayContainer.py
diff --git a/src/pymatlib/core/cpp/test_interpolation.cpp b/tests/cpp/test_interpolation.cpp
similarity index 100%
rename from src/pymatlib/core/cpp/test_interpolation.cpp
rename to tests/cpp/test_interpolation.cpp
diff --git a/tests/legacy/__init__.py b/tests/legacy/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/pymatlib/core/test_interpolate_functs_perf.py b/tests/legacy/test_interpolate_functs_perf.py
similarity index 95%
rename from src/pymatlib/core/test_interpolate_functs_perf.py
rename to tests/legacy/test_interpolate_functs_perf.py
index dc575586f8d67d55d4c5cf1b62c19e9aa852e4a7..4cb35850857fa8e3785bfb5229bf32687dafefbe 100644
--- a/src/pymatlib/core/test_interpolate_functs_perf.py
+++ b/tests/legacy/test_interpolate_functs_perf.py
@@ -1,3 +1,20 @@
+"""
+DEPRECATED: This file contains the original implementation of interpolation comparision.
+It is kept for historical reference only and should not be used in new code.
+
+As of March 2025, please use the new dictionary-based implementation in
+pymatlib.core.interpolators instead.
+"""
+
+import warnings
+
+warnings.warn(
+    "This module is deprecated and will be removed in a future version. "
+    "Please use pymatlib.core.interpolators instead.",
+    DeprecationWarning,
+    stacklevel=2
+)
+
 import time
 import numpy as np
 import sympy as sp
diff --git a/tests/python/__init__.py b/tests/python/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/tests/test_alloy.py b/tests/python/test_alloy.py
similarity index 100%
rename from tests/test_alloy.py
rename to tests/python/test_alloy.py
diff --git a/tests/test_assignment_converter.py b/tests/python/test_assignment_converter.py
similarity index 100%
rename from tests/test_assignment_converter.py
rename to tests/python/test_assignment_converter.py
diff --git a/tests/test_data_handler.py b/tests/python/test_data_handler.py
similarity index 100%
rename from tests/test_data_handler.py
rename to tests/python/test_data_handler.py
diff --git a/src/pymatlib/core/cpp/test_interpolation.py b/tests/python/test_interpolation.py
similarity index 100%
rename from src/pymatlib/core/cpp/test_interpolation.py
rename to tests/python/test_interpolation.py
diff --git a/tests/test_interpolators.py b/tests/python/test_interpolators.py
similarity index 100%
rename from tests/test_interpolators.py
rename to tests/python/test_interpolators.py
diff --git a/tests/test_models.py b/tests/python/test_models.py
similarity index 100%
rename from tests/test_models.py
rename to tests/python/test_models.py
diff --git a/tests/test_typedefs.py b/tests/python/test_typedefs.py
similarity index 100%
rename from tests/test_typedefs.py
rename to tests/python/test_typedefs.py
diff --git a/src/pymatlib/data/alloys/SS316L/test_yaml_config.py b/tests/python/test_yaml_config.py
similarity index 87%
rename from src/pymatlib/data/alloys/SS316L/test_yaml_config.py
rename to tests/python/test_yaml_config.py
index 6eb31369f22497c15d4b95fa6da39dcd9d47261a..d3709938f4a8979fe083d08467fae474c1986405 100644
--- a/src/pymatlib/data/alloys/SS316L/test_yaml_config.py
+++ b/tests/python/test_yaml_config.py
@@ -12,8 +12,13 @@ def print_property_value(prop, T, temp):
 # Create symbolic temperature variable
 T = sp.Symbol('T')
 
+# Get the path to the YAML file
+current_file = Path(__file__)
+yaml_path = current_file.parent.parent.parent / "src" / "pymatlib" / "data" / "alloys" / "SS316L" / "SS304L.yaml"
+
 # Create alloy from YAML
-ss316l = create_alloy_from_yaml("SS304L.yaml", T)
+ss316l = create_alloy_from_yaml(yaml_path, T)
+#ss316l_1 = create_alloy_from_yaml("SS304L_1.yaml", T)
 
 # Test various properties
 print(f"Elements: {ss316l.elements}")