diff --git a/src/pymatlib/core/cpp/test_interpolation.cpp b/src/pymatlib/core/cpp/test_interpolation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d4cab14a476403a98ec8ce23fb3d52b15a6d5897
--- /dev/null
+++ b/src/pymatlib/core/cpp/test_interpolation.cpp
@@ -0,0 +1,57 @@
+#include "temperature_from_energy_density_array_cpp.h"
+#include <vector>
+#include <array>
+#include <iostream>
+#include <cassert>
+
+void run_cpp_tests() {
+    // Test with std::vector
+    std::vector<double> temperatures = {3273.15, 3263.15, 3253.15, 3243.15};
+    std::vector<double> energy_densities = {1.71e10, 1.70e10, 1.69e10, 1.68e10};
+
+
+    try {
+        // Test normal interpolation
+        double h_in = 1.695e10;
+        double result = temperature_from_energy_density_array_cpp(
+            temperatures, h_in, energy_densities);
+        std::cout << "Vector test temperature: " << result << std::endl;
+        assert(result > 3243.15 && result < 3273.15);
+
+        // Test with std::array
+        std::array<double, 4> temp_arr = {3273.15, 3263.15, 3253.15, 3243.15};
+        std::array<double, 4> energy_arr = {1.71e10, 1.70e10, 1.69e10, 1.68e10};
+        result = temperature_from_energy_density_array_cpp(temp_arr, h_in, energy_arr);
+        std::cout << "Array test temperature: " << result << std::endl;
+        assert(result > 3243.15 && result < 3273.15);
+
+        // Test boundary values
+        result = temperature_from_energy_density_array_cpp(
+            temperatures, 1.67e10, energy_densities);
+        assert(std::abs(result - 3243.15) < 1e-6);
+
+        result = temperature_from_energy_density_array_cpp(
+            temperatures, 1.72e10, energy_densities);
+        assert(std::abs(result - 3273.15) < 1e-6);
+
+        // Test error cases
+        std::vector<double> wrong_size = {1.0, 2.0};
+        try {
+            temperature_from_energy_density_array_cpp(temperatures, h_in, wrong_size);
+            std::cerr << "Expected size mismatch error" << std::endl;
+            assert(false);
+        } catch (const std::runtime_error&) {
+            std::cout << "Size mismatch test passed" << std::endl;
+        }
+
+    } catch (const std::exception& e) {
+        std::cerr << "Test error: " << e.what() << std::endl;
+        assert(false);
+    }
+}
+
+int main() {
+    run_cpp_tests();
+    std::cout << "All C++ tests completed successfully" << std::endl;
+    return 0;
+}
diff --git a/src/pymatlib/core/cpp/test_interpolation.py b/src/pymatlib/core/cpp/test_interpolation.py
new file mode 100644
index 0000000000000000000000000000000000000000..69ba4bc24d24abedd96ae78fb556a9c401cc0d86
--- /dev/null
+++ b/src/pymatlib/core/cpp/test_interpolation.py
@@ -0,0 +1,30 @@
+import numpy as np
+import pytest
+from pymatlib.core.cpp.fast_interpolation import temperature_from_energy_density_array
+
+def test_numpy_implementation():
+    # Test data
+    temperatures = np.array([3273.15, 3263.15, 3253.15, 3243.15])
+    energy_densities = np.array([1.71e10, 1.70e10, 1.69e10, 1.68e10])
+    h_in = 1.695e10
+
+    # Test normal interpolation
+    result = temperature_from_energy_density_array(temperatures, h_in, energy_densities)
+    assert 3243.15 < result < 3273.15
+    print(f"NumPy test temperature: {result}")
+
+    # Test boundary values
+    result_min = temperature_from_energy_density_array(temperatures, 1.67e10, energy_densities)
+    assert abs(result_min - 3243.15) < 1e-6
+
+    result_max = temperature_from_energy_density_array(temperatures, 1.72e10, energy_densities)
+    assert abs(result_max - 3273.15) < 1e-6
+
+    # Test error cases
+    with pytest.raises(RuntimeError):
+        wrong_size = np.array([1.0, 2.0])
+        temperature_from_energy_density_array(temperatures, h_in, wrong_size)
+
+if __name__ == "__main__":
+    test_numpy_implementation()
+    print("All Python tests completed successfully")