diff --git a/src/pymatlib/core/cpp/test_interpolation.cpp b/src/pymatlib/core/cpp/test_interpolation.cpp
index fb9ef34f4cb7a58513872c09f2ed716c3e7c31b7..ec001f0309bed1176e31094bf3676f9ed0829443 100644
--- a/src/pymatlib/core/cpp/test_interpolation.cpp
+++ b/src/pymatlib/core/cpp/test_interpolation.cpp
@@ -1,604 +1,288 @@
 #include "interpolate_binary_search_cpp.h"
 #include "interpolate_double_lookup_cpp.h"
 #include <array>
-#include <vector>
 #include <random>
 #include <chrono>
 #include <iomanip>
 #include <cassert>
 #include <iostream>
+#include <numeric>
 #include "TestArrayContainer.hpp"
 
-/*
-void run_comprehensive_tests1() {
-    std::cout << "\nRunning comprehensive interpolation tests...\n" << std::endl;
-
-    // Test Case 1: Basic ascending arrays
-    {
-        std::vector<double> T = {100.0, 200.0, 300.0, 400.0, 500.0};
-        std::vector<double> E = {1.0, 2.0, 3.0, 4.0, 5.0};
-
-        // Test middle interpolation
-        double result = interpolate_binary_search_cpp(T, 2.5, E);
-        assert(std::abs(result - 250.0) < 1e-8);
-        std::cout << "Basic ascending test passed" << std::endl;
-    }
-
-    // Test Case 2: Descending arrays
-    {
-        std::vector<double> T = {500.0, 400.0, 300.0, 200.0, 100.0};
-        std::vector<double> E = {5.0, 4.0, 3.0, 2.0, 1.0};
-
-        double result = interpolate_binary_search_cpp(T, 2.5, E);
-        assert(std::abs(result - 250.0) < 1e-8);
-        std::cout << "Descending array test passed" << std::endl;
-    }
-
-    // Test Case 3: Edge cases
-    {
-        std::vector<double> T = {100.0, 200.0, 300.0};
-        std::vector<double> E = {1.0, 2.0, 3.0};
-
-        // Test below minimum
-        double result = interpolate_binary_search_cpp(T, 0.5, E);
-        assert(std::abs(result - 100.0) < 1e-8);
-
-        // Test above maximum
-        result = interpolate_binary_search_cpp(T, 3.5, E);
-        assert(std::abs(result - 300.0) < 1e-8);
-
-        std::cout << "Edge cases test passed" << std::endl;
-    }
-
-    // Test Case 4: Double lookup specific tests
-    {
-        std::vector<double> T_eq = {100.0, 200.0, 300.0, 400.0, 500.0};
-        std::vector<double> E_neq = {1.0, 2.1, 3.4, 4.2, 5.0};
-        std::vector<double> E_eq = {1.0, 2.0, 3.0, 4.0, 5.0};
-        double inv_delta_E_eq = 1.0;
-        std::vector<double> idx_map = {0, 0, 1, 2, 3};
-
-        double result = interpolate_double_lookup_cpp(2.5, T_eq, E_neq, E_eq,
-                                                    inv_delta_E_eq, idx_map);
-        std::cout << result << std::endl;
-        // assert(std::abs(result - 230.769) < 1e-6);
-        std::cout << "Double lookup basic test passed" << std::endl;
-    }
-
-    // Test Case 5: Performance comparison
-    {
-        const size_t size = 1000000;
-        std::vector<double> T(size);
-        std::vector<double> E(size);
-        std::vector<double> E_eq(size);
-        std::vector<double> idx_map(size);
-
-        // Initialize arrays
-        for(size_t i = 0; i < size; ++i) {
-            T[i] = static_cast<double>(i);
-            E[i] = static_cast<double>(i) * 1.5;
-            E_eq[i] = static_cast<double>(i);
-            idx_map[i] = i;
-        }
-
-        double inv_delta_E_eq = 1.0;
-
-        // Time binary search
-        auto start = std::chrono::high_resolution_clock::now();
-        double result1 = interpolate_binary_search_cpp(T, 500000.5, E);
-        std::cout << result1 << std::endl;
-        auto end = std::chrono::high_resolution_clock::now();
-        auto binary_duration = std::chrono::duration_cast<std::chrono::nanoseconds>
-                             (end - start).count();
-
-        // Time double lookup
-        start = std::chrono::high_resolution_clock::now();
-        double result2 = interpolate_double_lookup_cpp(500000.5, T, E, E_eq,
-                                                     inv_delta_E_eq, idx_map);
-        std::cout << result2 << std::endl;
-        end = std::chrono::high_resolution_clock::now();
-        auto double_duration = std::chrono::duration_cast<std::chrono::nanoseconds>
-                             (end - start).count();
-
-        std::cout << "\nPerformance test results:" << std::endl;
-        std::cout << "Binary search: " << binary_duration << " ns" << std::endl;
-        std::cout << "Double lookup: " << double_duration << " ns" << std::endl;
-    }
-
-    // Test Case 6: Error handling
-    {
-        try {
-            std::vector<double> T = {100.0};
-            std::vector<double> E = {1.0};
-            interpolate_binary_search_cpp(T, 1.5, E);
-            assert(false && "Should throw exception for size < 2");
-        } catch (const std::runtime_error&) {
-            std::cout << "Size validation test passed" << std::endl;
-        }
-    }
-
-    // Test Case 7: Non-uniform spacing
-    {
-        std::vector<double> T = {100.0, 150.0, 300.0, 320.0, 500.0};
-        std::vector<double> E = {1.0, 1.5, 3.0, 3.2, 5.0};
-
-        double result = interpolate_binary_search_cpp(T, 2.0, E);
-        std::cout << "Non-uniform spacing test passed" << std::endl;
-    }
-}
-*/
-
-/*
-void run_comprehensive_tests2() {
-    std::cout << "\nRunning comprehensive interpolation tests...\n" << std::endl;
-
-    // Test Case 1: Basic tests with both std::vector and std::array
-    {
-        // Vector version
-        std::vector<double> T_vec = {100.0, 200.0, 300.0, 400.0};
-        std::vector<double> E_vec = {1.0, 2.0, 3.0, 4.0};
-        double result_vec = interpolate_binary_search_cpp(T_vec, 2.5, E_vec);
-        assert(std::abs(result_vec - 250.0) < 1e-6);
-
-        // Array version
-        std::array<double, 4> T_arr = {100.0, 200.0, 300.0, 400.0};
-        std::array<double, 4> E_arr = {1.0, 2.0, 3.0, 4.0};
-        double result_arr = interpolate_binary_search_cpp(T_arr, 2.5, E_arr);
-        assert(std::abs(result_arr - 250.0) < 1e-6);
-
-        std::cout << "Basic vector/array tests passed" << std::endl;
-    }
-
-    // Test Case 2: Descending order with both containers
-    {
-        // Vector version
-        std::vector<double> T_vec = {400.0, 300.0, 200.0, 100.0};
-        std::vector<double> E_vec = {4.0, 3.0, 2.0, 1.0};
-        double result_vec = interpolate_binary_search_cpp(T_vec, 2.5, E_vec);
-
-        // Array version
-        std::array<double, 4> T_arr = {400.0, 300.0, 200.0, 100.0};
-        std::array<double, 4> E_arr = {4.0, 3.0, 2.0, 1.0};
-        double result_arr = interpolate_binary_search_cpp(T_arr, 2.5, E_arr);
-
-        assert(std::abs(result_vec - result_arr) < 1e-8);
-        std::cout << "Descending order vector/array tests passed" << std::endl;
-    }
-
-    // Test Case 3: Edge cases for both containers
-    {
-        // Vector version
-        std::vector<double> T_vec = {3273.15, 3263.15, 3253.15, 3243.15};
-        std::vector<double> E_vec = {1.71e10, 1.70e10, 1.69e10, 1.68e10};
-
-        double result_vec_min = interpolate_binary_search_cpp(T_vec, 1.67e10, E_vec);
-        double result_vec_max = interpolate_binary_search_cpp(T_vec, 1.72e10, E_vec);
-
-        // Array version
-        std::array<double, 4> T_arr = {3273.15, 3263.15, 3253.15, 3243.15};
-        std::array<double, 4> E_arr = {1.71e10, 1.70e10, 1.69e10, 1.68e10};
-
-        double result_arr_min = interpolate_binary_search_cpp(T_arr, 1.67e10, E_arr);
-        double result_arr_max = interpolate_binary_search_cpp(T_arr, 1.72e10, E_arr);
-
-        assert(std::abs(result_vec_min - 3243.15) < 1e-6);
-        assert(std::abs(result_arr_min - 3243.15) < 1e-6);
-        assert(std::abs(result_vec_max - 3273.15) < 1e-6);
-        assert(std::abs(result_arr_max - 3273.15) < 1e-6);
-
-        std::cout << "Edge cases vector/array tests passed" << std::endl;
-    }
-
-    // Test Case 4: Double lookup with both containers
-    {
-        // Vector version
-        std::vector<double> T_vec = {100.0, 200.0, 300.0, 400.0};
-        std::vector<double> E_neq_vec = {1.0, 2.1, 3.4, 4.2};
-        std::vector<double> E_eq_vec = {1.0, 2.0, 3.0, 4.0};
-        std::vector<double> idx_map_vec = {0, 0, 1, 2};
-
-        // Array version
-        std::array<double, 4> T_arr = {100.0, 200.0, 300.0, 400.0};
-        std::array<double, 4> E_neq_arr = {1.0, 2.1, 3.4, 4.2};
-        std::array<double, 4> E_eq_arr = {1.0, 2.0, 3.0, 4.0};
-        std::array<double, 4> idx_map_arr = {0, 0, 1, 2};
-
-        double inv_delta_E_eq = 1.0;
-
-        double result_vec = interpolate_double_lookup_cpp(2.5, T_vec, E_neq_vec,
-                                                        E_eq_vec, inv_delta_E_eq,
-                                                        idx_map_vec);
-
-        double result_arr = interpolate_double_lookup_cpp(2.5, T_arr, E_neq_arr,
-                                                        E_eq_arr, inv_delta_E_eq,
-                                                        idx_map_arr);
-
-        assert(std::abs(result_vec - result_arr) < 1e-6);
-        std::cout << "Double lookup vector/array tests passed" << std::endl;
-    }
-
-    // Test Case 5: Error handling for both containers
-    {
-        try {
-            std::vector<double> T_vec = {100.0};
-            std::vector<double> E_vec = {1.0};
-            interpolate_binary_search_cpp(T_vec, 1.5, E_vec);
-            assert(false && "Should throw exception for vector size < 2");
-        } catch (const std::runtime_error&) {
-            std::cout << "Vector size validation passed" << std::endl;
-        }
-
-        try {
-            std::array<double, 1> T_arr = {100.0};
-            std::array<double, 1> E_arr = {1.0};
-            interpolate_binary_search_cpp(T_arr, 1.5, E_arr);
-            assert(false && "Should throw exception for array size < 2");
-        } catch (const std::runtime_error&) {
-            std::cout << "Array size validation passed" << std::endl;
-        }
-    }
-
-    // Test Case 6: Performance comparison between vector and array
-    {
-        const size_t size = 1000000;
-        std::vector<double> T_vec(size);
-        std::vector<double> E_vec(size);
-        std::array<double, 1000> T_arr;
-        std::array<double, 1000> E_arr;
-
-        // Initialize containers
-        for(size_t i = 0; i < size; ++i) {
-            T_vec[i] = static_cast<double>(i);
-            E_vec[i] = static_cast<double>(i) * 1.5;
-        }
-        for(size_t i = 0; i < 1000; ++i) {
-            T_arr[i] = static_cast<double>(i);
-            E_arr[i] = static_cast<double>(i) * 1.5;
-        }
-
-        // Time vector version
-        auto start = std::chrono::high_resolution_clock::now();
-        double result_vec = interpolate_binary_search_cpp(T_vec, 500000.5, E_vec);
-        auto vector_duration = std::chrono::duration_cast<std::chrono::nanoseconds>
-                             (std::chrono::high_resolution_clock::now() - start).count();
-
-        // Time array version
-        start = std::chrono::high_resolution_clock::now();
-        double result_arr = interpolate_binary_search_cpp(T_arr, 500.5, E_arr);
-        auto array_duration = std::chrono::duration_cast<std::chrono::nanoseconds>
-                            (std::chrono::high_resolution_clock::now() - start).count();
-
-        std::cout << "\nPerformance comparison:" << std::endl;
-        std::cout << "Vector duration: " << vector_duration << " ns" << std::endl;
-        std::cout << "Array duration:  " << array_duration << " ns" << std::endl;
-    }
-}
-*/
-
-
 // Helper function to compare floating point numbers
 bool is_equal(double a, double b, double tolerance = 1e-10) {
     return std::abs(a - b) < tolerance;
 }
 
-
-void test_basic_functionality1() {
-    // Use realistic temperature-enthalpy values from search results
-    const std::vector<double> T_vec = {2400.0, 2500.0, 2600.0, 2700.0};
-    const std::vector<double> E_vec = {1.68e10, 1.69e10, 1.70e10, 1.71e10};
-
-    // Test point should be within the data range
-    const double test_E = 1.695e10;
-    const double expected_T = 2550.0;  // Midpoint interpolation
-
-    // Add tolerance appropriate for temperature values
-    const double tolerance = 1e-6;
-
-    DoubleLookupTests tests;
-    double result_vec = interpolate_binary_search_cpp(test_E, tests);
-    assert(std::abs(result_vec - expected_T) < tolerance);
-}
-
-
 void test_basic_functionality() {
     std::cout << "\nTesting basic functionality..." << std::endl;
 
-    // Setup test data with precise values
-    // Ascending order vectors
-    const std::vector<double> T_vec = {3243.15, 3253.15, 3263.15, 3273.15};
-    const std::vector<double> E_vec = {1.68e10, 1.69e10, 1.70e10, 1.71e10};
-    // Ascending order arrays
-    const std::array<double, 4> T_arr = {3243.15, 3253.15, 3263.15, 3273.15};
-    const std::array<double, 4> E_arr = {1.68e10, 1.69e10, 1.70e10, 1.71e10};
-
     // Test middle point interpolation
-    const double test_E = 1.695e10;
-    const double expected_T = 3258.15;
+    const double test_E = 16950000000.0;
 
     // Binary Search Tests
     {
-        DoubleLookupTests tests;
-        double result_arr = interpolate_binary_search_cpp(test_E, tests);
-        assert(is_equal(result_arr, expected_T));
-
-        // Verify vector and array results match
-        assert(is_equal(result_vec, result_arr));
+        BinarySearchTests tests;
+        double result = tests.interpolateBS(test_E);
+        // Expected value based on T_bs array
+        double expected_T = 3253.15;
+        assert(is_equal(result, expected_T));
+        std::cout << "Basic binary search interpolation passed" << std::endl;
     }
 
     // Double Lookup Tests
     {
-        // Setup double lookup specific data
-        // std::vector<double> E_eq_vec = {1.68e10, 1.69e10, 1.70e10, 1.71e10};
-        std::array<double, 4> E_eq_arr = {1.68e10, 1.69e10, 1.70e10, 1.71e10};
-        // std::vector<double> idx_map_vec = {0, 1, 2, 3};
-        std::array<double, 4> idx_map_arr = {0, 1, 2, 3};
-        const double inv_delta_E_eq = 1.0 / (1e9);
-
-        /*double result_vec = interpolate_double_lookup_cpp(
-            test_E, T_vec, E_vec, E_eq_vec, inv_delta_E_eq, idx_map_vec);
-        assert(is_equal(result_vec, expected_T));*/
-
         DoubleLookupTests tests;
-        double result_arr = tests.interpolateDL(test_E);
-        /*double result_arr = interpolate_double_lookup_cpp(
-            test_E, T_arr, E_arr, E_eq_arr, inv_delta_E_eq, idx_map_arr);*/
-        assert(is_equal(result_arr, expected_T));
-
-        // Verify vector and array results match
-        // assert(is_equal(result_vec, result_arr));
+        double result = tests.interpolateDL(test_E);
+        // Expected value based on T_eq array
+        double expected_T = 3258.15;
+        assert(is_equal(result, expected_T));
+        std::cout << "Basic double lookup interpolation passed" << std::endl;
     }
 }
 
-
 void test_edge_cases_and_errors() {
-    std::cout << "Testing edge cases and error handling..." << std::endl;
-
-    // Test data for vectors
-    const std::vector<double> T_vec = {3243.15, 3253.15, 3263.15, 3273.15};
-    const std::vector<double> E_vec = {1.68e10, 1.69e10, 1.70e10, 1.71e10};
+    std::cout << "\nTesting edge cases and error handling..." << std::endl;
 
-    // Test data for arrays
-    const std::array<double, 4> T_arr = {3243.15, 3253.15, 3263.15, 3273.15};
-    const std::array<double, 4> E_arr = {1.68e10, 1.69e10, 1.70e10, 1.71e10};
-
-    // Test boundary values
+    // Test boundary values for binary search
     {
-        // Vector tests
-        /*double result_vec_min = interpolate_binary_search_cpp(T_vec, 1.67e10, E_vec);
-        double result_vec_max = interpolate_binary_search_cpp(T_vec, 1.72e10, E_vec);
-        assert(is_equal(result_vec_min, 3243.15));
-        assert(is_equal(result_vec_max, 3273.15));*/
+        BinarySearchTests tests;
 
-        // Array tests
-        DoubleLookupTests tests;
-        double result_arr_min = interpolate_binary_search_cpp(1.67e10, tests);
-        double result_arr_max = interpolate_binary_search_cpp(1.72e10, tests);
-        assert(is_equal(result_arr_min, 3243.15));
-        assert(is_equal(result_arr_max, 3273.15));
-    }
+        // Test below minimum
+        double result_min = tests.interpolateBS(16750000000.0);
+        assert(is_equal(result_min, 3243.15));
 
-    // Error cases
-    /*{
-        // Test vector size mismatch
-        std::vector<double> wrong_size_vec = {1.0};
-        bool caught_vector_error = false;
-        try {
-            interpolate_binary_search_cpp(T_vec, 1.70e10, wrong_size_vec);
-        } catch (const std::runtime_error&) {
-            caught_vector_error = true;
-        }
-        assert(caught_vector_error);
-
-        // For arrays, test with non-monotonic values
-        std::array<double, 4> non_monotonic_arr = {1.71e10, 1.69e10, 1.70e10, 1.68e10};  // Values not in order
-        bool caught_array_error = false;
-        try {
-            interpolate_binary_search_cpp(T_arr, 1.70e10, non_monotonic_arr);
-        } catch (const std::runtime_error&) {
-            caught_array_error = true;
-        }
-        assert(caught_array_error);
+        // Test above maximum
+        double result_max = tests.interpolateBS(17150000000.0);
+        assert(is_equal(result_max, 3278.15));
 
-        // Compile-time size check for arrays
-        static_assert(T_arr.size() >= 2, "Array size must be at least 2");
-        static_assert(E_arr.size() >= 2, "Array size must be at least 2");
-    }*/
+        std::cout << "Edge cases for binary search passed" << std::endl;
+    }
 
     // Double Lookup edge cases
     {
-        std::vector<double> E_eq_vec = {1.68e10, 1.69e10, 1.70e10, 1.71e10};
-        std::vector<double> idx_map_vec = {0, 1, 2, 3};
-        std::array<double, 4> E_eq_arr = {1.68e10, 1.69e10, 1.70e10, 1.71e10};
-        std::array<double, 4> idx_map_arr = {0, 1, 2, 3};
-        const double inv_delta_E_eq = 1.0 / (1e9);
+        DoubleLookupTests tests;
 
+        // Test below minimum
+        double result_min = tests.interpolateDL(16750000000.0);
+        assert(is_equal(result_min, 3243.15));
 
+        // Test above maximum
+        double result_max = tests.interpolateDL(17150000000.0);
+        assert(is_equal(result_max, 3273.15));
 
-        // Array tests
-        DoubleLookupTests tests;
-        double result_arr_min = tests.interpolateDL(1.67e10);
-        assert(is_equal(result_arr_min, 3243.15));
+        std::cout << "Edge cases for double lookup passed" << std::endl;
     }
 }
 
-
 void test_interpolation_accuracy() {
     std::cout << "\nTesting interpolation accuracy..." << std::endl;
 
     // Test data with precise values
     struct TestCase {
         double input_E;
-        double expected_T;
+        double expected_T_bs;
+        double expected_T_dl;
         double tolerance;
         std::string description;
     };
 
     const std::vector<TestCase> test_cases = {
-        // For E=1.705e10 (between 1.70e10 and 1.71e10)
-        {1.705e10, 3268.15, 1e-10, "Mid-point interpolation"},
-        // For E=1.695e10 (between 1.69e10 and 1.70e10)
-        {1.695e10, 3258.15, 1e-10, "Quarter-point interpolation"},
-        // For E=1.685e10 (between 1.68e10 and 1.69e10)
-        {1.685e10, 3248.15, 1e-10, "Eighth-point interpolation"},
+        // Test cases with expected interpolated values for both methods
+        {16850000000.0, 3245.65, 3248.15, 1e-10, "Quarter-point interpolation"},
+        {16950000000.0, 3253.15, 3258.15, 1e-10, "Mid-point interpolation"},
+        {17050000000.0, 3268.15, 3268.15, 1e-10, "Three-quarter-point interpolation"},
         // Near exact points
-        {1.70e10, 3263.15, 1e-8, "Near-exact point"},
-        {1.69e10, 3253.15, 1e-8, "Near-exact point"}
-    };
-
-    // Setup arrays with high precision values
-    // For std::vector
-    const std::vector<double> T_vec = {
-        3243.15000000, 3253.15000000, 3263.15000000, 3273.15000000
-    };
-    const std::vector<double> E_vec = {
-        1.68000000e10, 1.69000000e10, 1.70000000e10, 1.71000000e10
-    };
-
-    // For std::array
-    const std::array<double, 4> T_arr = {
-        3243.15000000, 3253.15000000, 3263.15000000, 3273.15000000
-    };
-    const std::array<double, 4> E_arr = {
-        1.68000000e10, 1.69000000e10, 1.70000000e10, 1.71000000e10
+        {16900000000.0, 3248.15, 3253.15, 1e-8, "Near-exact point"},
+        {17000000000.0, 3258.15, 3263.15, 1e-8, "Near-exact point"}
     };
 
-    // Setup for double lookup
-    std::vector<double> E_eq_vec = {
-        1.68000000e10, 1.69000000e10, 1.70000000e10, 1.71000000e10
-    };
-    std::array<double, 4> E_eq_arr = {
-        1.68000000e10, 1.69000000e10, 1.70000000e10, 1.71000000e10
-    };
-    std::vector<double> idx_map_vec = {0, 1, 2, 3};
-    std::array<double, 4> idx_map_arr = {0, 1, 2, 3};
-    const double inv_delta_E_eq = 1.0 / (1e9);
-
-    DoubleLookupTests tests;
     for (const auto& test : test_cases) {
         // Binary Search Tests
-        double result_bin_arr = interpolate_binary_search_cpp(test.input_E, tests);
-
-        assert(is_equal(result_bin_vec, test.expected_T, test.tolerance));
-        assert(is_equal(result_bin_arr, test.expected_T, test.tolerance));
-        assert(is_equal(result_bin_vec, result_bin_arr, test.tolerance));
+        BinarySearchTests bs_tests;
+        double result_bs = bs_tests.interpolateBS(test.input_E);
+        assert(is_equal(result_bs, test.expected_T_bs, test.tolerance));
 
         // Double Lookup Tests
-        double result_dl_arr = tests.interpolateDL(test.input_E);
+        DoubleLookupTests dl_tests;
+        double result_dl = dl_tests.interpolateDL(test.input_E);
+        assert(is_equal(result_dl, test.expected_T_dl, test.tolerance));
 
-        // assert(is_equal(result_dl_vec, test.expected_T, test.tolerance));
-        assert(is_equal(result_dl_arr, test.expected_T, test.tolerance));
-        // assert(is_equal(result_dl_vec, result_dl_arr, test.tolerance));
+        std::cout << "Accuracy test passed for " << test.description << std::endl;
     }
 }
 
+void test_consistency() {
+    std::cout << "\nTesting interpolation consistency..." << std::endl;
 
-void test_performance() {
-    std::cout << "\nTesting performance..." << std::endl;
+    // Test that small changes in input produce correspondingly small changes in output
+    BinarySearchTests bs_tests;
+    DoubleLookupTests dl_tests;
 
-    // Test with different sizes
-    const std::vector<size_t> test_sizes = {100, 1000, 10000, 100000};
+    const double base_E = 16900000000.0;
+    const double delta_E = 1000000.0;  // Small change in energy
 
-    for (size_t size : test_sizes) {
-        std::cout << "\nTesting with array size: " << size << std::endl;
+    double base_T_bs = bs_tests.interpolateBS(base_E);
+    double delta_T_bs = bs_tests.interpolateBS(base_E + delta_E) - base_T_bs;
 
-        // Generate test data
-        std::vector<double> T_vec(size);
-        std::vector<double> E_vec(size);
-        std::vector<double> E_eq_vec(size);
-        std::vector<double> idx_map_vec(size);
+    double base_T_dl = dl_tests.interpolateDL(base_E);
+    double delta_T_dl = dl_tests.interpolateDL(base_E + delta_E) - base_T_dl;
 
-        // Fill with realistic values
-        for (size_t i = 0; i < size; ++i) {
-            T_vec[i] = 3273.15 + i * 0.01;
-            E_vec[i] = 1.71e10 + i * 1e6 + 1e3 * i/2.;
-            E_eq_vec[i] = 1.71e10 + i * 1e6 + 1e3 * i/2.;
-            idx_map_vec[i] = i;
-        }
+    // Check that changes are reasonable (not too large for small input change)
+    assert(std::abs(delta_T_bs) < 1.0);
+    assert(std::abs(delta_T_dl) < 1.0);
 
-        // Create smaller arrays for comparison
-        constexpr size_t arr_size = 1000;
-        std::array<double, arr_size> T_arr;
-        std::array<double, arr_size> E_arr;
-        std::array<double, arr_size> E_eq_arr;
-        std::array<double, arr_size> idx_map_arr;
-
-        for (size_t i = 0; i < arr_size; ++i) {
-            T_arr[i] = T_vec[i];
-            E_arr[i] = E_vec[i];
-            E_eq_arr[i] = E_eq_vec[i];
-            idx_map_arr[i] = idx_map_vec[i];
-        }
+    std::cout << "Consistency test passed" << std::endl;
+}
 
-        const double inv_delta_E_eq = 1.0 / (1e6);
-        const int num_tests = 1000;
+void test_stress() {
+    std::cout << "\nPerforming stress testing..." << std::endl;
 
-        // Test points
-        /*std::vector<double> test_points(num_tests);
-        for (int i = 0; i < num_tests; ++i) {
-            test_points[i] = 1.70e10 + (i * 1e6);
-        }*/
+    BinarySearchTests bs_tests;
+    DoubleLookupTests dl_tests;
 
-        // Generate random monotonically increasing test points
-        std::vector<double> test_points(num_tests);
-        test_points[0] = 1.70e10;  // Start value
-        std::random_device rd;
-        std::mt19937 gen(rd());
-        std::uniform_real_distribution<double> dist(1e5, 1e6);  // Random increments between 1e5 and 1e6
+    // Test with extremely large values
+    double large_E = 1.0e20;
+    double result_large_bs = bs_tests.interpolateBS(large_E);
+    double result_large_dl = dl_tests.interpolateDL(large_E);
 
-        for (int i = 1; i < num_tests; ++i) {
-            test_points[i] = test_points[i-1] + dist(gen);
-        }
+    // Test with extremely small values
+    double small_E = 1.0e-20;
+    double result_small_bs = bs_tests.interpolateBS(small_E);
+    double result_small_dl = dl_tests.interpolateDL(small_E);
 
-        // Measure binary search performance
-        {
-            DoubleLookupTests tests;
-            auto start = std::chrono::high_resolution_clock::now();
-            for (const double& E : test_points) {
-                volatile double result = interpolate_binary_search_cpp(E, tests);
-            }
-            auto end = std::chrono::high_resolution_clock::now();
-            auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(
-                end - start).count();
-            std::cout << "Binary Search (vector) time: " << duration << " ns" << std::endl;
-        }
+    // For extreme values, we should get boundary values
+    assert(is_equal(result_large_bs, 3278.15));
+    assert(is_equal(result_large_dl, 3273.15));
+    assert(is_equal(result_small_bs, 3243.15));
+    assert(is_equal(result_small_dl, 3243.15));
 
-        {
-            DoubleLookupTests tests;
-            auto start = std::chrono::high_resolution_clock::now();
-            for (const double& E : test_points) {
-                volatile double result = interpolate_binary_search_cpp(E, tests);
-            }
-            auto end = std::chrono::high_resolution_clock::now();
-            auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(
-                end - start).count();
-            std::cout << "Binary Search (array) time: " << duration << " ns" << std::endl;
+    std::cout << "Stress tests passed" << std::endl;
+}
+
+void test_random_validation() {
+    std::cout << "\nTesting with random inputs..." << std::endl;
+
+    BinarySearchTests bs_tests;
+    DoubleLookupTests dl_tests;
+
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<double> dist(16800000000.0, 17100000000.0);
+
+    const int num_tests = 1000;
+    for (int i = 0; i < num_tests; ++i) {
+        double random_E = dist(gen);
+
+        double result_bs = bs_tests.interpolateBS(random_E);
+        double result_dl = dl_tests.interpolateDL(random_E);
+
+        // Results should be within the temperature range
+        assert(result_bs >= 3243.15 && result_bs <= 3278.15);
+        assert(result_dl >= 3243.15 && result_dl <= 3273.15);
+    }
+
+    std::cout << "Random input validation passed" << std::endl;
+}
+
+
+void test_performance() {
+    std::cout << "\nTesting performance..." << std::endl;
+
+    // Configuration parameters
+    constexpr int warmupSteps = 5;
+    constexpr int outerIterations = 10;
+    constexpr int numCells = 64*64*64;
+
+    // Setup test data
+    SS316L test;
+    std::vector<double> random_energies(numCells);
+
+    // Generate random values
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    const double E_min = SS316L::E_neq.front() * 0.8;
+    const double E_max = SS316L::E_neq.back() * 1.2;
+    std::uniform_real_distribution<double> dist(E_min, E_max);
+
+    // Fill random energies
+    for(auto& E : random_energies) {
+        E = dist(gen);
+    }
+
+    // Warmup runs
+    std::cout << "Performing warmup steps..." << std::endl;
+    for(int i = 0; i < warmupSteps; ++i) {
+        for(const double& E : random_energies) {
+            // volatile double result = interpolate_double_lookup_cpp(E, test);
+            volatile double result = test.interpolateDL(E);
+        }
+    }
+    for(int i = 0; i < warmupSteps; ++i) {
+        for(const double& E : random_energies) {
+            volatile double result = test.interpolateBS(E);
         }
+    }
 
-        // Measure double lookup performance
+    // Performance measurement
+    std::cout << "\nStarting performance measurement..." << std::endl;
+    std::vector<double> timings_binary;
+    std::vector<double> timings_double_lookup;
+
+    for(int iter = 0; iter < outerIterations; ++iter) {
+        std::cout << "\nIteration " << iter + 1 << "/" << outerIterations << std::endl;
+
+        // Double Lookup timing
         {
-            DoubleLookupTests tests;
-            auto start = std::chrono::high_resolution_clock::now();
-            for (const double& E : test_points) {
-                volatile double result = tests.interpolateDL(E);
+            const auto start1 = std::chrono::high_resolution_clock::now();
+            for(const double& E : random_energies) {
+                volatile double result = test.interpolateDL(E);
             }
-            auto end = std::chrono::high_resolution_clock::now();
-            auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(
-                end - start).count();
-            std::cout << "Double Lookup (vector) time: " << duration << " ns" << std::endl;
+            const auto end1 = std::chrono::high_resolution_clock::now();
+            const auto duration1 = std::chrono::duration_cast<std::chrono::nanoseconds>(
+                end1 - start1).count();
+            timings_double_lookup.push_back(static_cast<double>(duration1));
+
+            std::cout << "Double Lookup - Iteration time: " << duration1 << " ns" << std::endl;
         }
 
+        // Binary Search timing
         {
-            auto start = std::chrono::high_resolution_clock::now();
-            for (const double& E : test_points) {
-                DoubleLookupTests tests;
-                volatile double result = tests.interpolateDL(E);
+            const auto start2 = std::chrono::high_resolution_clock::now();
+            for(const double& E : random_energies) {
+                volatile double result = test.interpolateBS(E);
             }
-            auto end = std::chrono::high_resolution_clock::now();
-            auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(
-                end - start).count();
-            std::cout << "Double Lookup (array) time: " << duration << " ns" << std::endl;
+            const auto end2 = std::chrono::high_resolution_clock::now();
+            const auto duration2 = std::chrono::duration_cast<std::chrono::nanoseconds>(
+                end2 - start2).count();
+            timings_binary.push_back(static_cast<double>(duration2));
+
+            std::cout << "Binary Search - Iteration time: " << duration2 << " ns" << std::endl;
         }
     }
+
+    // Calculate and print statistics
+    auto calc_stats = [](const std::vector<double>& timings) {
+        double sum = std::accumulate(timings.begin(), timings.end(), 0.0);
+        double mean = sum / static_cast<double>(timings.size());
+        double sq_sum = std::inner_product(timings.begin(), timings.end(),
+                                         timings.begin(), 0.0);
+        double stdev = std::sqrt(sq_sum / static_cast<double>(timings.size()) - mean * mean);
+        return std::make_pair(mean, stdev);
+    };
+
+    auto [binary_mean, binary_stdev] = calc_stats(timings_binary);
+    auto [lookup_mean, lookup_stdev] = calc_stats(timings_double_lookup);
+
+    std::cout << "\nPerformance Results (" << numCells << " cells, "
+              << outerIterations << " iterations):" << std::endl;
+    std::cout << "Binary Search:" << std::endl;
+    std::cout << "  Mean time: " << binary_mean << " ± " << binary_stdev << " ns" << std::endl;
+    std::cout << "  Per cell: " << binary_mean/numCells << " ns" << std::endl;
+
+    std::cout << "Double Lookup:" << std::endl;
+    std::cout << "  Mean time: " << lookup_mean << " ± " << lookup_stdev << " ns" << std::endl;
+    std::cout << "  Per cell: " << lookup_mean/numCells << " ns" << std::endl;
 }
 
 
@@ -607,6 +291,9 @@ int main() {
         test_basic_functionality();
         test_edge_cases_and_errors();
         test_interpolation_accuracy();
+        test_consistency();
+        test_stress();
+        test_random_validation();
         test_performance();
 
         std::cout << "\nAll tests completed successfully!" << std::endl;