From 2cd79231db0559c2840ccc0afa799609923d73f3 Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Fri, 17 Feb 2023 14:57:33 +0100
Subject: [PATCH] export omega5 interpolation: use new initial conditions

---
 final_plots.py                                | 209 +++++++++++-------
 plot.py                                       |  31 +++
 tikz/current_convergence_all.tex              | 179 +++++++++++++++
 tikz/current_convergence_bad.tex              |  16 +-
 tikz/current_convergence_good.tex             | 114 ++++++++++
 ...rgence.tex => current_convergence_old.tex} |  10 +-
 6 files changed, 461 insertions(+), 98 deletions(-)
 create mode 100644 tikz/current_convergence_all.tex
 create mode 100644 tikz/current_convergence_good.tex
 rename tikz/{current_convergence.tex => current_convergence_old.tex} (92%)

diff --git a/final_plots.py b/final_plots.py
index 8d87f84..95ce3b3 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -317,6 +317,8 @@ def export_interp(
             o3a = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             o3p = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
             )
+    i_flags_good = DataManager.SOLVER_FLAGS["improved_initial_conditions"]
+    g_flags_bad = DataManager.SOLVER_FLAGS["improved_initial_conditions"]
 
     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"]]
@@ -348,7 +350,6 @@ def export_interp(
         results = {}
         if reduced_data.shape[0] < 100:
             return results
-        settings.logger.info(f"Starting with {order}{suffix} {method}, found {reduced_data.shape[0]} data points")
         if extend_vdc or extend_vac:
             extended_data = [reduced_data]
             if extend_vdc:
@@ -370,62 +371,78 @@ def export_interp(
                     del acdc_mirrored
                 del ac_mirrored
             reduced_data = pd.concat(extended_data)
-        xy_data = (x_func(reduced_data), y_func(reduced_data))
+        g_data = reduced_data.loc[reduced_data.solver_flags & g_flags_bad == 0]
+        i_data = reduced_data.loc[reduced_data.solver_flags & i_flags_good == i_flags_good]
+        settings.logger.info(f"Starting with {order}{suffix} {method}, using {g_data.shape[0]}/{i_data.shape[0]} data points")
+        xy_data = (x_func(g_data), y_func(g_data))
+        has_current = i_data.shape[0] >= 100
+        if has_current:
+            xy_data_i = (x_func(i_data), y_func(i_data))
         if order == "o2":
             xy_data = (xy_data[0]/o2_scale_x, xy_data[1]/o2_scale_y)
+            if has_current:
+                xy_data_i = (xy_data_i[0]/o2_scale_x, xy_data_i[1]/o2_scale_y)
         elif order == "o3r":
             xy_data = (xy_data[0]/o3r_scale_x, xy_data[1]/o3r_scale_y)
+            if has_current:
+                xy_data_i = (xy_data_i[0]/o3r_scale_x, xy_data_i[1]/o3r_scale_y)
         try:
             results[f"gdc_{method}_{order}{suffix}"] = griddata(
                     xy_data,
-                    reduced_data.dc_conductance,
+                    g_data.dc_conductance,
                     mesh,
                     method="cubic")
-            #gdc_tck = bisplrep(reduced_data[y], reduced_data[x], reduced_data.dc_conductance, s=s, kx=kxorder, ky=kyorder)
+            #gdc_tck = bisplrep(g_data[y], g_data[x], g_data.dc_conductance, s=s, kx=kxorder, ky=kyorder)
             #results[f"gdc_{method}_{order}{suffix}"] = bisplev(y_arr, x_arr, gdc_tck)
         except:
             settings.logger.exception(f"in griddata interpolation of Gdc for {order}{suffix} {method}")
         try:
-            results[f"idc_{method}_{order}{suffix}"] = griddata(
-                    xy_data,
-                    reduced_data.dc_current,
-                    mesh,
-                    method="cubic")
-        except:
-            settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method}")
-        try:
-            idc_tck = bisplrep(*xy_data[::-1], reduced_data.dc_current, s=s, kx=kxorder, ky=kyorder)
-            results[f"idc_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, idc_tck)
-            results[f"gdc_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, idc_tck, dy=1)
-        except:
-            settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method}")
-        try:
-            results[f"iac_{method}_{order}{suffix}"] = griddata(
-                    xy_data,
-                    2*reduced_data.ac_current_abs,
-                    mesh,
-                    method="cubic")
-        except:
-            settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
-        try:
-            iac_tck = bisplrep(*xy_data[::-1], 2*reduced_data.ac_current_abs, s=s, kx=kxorder, ky=kyorder)
-            results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, iac_tck)
-            results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
+            idc_tck_old = bisplrep(*xy_data[::-1], g_data.dc_current, s=s, kx=kxorder, ky=kyorder)
+            results[f"gdc_{method}_{order}{suffix}_ispline_old"] = bisplev(y_arr, x_arr, idc_tck_old, dy=1)
         except:
-            settings.logger.exception(f"in spline interpolation of Iac for {order}{suffix} {method}")
-        try:
-            sel = reduced_data.vac > 0
-            results[f"phase_{method}_{order}{suffix}"] = griddata(
-                    (xy_data[0][sel], xy_data[1][sel]),
-                    reduced_data.ac_current_phase[sel],
-                    mesh,
-                    method="cubic")
-        except:
-            settings.logger.exception(f"in griddata interpolation of phase for {order}{suffix} {method}")
+            settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method} (old)")
+        if has_current:
+            try:
+                results[f"idc_{method}_{order}{suffix}"] = griddata(
+                        xy_data_i,
+                        i_data.dc_current,
+                        mesh,
+                        method="cubic")
+            except:
+                settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method}")
+            try:
+                idc_tck = bisplrep(*xy_data_i[::-1], i_data.dc_current, s=s, kx=kxorder, ky=kyorder)
+                results[f"idc_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, idc_tck)
+                results[f"gdc_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, idc_tck, dy=1)
+            except:
+                settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method}")
+            try:
+                results[f"iac_{method}_{order}{suffix}"] = griddata(
+                        xy_data_i,
+                        2*i_data.ac_current_abs,
+                        mesh,
+                        method="cubic")
+            except:
+                settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
+            try:
+                iac_tck = bisplrep(*xy_data_i[::-1], 2*i_data.ac_current_abs, s=s, kx=kxorder, ky=kyorder)
+                results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, iac_tck)
+                results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
+            except:
+                settings.logger.exception(f"in spline interpolation of Iac for {order}{suffix} {method}")
+            try:
+                sel = i_data.vac > 0
+                results[f"phase_{method}_{order}{suffix}"] = griddata(
+                        (xy_data_i[0][sel], xy_data_i[1][sel]),
+                        i_data.ac_current_phase[sel],
+                        mesh,
+                        method="cubic")
+            except:
+                settings.logger.exception(f"in griddata interpolation of phase for {order}{suffix} {method}")
         try:
             results[f"gamma_{method}_{order}{suffix}"] = griddata(
                     xy_data,
-                    reduced_data.gamma,
+                    g_data.gamma,
                     mesh,
                     method="cubic")
         except:
@@ -550,6 +567,9 @@ def export_interp_relative_o3a(
             o3p = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
             o3a_vb7 = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             )
+    i_flags_good = DataManager.SOLVER_FLAGS["improved_initial_conditions"]
+    g_flags_bad = DataManager.SOLVER_FLAGS["improved_initial_conditions"]
+
     voltage_branches = dict(
             o3 = 4,
             o3a = 4,
@@ -582,7 +602,6 @@ def export_interp_relative_o3a(
     def gen_data(order, method, padding=False, s=spline_s, extend_vdc=10, extend_vac=10, **kwargs):
         suffix = "_p" if padding else ""
         reduced_data = selection(order, method, padding, voltage_branches[order], **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
@@ -610,56 +629,63 @@ def export_interp_relative_o3a(
                     del acdc_mirrored
                 del ac_mirrored
             reduced_data = pd.concat(extended_data)
-        xy_data = (x_func(reduced_data), y_func(reduced_data))
+        g_data = reduced_data.loc[(reduced_data.solver_flags_ref | reduced_data.solver_flags) & g_flags_bad == 0]
+        i_data = reduced_data.loc[reduced_data.solver_flags_ref & reduced_data.solver_flags & i_flags_good == i_flags_good]
+        settings.logger.info(f"Starting with {order}{suffix} {method}, using {g_data.shape[0]}/{i_data.shape[0]} data points")
+        xy_data = (x_func(g_data), y_func(g_data))
+        has_current = i_data.shape[0] >= 100
+        if has_current:
+            xy_data_i = (x_func(i_data), y_func(i_data))
         try:
             results[f"gdc_{method}_{order}{suffix}"] = griddata(
                     xy_data,
-                    reduced_data.dc_conductance - reduced_data.dc_conductance_ref,
+                    g_data.dc_conductance - g_data.dc_conductance_ref,
                     mesh,
                     method="cubic")
         except:
             settings.logger.exception(f"in griddata interpolation of Gdc for {order}{suffix} {method}")
-        try:
-            results[f"idc_{method}_{order}{suffix}"] = griddata(
-                    xy_data,
-                    reduced_data.dc_current - reduced_data.dc_current_ref,
-                    mesh,
-                    method="cubic")
-        except:
-            settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method}")
-        try:
-            idc_tck = bisplrep(*xy_data[::-1], reduced_data.dc_current - reduced_data.dc_current_ref, s=s, kx=kxorder, ky=kyorder)
-            results[f"idc_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, idc_tck)
-            results[f"gdc_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, idc_tck, dy=1)
-        except:
-            settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method}")
-        try:
-            results[f"iac_{method}_{order}{suffix}"] = griddata(
-                    xy_data,
-                    2*(reduced_data.ac_current_abs - reduced_data.ac_current_abs_ref),
-                    mesh,
-                    method="cubic")
-        except:
-            settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
-        try:
-            iac_tck = bisplrep(*xy_data[::-1], 2*(reduced_data.ac_current_abs - reduced_data.ac_current_abs_ref), s=s, kx=kxorder, ky=kyorder)
-            results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, iac_tck)
-            results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
-        except:
-            settings.logger.exception(f"in spline interpolation of Iac for {order}{suffix} {method}")
-        try:
-            sel = reduced_data.vac > 0
-            results[f"phase_{method}_{order}{suffix}"] = griddata(
-                    (xy_data[0][sel], xy_data[1][sel]),
-                    reduced_data.ac_current_phase[sel] - reduced_data.ac_current_phase_ref[sel],
-                    mesh,
-                    method="cubic")
-        except:
-            settings.logger.exception(f"in griddata interpolation of phase for {order}{suffix} {method}")
+        if has_current:
+            try:
+                results[f"idc_{method}_{order}{suffix}"] = griddata(
+                        xy_data_i,
+                        i_data.dc_current - i_data.dc_current_ref,
+                        mesh,
+                        method="cubic")
+            except:
+                settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method}")
+            try:
+                idc_tck = bisplrep(*xy_data_i[::-1], i_data.dc_current - i_data.dc_current_ref, s=s, kx=kxorder, ky=kyorder)
+                results[f"idc_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, idc_tck)
+                results[f"gdc_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, idc_tck, dy=1)
+            except:
+                settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method}")
+            try:
+                results[f"iac_{method}_{order}{suffix}"] = griddata(
+                        xy_data_i,
+                        2*(i_data.ac_current_abs - i_data.ac_current_abs_ref),
+                        mesh,
+                        method="cubic")
+            except:
+                settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
+            try:
+                iac_tck = bisplrep(*xy_data_i[::-1], 2*(i_data.ac_current_abs - i_data.ac_current_abs_ref), s=s, kx=kxorder, ky=kyorder)
+                results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, iac_tck)
+                results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
+            except:
+                settings.logger.exception(f"in spline interpolation of Iac for {order}{suffix} {method}")
+            try:
+                sel = i_data.vac > 0
+                results[f"phase_{method}_{order}{suffix}"] = griddata(
+                        (xy_data_i[0][sel], xy_data_i[1][sel]),
+                        i_data.ac_current_phase[sel] - i_data.ac_current_phase_ref[sel],
+                        mesh,
+                        method="cubic")
+            except:
+                settings.logger.exception(f"in griddata interpolation of phase for {order}{suffix} {method}")
         try:
             results[f"gamma_{method}_{order}{suffix}"] = griddata(
                     xy_data,
-                    reduced_data.gamma - reduced_data.gamma_ref,
+                    g_data.gamma - g_data.gamma_ref,
                     mesh,
                     method="cubic")
         except:
@@ -828,6 +854,8 @@ def prepare_RGconvergence():
     real_cmap = lambda x: seismic((x+1)/2)[...,:3]
     imwrite(f"figdata/colorbar_seismic.png", np.array(0xffff*seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((1,0x1000,3)), format="PNG-FI")
     imwrite(f"figdata/colorbar_seismic_vertical.png", np.array(0xffff*seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((0x1000,1,3)), format="PNG-FI")
+    imwrite(f"figdata/colorbar_viridis_horizontal.png", np.array(0xffff*viridis(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((1,0x1000,3)), format="PNG-FI")
+    imwrite(f"figdata/colorbar_viridis.png", np.array(0xffff*viridis(np.linspace(0, 1, 0x1000))[::-1,:3], dtype=np.uint16).reshape((0x1000,1,3)), format="PNG-FI")
     def export_img(name, array):
         imwrite(f"figdata/{name}.png", np.array(0xffff*real_cmap(array[::-1]), dtype=np.uint16), format="PNG-FI")
 
@@ -1633,7 +1661,7 @@ def export_harmonic_modes(
     plt.show()
 
 
-def export_convergence_pgfplots(simplified_initial_conditions = False):
+def export_convergence_pgfplots(improved_initial_conditions=True, simplified_initial_conditions=False):
     """
     Plot convergence of current as function of D (Λ₀) for a fine
     grid of dc and ac voltage at fixed omega.
@@ -1653,17 +1681,29 @@ def export_convergence_pgfplots(simplified_initial_conditions = False):
                     | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
                     | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
                     | DataManager.SOLVER_FLAGS["deleted"],
-            good_flags = DataManager.SOLVER_FLAGS["include_Ga"],
+            good_flags = DataManager.SOLVER_FLAGS["include_Ga"] \
+                    | DataManager.SOLVER_FLAGS["improved_initial_conditions"]
                 )
     if simplified_initial_conditions:
         parameters["good_flags"] ^= DataManager.SOLVER_FLAGS["simplified_initial_conditions"]
         parameters["bad_flags"] ^= DataManager.SOLVER_FLAGS["simplified_initial_conditions"]
+    if not improved_initial_conditions:
+        parameters["good_flags"] ^= DataManager.SOLVER_FLAGS["improved_initial_conditions"]
+        parameters["bad_flags"] ^= DataManager.SOLVER_FLAGS["improved_initial_conditions"]
 
     vac_num = 26
     vdc_num = 26
     vdc_max = 82.686
     vac_max = 82.686
-    exponent = 2 if parameters["good_flags"] & DataManager.SOLVER_FLAGS["simplified_initial_conditions"] else 3
+    if simplified_initial_conditions:
+        exponent = 2
+        suffix = "_bad"
+    elif improved_initial_conditions:
+        exponent = 4
+        suffix = "_good"
+    else:
+        exponent = 3
+        suffix = "_old"
     d_num = 5
     log10d_min = 5
     log10d_max = 9
@@ -1699,7 +1739,6 @@ def export_convergence_pgfplots(simplified_initial_conditions = False):
     #print(f"Extent: xmin={extent_dc[0]/TK_VOLTAGE:.6g}, xmax={extent_dc[1]/TK_VOLTAGE:,.6g}, ymin={extent_dc[2]/TK_VOLTAGE:.6g}, ymax={extent_dc[3]/TK_VOLTAGE:,.6g}")
     #print(f"Extent: xmin={extent_ac[0]/TK_VOLTAGE:.6g}, xmax={extent_ac[1]/TK_VOLTAGE:,.6g}, ymin={extent_ac[2]/TK_VOLTAGE:.6g}, ymax={extent_ac[3]/TK_VOLTAGE:,.6g}")
 
-    suffix = "_bad" if simplified_initial_conditions else "_good"
     def export_img(name, array):
         imwrite(f"figdata/{name}{suffix}.png", np.array(0xffff*viridis(array[::-1])[...,:3], dtype=np.uint16), format="PNG-FI")
     idc_max = np.nanmax(idc_fit_a)
diff --git a/plot.py b/plot.py
index 96a0c45..c6b6287 100644
--- a/plot.py
+++ b/plot.py
@@ -2489,5 +2489,36 @@ def compare_resonant(dm, **parameters):
     plt.show()
 
 
+def compare_current_initial_conditions(dm, **parameters):
+    data = dm.list(**parameters)
+    global_bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
+            | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
+            | DataManager.SOLVER_FLAGS["deleted"]
+    data = data.loc[data.solver_flags & global_bad_flags == 0]
+    data = data.sort_values(["vdc", "vac"])
+    data.vdc = np.round(data.vdc, 8)
+    data.vac = np.round(data.vac, 8)
+    cmp_good_flags = DataManager.SOLVER_FLAGS["include_Ga"]
+    cmp_bad_flags = DataManager.SOLVER_FLAGS["improved_initial_conditions"] | DataManager.SOLVER_FLAGS["simplified_initial_conditions"]
+    ref_good_flags = DataManager.SOLVER_FLAGS["improved_initial_conditions"] | DataManager.SOLVER_FLAGS["include_Ga"]
+    ref_bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"]
+    cmp_flags = cmp_good_flags | cmp_bad_flags
+    ref_flags = ref_good_flags | ref_bad_flags
+
+    ref_data = data.loc[data.solver_flags & ref_flags == ref_good_flags]
+    cmp_data = data.loc[data.solver_flags & cmp_flags == cmp_good_flags]
+    merged = pd.merge(
+            data,
+            ref_data,
+            how="inner",
+            on=["vdc","vac","d","omega","energy_im","energy_re","solver_tol_rel","solver_tol_abs","xL","lazy_inverse_factor","resonant_dc_shift"],
+            suffixes=("", "_ref"))
+
+    s = plt.scatter(merged.vdc, merged.vac, c=(merged.dc_current_ref-merged.dc_current)/merged.dc_current_ref, cmap=plt.cm.seismic)
+    plt.colorbar(s)
+    plt.show()
+
+
+
 if __name__ == '__main__':
     main()
diff --git a/tikz/current_convergence_all.tex b/tikz/current_convergence_all.tex
new file mode 100644
index 0000000..8857e16
--- /dev/null
+++ b/tikz/current_convergence_all.tex
@@ -0,0 +1,179 @@
+\tikzsetnextfilename{current_convergence_all}%
+\newcommand\imgsize{32mm}%
+\begin{tikzpicture}
+  \begin{axis}[
+      name = idc,
+      title = $I_\mathrm{avg}$,
+      title style = {yshift = -2mm},
+      height = \imgsize,
+      width = \imgsize,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xticklabels = {},
+      ylabel = $\vac~(k_B \tkv/e)$,
+      axis on top,
+    ]
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/idc_converged_good.png};
+  \end{axis}
+  \begin{axis}[
+      name = iac,
+      at = (idc.south),
+      anchor = north,
+      title = $I_\mathrm{osc}$,
+      title style = {yshift = -2mm},
+      yshift = -8mm,
+      height = \imgsize,
+      width = \imgsize,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xlabel = $\vdc~(k_B \tkv/e)$,
+      ylabel = $\vac~(k_B \tkv/e)$,
+      axis on top,
+    ]
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/iac_converged_good.png};
+  \end{axis}
+  %
+  \begin{axis}[
+      name = idc_diff,
+      at = (idc.east),
+      anchor = west,
+      xshift = 4mm,
+      %title = {$|I_\mathrm{avg}(\Lambda_0=10^9\tkrg) - I_\mathrm{avg}|/I_\mathrm{avg}$},
+      %title = {relative deviation at $\Lambda_0=10^9\tkrg$},
+      title = {good initial condition},
+      title style = {yshift = -2mm},
+      height = \imgsize,
+      width = \imgsize,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xticklabels = {},
+      yticklabels = {},
+      axis on top,
+    ]
+    \addplot graphics[xmin=0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/idc_diff_good.png};
+  \end{axis}
+  \begin{axis}[
+      name = iac_diff,
+      at = (iac.east),
+      anchor = west,
+      xshift = 4mm,
+      %title = {$|I_\mathrm{osc}(\Lambda_0=10^9\tkrg) - I_\mathrm{osc}|/I_\mathrm{osc}$},
+      %title style = {yshift = -2mm},
+      height = \imgsize,
+      width = \imgsize,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xlabel = $\vdc~(k_B \tkv/e)$,
+      yticklabels = {},
+      axis on top,
+    ]
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=0.480379, ymax=24.4993] {../figdata/iac_diff_good.png};
+  \end{axis}
+  %
+  \begin{axis}[
+      name = idc_diff2,
+      at = (idc_diff.east),
+      anchor = west,
+      xshift = 4mm,
+      title = {bad initial condition},
+      title style = {yshift = -1.2mm},
+      height = \imgsize,
+      width = \imgsize,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xticklabels = {},
+      yticklabels = {},
+      axis on top,
+    ]
+    \addplot graphics[xmin=0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/idc_diff_bad.png};
+  \end{axis}
+  \begin{axis}[
+      name = iac_diff2,
+      at = (iac_diff.east),
+      anchor = west,
+      xshift = 4mm,
+      title style = {yshift = -2mm},
+      height = \imgsize,
+      width = \imgsize,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xlabel = $\vdc~(k_B \tkv/e)$,
+      yticklabels = {},
+      axis on top,
+    ]
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=0.480379, ymax=24.4993] {../figdata/iac_diff_bad.png};
+  \end{axis}
+  %
+  \begin{axis}[
+      name = cb1,
+      at = (iac.south),
+      anchor = north,
+      xmin = 0,
+      xmax = 1.54839,
+      yshift = -17mm,
+      enlargelimits = false,
+      axis on top,
+      height = .75em,
+      width = \imgsize,
+      scale only axis,
+      tickpos = bottom,
+      tick align = outside,
+      ytick = \empty,
+      xlabel = {current $(e k_B \tkv/\hbar)$},
+      xlabel shift = -14.25mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=0, xmax=1.54775] {../figdata/colorbar_viridis_horizontal.png};
+  \end{axis}
+  \begin{axis}[
+      name = cb2,
+      at = (iac_diff.south),
+      anchor = north,
+      xmin = 0,
+      xmax = 0.000310975,
+      yshift = -17mm,
+      enlargelimits = false,
+      axis on top,
+      height = .75em,
+      width = \imgsize,
+      scale only axis,
+      tickpos = bottom,
+      tick align = outside,
+      ytick = \empty,
+      xlabel = {relative deviation},
+      scaled ticks = false,
+      xtick = {0, 2e-4},
+      xlabel shift = -15mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=0, xmax=0.000310975] {../figdata/colorbar_viridis_horizontal.png};
+  \end{axis}
+  \begin{axis}[
+      name = cb3,
+      at = (iac_diff2.south),
+      anchor = north,
+      xmin = 0,
+      xmax = 0.0534487,
+      yshift = -17mm,
+      enlargelimits = false,
+      axis on top,
+      height = .75em,
+      width = \imgsize,
+      scale only axis,
+      tickpos = bottom,
+      tick align = outside,
+      ytick = \empty,
+      xlabel = {relative deviation},
+      xtick = {0, 0.02, 0.04},
+      xticklabels = {$0$, $2\%$, $4\%$},
+      scaled ticks = false,
+      xlabel shift = -14.17mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=0, xmax=0.0534487] {../figdata/colorbar_viridis_horizontal.png};
+  \end{axis}
+  \node[xshift=2mm, yshift=9mm] at (idc_diff.north east) {relative deviation at $\Lambda_0=10^9\tkrg$ from $\Lambda_0=\infty$};
+\end{tikzpicture}
diff --git a/tikz/current_convergence_bad.tex b/tikz/current_convergence_bad.tex
index 7c554a5..9102afa 100644
--- a/tikz/current_convergence_bad.tex
+++ b/tikz/current_convergence_bad.tex
@@ -7,8 +7,8 @@
       height = 4.5cm,
       width = 4.5cm,
       scale only axis,
-      xmin = -0.5, xmax = 24.5,
-      ymin = -0.5, ymax = 24.5,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
       xticklabels = {},
       ylabel = $\vac~(k_B \tkv/e)$,
       axis on top,
@@ -25,8 +25,8 @@
       height = 4.5cm,
       width = 4.5cm,
       scale only axis,
-      xmin = -0.5, xmax = 24.5,
-      ymin = -0.5, ymax = 24.5,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
       xticklabels = {},
       yticklabels = {},
       axis on top,
@@ -43,8 +43,8 @@
       height = 4.5cm,
       width = 4.5cm,
       scale only axis,
-      xmin = -0.5, xmax = 24.5,
-      ymin = -0.5, ymax = 24.5,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
       xlabel = $\vdc~(k_B \tkv/e)$,
       ylabel = $\vac~(k_B \tkv/e)$,
       axis on top,
@@ -61,8 +61,8 @@
       height = 4.5cm,
       width = 4.5cm,
       scale only axis,
-      xmin = -0.5, xmax = 24.5,
-      ymin = -0.5, ymax = 24.5,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
       xlabel = $\vdc~(k_B \tkv/e)$,
       yticklabels = {},
       axis on top,
diff --git a/tikz/current_convergence_good.tex b/tikz/current_convergence_good.tex
new file mode 100644
index 0000000..9a3f0ab
--- /dev/null
+++ b/tikz/current_convergence_good.tex
@@ -0,0 +1,114 @@
+\tikzsetnextfilename{current_convergence_good}%
+\begin{tikzpicture}
+  \begin{axis}[
+      name = idc,
+      title = $I_\mathrm{avg}$,
+      title style = {yshift = -2mm},
+      height = 4.5cm,
+      width = 4.5cm,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xticklabels = {},
+      ylabel = $\vac~(k_B \tkv/e)$,
+      axis on top,
+    ]
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/idc_converged_good.png};
+  \end{axis}
+  \begin{axis}[
+      name = iac,
+      at = (idc.east),
+      anchor = west,
+      title = $I_\mathrm{osc}$,
+      title style = {yshift = -2mm},
+      xshift = 4mm,
+      height = 4.5cm,
+      width = 4.5cm,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xticklabels = {},
+      yticklabels = {},
+      axis on top,
+    ]
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/iac_converged_good.png};
+  \end{axis}
+  \begin{axis}[
+      name = idc_diff,
+      at = (idc.south),
+      anchor = north,
+      yshift = -9mm,
+      title = {$|I_\mathrm{avg}(\Lambda_0=10^9\tkrg) - I_\mathrm{avg}|/I_\mathrm{avg}$},
+      title style = {yshift = -2mm},
+      height = 4.5cm,
+      width = 4.5cm,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xlabel = $\vdc~(k_B \tkv/e)$,
+      ylabel = $\vac~(k_B \tkv/e)$,
+      axis on top,
+    ]
+    \addplot graphics[xmin=0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/idc_diff_good.png};
+  \end{axis}
+  \begin{axis}[
+      name = iac_diff,
+      at = (idc_diff.east),
+      anchor = west,
+      xshift = 4mm,
+      title = {$|I_\mathrm{osc}(\Lambda_0=10^9\tkrg) - I_\mathrm{osc}|/I_\mathrm{osc}$},
+      title style = {yshift = -2mm},
+      height = 4.5cm,
+      width = 4.5cm,
+      scale only axis,
+      xmin = -0.480379, xmax = 24.4993,
+      ymin = -0.480379, ymax = 24.4993,
+      xlabel = $\vdc~(k_B \tkv/e)$,
+      yticklabels = {},
+      axis on top,
+    ]
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=0.480379, ymax=24.4993] {../figdata/iac_diff_good.png};
+  \end{axis}
+  %
+  \begin{axis}[
+      name = colorbar,
+      at = (iac.east),
+      anchor = west,
+      ymin = 0,
+      ymax = 1.54775,
+      xmin = 0,
+      xmax = 1,
+      xshift = 5mm,
+      enlargelimits = false,
+      axis on top,
+      width = .75em,
+      height = 4.5cm,
+      scale only axis,
+      tickpos = right,
+      tick align = outside,
+      xtick = \empty,
+      ylabel = {current $(e k_B \tkv/\hbar)$},
+    ]
+    \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=1.54839] {../figdata/colorbar_viridis.png};
+  \end{axis}
+  \begin{axis}[
+      name = colorbar_diff,
+      at = (iac_diff.east),
+      anchor = west,
+      ymin = 0,
+      ymax = 0.000310975,
+      xshift = 5mm,
+      enlargelimits = false,
+      axis on top,
+      width = .75em,
+      height = 4.5cm,
+      scale only axis,
+      tickpos = right,
+      tick align = outside,
+      xtick = \empty,
+      scaled ticks = false,
+      ylabel = {relative deviation},
+    ]
+    \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=0.000310975] {../figdata/colorbar_viridis.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/current_convergence.tex b/tikz/current_convergence_old.tex
similarity index 92%
rename from tikz/current_convergence.tex
rename to tikz/current_convergence_old.tex
index b4acc9e..008405b 100644
--- a/tikz/current_convergence.tex
+++ b/tikz/current_convergence_old.tex
@@ -1,4 +1,4 @@
-\tikzsetnextfilename{current_convergence}%
+\tikzsetnextfilename{current_convergence_old}%
 \begin{tikzpicture}
   \begin{axis}[
       name = idc,
@@ -13,7 +13,7 @@
       ylabel = $\vac~(k_B \tkv/e)$,
       axis on top,
     ]
-    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/idc_converged_good.png};
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/idc_converged_old.png};
   \end{axis}
   \begin{axis}[
       name = iac,
@@ -31,7 +31,7 @@
       yticklabels = {},
       axis on top,
     ]
-    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/iac_converged_good.png};
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/iac_converged_old.png};
   \end{axis}
   \begin{axis}[
       name = idc_diff,
@@ -49,7 +49,7 @@
       ylabel = $\vac~(k_B \tkv/e)$,
       axis on top,
     ]
-    \addplot graphics[xmin=0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/idc_diff_good.png};
+    \addplot graphics[xmin=0.480379, xmax=24.4993, ymin=-0.480379, ymax=24.4993] {../figdata/idc_diff_old.png};
   \end{axis}
   \begin{axis}[
       name = iac_diff,
@@ -67,7 +67,7 @@
       yticklabels = {},
       axis on top,
     ]
-    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=0.480379, ymax=24.4993] {../figdata/iac_diff_good.png};
+    \addplot graphics[xmin=-0.480379, xmax=24.4993, ymin=0.480379, ymax=24.4993] {../figdata/iac_diff_old.png};
   \end{axis}
   %
   \begin{axis}[
-- 
GitLab