diff --git a/final_plots.py b/final_plots.py
index d770e7c47829d6ae69faca085102fc9408f08149..b48fab87c5e937539a8ca822516c0040ea1dd9ec 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -24,7 +24,14 @@ from numbers import Number
 # temperature Tkrg, which is an integration constant of the E-flow RG
 # equations. The more conventional definition of the Kondo temperature is
 # G(V=Tk)=G(V=0)/2=e²/h. The ratio Tk/Tkrg is:
-TK_VOLTAGE = 3.44249
+
+#TK_VOLTAGE = 3.44249 # for rtol=1e-10, atol=1e-12, voltage_branches=10
+TK_VOLTAGE = 3.4425351 # for rtol=1e-9, atol=1e-11, voltage_branches=10
+#TK_VOLTAGE = 3.44334 # with correction
+TK_VOLTAGE_O3P = 3.458524 # for rtol=1e-9, atol=1e-11, voltage_branches=10
+#TK_VOLTAGE_O3P = 3.4593 # with correction
+TK_VOLTAGE_ORDER2 = 10.1368086 # for rtol=1e-9, atol=1e-11, voltage_branches=10
+#TK_VOLTAGE_ORDER2 = 10.13754 # with correction
 
 def save_overview(
         omega = 16.5372,
@@ -250,6 +257,9 @@ def export_interp(
         y_func = None,
         korder = 2,
         special_selection = None,
+        o2_scale_x = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
+        o2_scale_y = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
+        o2_scale_fixed = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
         **parameters
         ):
     """
@@ -330,28 +340,35 @@ def export_interp(
             o3p = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
             )
 
+    o2_data = dm.list(**{fixed_parameter:fixed_value*o2_scale_fixed}, **parameters, min_version=(14,15,-1,-1))
+    o2_data = o2_data.loc[o2_data.solver_flags & (global_bad_flags | bad_flags["o2"] | required_flags["o2"]) == required_flags["o2"]]
+    o2_data = o2_data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
+
     def selection(order, method, padding=False, **kwargs):
-        selected = (data.solver_flags & required_flags[order] == required_flags[order]) \
-                & (data.solver_flags & bad_flags[order] == 0) \
-                & (data.method == method)
+        thedata = o2_data if order == "o2" else data
+        selected = (thedata.solver_flags & required_flags[order] == required_flags[order]) \
+                & (thedata.solver_flags & bad_flags[order] == 0) \
+                & (thedata.method == method)
         if padding:
-            selected &= (data["padding"] > 0) | (data["nmax"] == 0)
+            selected &= (thedata["padding"] > 0) | (thedata["nmax"] == 0)
         else:
-            selected &= data["padding"] == 0
+            selected &= thedata["padding"] == 0
         for key, value in kwargs.items():
-            selected &= (data[key] == value)
+            selected &= (thedata[key] == value)
         if special_selection is not None:
-            selected &= special_selection(data, order, method, padding, **kwargs)
-        return selected
+            selected &= special_selection(thedata, order, method, padding, **kwargs)
+        return thedata.loc[selected]
 
     def gen_data(order, method, padding=False, s=5e-5, **kwargs):
         suffix = "_p" if padding else ""
-        reduced_data = data.loc[selection(order, method, padding, **kwargs)]
+        reduced_data = selection(order, method, padding, **kwargs)
         settings.logger.info(f"Starting with {order}{suffix} {method}, found {reduced_data.shape[0]} data points")
         results = {}
         if reduced_data.shape[0] < 100:
             return results
         xy_data = (x_func(reduced_data), y_func(reduced_data))
+        if order == "o2":
+            xy_data = (xy_data[0]/o2_scale_x, xy_data[1]/o2_scale_y)
         try:
             results[f"gdc_{method}_{order}{suffix}"] = griddata(
                     xy_data,
diff --git a/plot.py b/plot.py
index 17deef0492efd785d13a35674843488a6f325462..f6e3ffdaddbbd9639c3752daeb2996aab7ec1411 100644
--- a/plot.py
+++ b/plot.py
@@ -38,7 +38,13 @@ REFERENCE_VALUES = {
             ),
         }
 
-TK_VOLTAGE = 3.44249
+#TK_VOLTAGE = 3.44249 # for rtol=1e-10, atol=1e-12, voltage_branches=10
+TK_VOLTAGE = 3.4425351 # for rtol=1e-9, atol=1e-11, voltage_branches=10
+#TK_VOLTAGE = 3.44334 # with correction
+TK_VOLTAGE_O3P = 3.458524 # for rtol=1e-9, atol=1e-11, voltage_branches=10
+#TK_VOLTAGE_O3P = 3.4593 # with correction
+TK_VOLTAGE_ORDER2 = 10.1368086 # for rtol=1e-9, atol=1e-11, voltage_branches=10
+#TK_VOLTAGE_ORDER2 = 10.13754 # with correction
 
 
 def main():
@@ -1568,8 +1574,8 @@ def compare_RGeq_interp(
         levels = np.linspace(-np.pi, np.pi, 41)
     else:
         raise ValueError(f"unknown observable: {observable}")
-    vdc = data["vdc"]
-    vac = data["vac"]
+    vdc = data["vdc"][0]
+    vac = data["vac"][:,0]
     extent = (1.5*vdc[0] - 0.5*vdc[1], 1.5*vdc[-1] - 0.5*vdc[-2], 1.5*vac[0] - 0.5*vac[1], 1.5*vac[-1] - 0.5*vac[-2])
     imshow_kwargs = dict(cmap=plt.cm.seismic, extent=extent, origin="lower", aspect="auto", zorder=0)
     contour_kwargs = dict(colors="#000000", alpha=0.7, linewidths=0.5, zorder=2)