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 d1f3b6a686fcbdc08ea841849f25ae5bcb8f78a1..a5770a1c6a0d5697665417ec34fabeb95065e753 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,36 +5,36 @@
 
 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
-    /*if (arrs.E_neq[idx_E_neq + 1] < E_target) {
-        ++idx_E_neq;
+    /*if (arrs.y_neq[idx_y_neq + 1] < y_target) {
+        ++idx_y_neq;
     }*/
-    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);
 }