From ef05c910e53309d08fb06a9b25182615cdb4e2ba Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Thu, 16 Feb 2023 20:48:22 +0100
Subject: [PATCH] more plots, added TK_VOLTAGE_O3

---
 final_plots.py                      | 115 ++++++--------------
 plot.py                             |  34 +++---
 tikz/contour_wrapper.tex            |  59 -----------
 tikz/convergence-o2-rel.tex         | 158 ++++++++++++++++++++++++++++
 tikz/convergence-o3-rel.tex         | 158 ++++++++++++++++++++++++++++
 tikz/convergence-o3a-rel.tex        | 158 ++++++++++++++++++++++++++++
 tikz/convergence-o3a_ispline.tex    |   2 +-
 tikz/g_overview_vdc_vac_wrapper.tex |  13 ---
 tikz/omega5-3d_wrapper.tex          |  59 -----------
 tikz/overview_wrapper.tex           |  59 -----------
 10 files changed, 527 insertions(+), 288 deletions(-)
 delete mode 100644 tikz/contour_wrapper.tex
 create mode 100644 tikz/convergence-o2-rel.tex
 create mode 100644 tikz/convergence-o3-rel.tex
 create mode 100644 tikz/convergence-o3a-rel.tex
 delete mode 100644 tikz/g_overview_vdc_vac_wrapper.tex
 delete mode 100644 tikz/omega5-3d_wrapper.tex
 delete mode 100644 tikz/overview_wrapper.tex

diff --git a/final_plots.py b/final_plots.py
index bfbf17a..8d87f84 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -28,6 +28,10 @@ from gen_pulse_data import fourier_coef_gauss_symmetric
 from plot_pulses import load_all_pulses_full, integrate_ft
 from kondo import solveTV0_scalar
 
+# color maps with better sampling
+viridis = plt.cm.viridis.resampled(0x10000)
+seismic = plt.cm.seismic.resampled(0x10000)
+
 # In this program all energies are given in units of the RTRG Kondo
 # temperature Tkrg, which is an integration constant of the E-flow RG
 # equations. The more conventional definition of the Kondo temperature is
@@ -36,6 +40,7 @@ from kondo import solveTV0_scalar
 #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_O3 = 3.307473 # for rtol=1e-9, atol=1e-11, voltage_branches=10
 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
@@ -207,74 +212,6 @@ def export_omega5_pgfplots(filename="figdata/omega5_interp.dat", dc_steps=3, ac_
     export_matrix_pgfplots(filename, vdc, vac, gdc, gac, idc, iac, phase, header="vdc vac gdc gac idc iac phase")
 
 
-def export_omega5():
-    "deprecated!"
-    omega = 16.5372
-    dm = DataManager()
-    # Full overview
-    data = dm.list(
-            omega = omega,
-            vdc = None,
-            vac = None,
-            method = "mu",
-            d = 1e9,
-            xL = 0.5,
-            solver_tol_rel = 1e-8,
-            solver_tol_abs = 1e-10,
-            voltage_branches = 4,
-            truncation_order = 3,
-            include_Ga = True,
-            solve_integral_exactly = False,
-            )
-    bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
-            | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
-            | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
-            | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
-            | DataManager.SOLVER_FLAGS["deleted"]
-    data = data.loc[data.solver_flags & bad_flags == 0]
-    vdc, vac, gdc, idc, iac, phase = filter_grid_data(
-            data,
-            omega = omega,
-            vac_min = 0,
-            vac_max = 165.372,
-            vac_num = 51,
-            vdc_min = 0,
-            vdc_max = 165.372,
-            vdc_num = 51,
-            )
-    with open("figdata/vdc_vac_omega16.5372.dat", "w") as file:
-        file.write("vac vdc gdc idc iac")
-        for i in range(51):
-            file.write("\n")
-            np.savetxt(file, np.array([vac[i]/omega, vdc[i]/omega, np.pi*gdc[i], idc[i], iac[i]]).T)
-    ## Lines at constant Vac
-    #for i in range(11):
-    #    with open(f"figdata/vdc_vac{i}_omega16.5372.dat", "w") as file:
-    #        file.write("vac vdc gdc idc iac\n")
-    #        np.savetxt(file, np.array([vac[5*i]/omega, vdc[5*i]/omega, np.pi*gdc[5*i], idc[5*i], iac[5*i]]).T)
-    ## Lines at constant Vdc
-    #for i in range(11):
-    #    with open(f"figdata/vac_vdc{i}_omega16.5372.dat", "w") as file:
-    #        file.write("vac vdc gdc idc iac\n")
-    #        np.savetxt(file, np.array([vac[:,5*i]/omega, vdc[:,5*i]/omega, np.pi*gdc[:,5*i], idc[:,5*i], iac[:,5*i]]).T)
-    # Zoom to primary Kondo peak
-    vdc, vac, gdc, idc, iac, phase = filter_grid_data(
-            data,
-            omega = omega,
-            vac_min = 0,
-            vac_max = 16.5372,
-            vac_num = 21,
-            vdc_min = 0,
-            vdc_max = 16.5372,
-            vdc_num = 21,
-            )
-    with open("figdata/vdc_vac_omega16.5372_zoom.dat", "w") as file:
-        file.write("vac vdc gdc idc iac")
-        for i in range(21):
-            file.write("\n")
-            np.savetxt(file, np.array([vac[i]/omega, vdc[i]/omega, np.pi*gdc[i], idc[i], iac[i]]).T)
-
-
 def export_interp(
         filename,
         fixed_parameter,
@@ -294,6 +231,9 @@ def export_interp(
         o2_scale_x = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
         o2_scale_y = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
         o2_scale_fixed = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
+        o3r_scale_x = TK_VOLTAGE_O3/TK_VOLTAGE,
+        o3r_scale_y = TK_VOLTAGE_O3/TK_VOLTAGE,
+        o3r_scale_fixed = TK_VOLTAGE_O3/TK_VOLTAGE,
         **parameters
         ):
     """
@@ -331,6 +271,7 @@ def export_interp(
     order:
         o2    2nd order RG equations
         o3    3rd order RG equations without Ga and with approximated integral
+        o3r   3rd order RG equations corrected Tk scale
         o3a   3rd order RG equations with approximated integral
         o3p   3rd order RG equations with full integral
 
@@ -365,12 +306,14 @@ def export_interp(
     required_flags = dict(
             o2 = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
             o3 = 0,
+            o3r = 0,
             o3a = DataManager.SOLVER_FLAGS["include_Ga"],
             o3p = DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             )
     bad_flags = dict(
             o2 = DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             o3 = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
+            o3r = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             o3a = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             o3p = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
             )
@@ -379,8 +322,13 @@ def export_interp(
     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])
 
+    o3r_data = dm.list(**{fixed_parameter:fixed_value*o3r_scale_fixed}, **parameters)
+    o3r_data = o3r_data.loc[o3r_data.solver_flags & (global_bad_flags | bad_flags["o3r"] | required_flags["o3r"]) == required_flags["o3r"]]
+    o3r_data = o3r_data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
+    all_data = dict(o2=o2_data, o3r=o3r_data, o3a=data, o3p=data, o3=data)
+
     def selection(order, method, padding=False, **kwargs):
-        thedata = o2_data if order == "o2" else data
+        thedata = all_data[order]
         selected = (thedata.solver_flags & required_flags[order] == required_flags[order]) \
                 & (thedata.solver_flags & bad_flags[order] == 0) \
                 & (thedata.method == method)
@@ -425,6 +373,8 @@ def export_interp(
         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)
+        elif order == "o3r":
+            xy_data = (xy_data[0]/o3r_scale_x, xy_data[1]/o3r_scale_y)
         try:
             results[f"gdc_{method}_{order}{suffix}"] = griddata(
                     xy_data,
@@ -487,14 +437,17 @@ def export_interp(
             **{fixed_parameter:fixed_value, x_parameter:x_mesh, y_parameter:y_mesh},
             **gen_data("o2", "mu"),
             **gen_data("o3", "mu"),
+            **gen_data("o3r", "mu"),
             **gen_data("o3a", "mu"),
             **gen_data("o3p", "mu", integral_method=-1),
             **gen_data("o2", "J", False),
             **gen_data("o3", "J", False),
+            **gen_data("o3r", "J", False),
             **gen_data("o3a", "J", False),
             **gen_data("o3p", "J", False, integral_method=-1),
             **gen_data("o2", "J", True),
             **gen_data("o3", "J", True),
+            **gen_data("o3r", "J", True),
             **gen_data("o3a", "J", True),
             **gen_data("o3p", "J", True, integral_method=-1),
             )
@@ -516,9 +469,6 @@ def export_interp_relative_o3a(
         kyorder = 2,
         spline_s = 2e-4,
         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
         ):
     """
@@ -554,7 +504,6 @@ def export_interp_relative_o3a(
         gamma Γ(E=0) (not an observable)
 
     order:
-        o2    2nd order RG equations
         o3    3rd order RG equations without Ga and with approximated integral
         o3a   3rd order RG equations with approximated integral
         o3p   3rd order RG equations with full integral
@@ -876,9 +825,9 @@ def prepare_RGconvergence():
     dc_indices = np.array([0, 9, 21, 30, 90, 165])
     ac_indices = np.array([32, 52, 80, 124, 172])
 
-    real_cmap = lambda x: plt.cm.seismic((x+1)/2)[...,:3]
-    imwrite(f"figdata/colorbar_seismic.png", np.array(0xffff*plt.cm.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*plt.cm.seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((0x1000,1,3)), format="PNG-FI")
+    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")
     def export_img(name, array):
         imwrite(f"figdata/{name}.png", np.array(0xffff*real_cmap(array[::-1]), dtype=np.uint16), format="PNG-FI")
 
@@ -888,9 +837,9 @@ def prepare_RGconvergence():
 
     file = open("tikz/colorscales.tex", "w")
 
-    # O2
     omega = 16.5372
-    omega_o2 = 48.695054
+    #omega_o2 = 48.695054
+    #omega_o3r = 15.8884
     vdc = data["vdc"][0] / omega
     vac = data["vac"][:,0] / omega
 
@@ -900,7 +849,7 @@ def prepare_RGconvergence():
 \\newcommand\\extentymin{{{1.5*vac[0]-0.5*vac[1]}}}
 \\newcommand\\extentymax{{{1.5*vac[-1]-0.5*vac[-2]}}}""", file=file)
 
-    for method, u, ref in (("o2","mu","o3p"), ("o3","mu","o3p"), ("o3a","mu","o3p"), ("o3a_ispline","mu","o3a"), ("o3a","J","o3a"), ("o3a_p","J","o3a")):
+    for method, u, ref in (("o2","mu","o3p"), ("o3","mu","o3p"), ("o3a","mu","o3p"), ("o3a_ispline","mu","o3a"), ("o3a","J","o3a"), ("o3a_p","J","o3a"), ("o3r","mu","o3p")):
         gdiff = data[f"gdc_{u}_{method}"] - data[f"gdc_mu_{ref}"]
         if u != "mu":
             method += "_" + u
@@ -915,7 +864,7 @@ def prepare_RGconvergence():
         print(f"Color scale for {method} relative: {normalize}")
         print(f"\\newcommand\\scale{tex_name}rel{{{normalize}}}% convergence_{method}_relative.png", file=file)
 
-    for method, u, sign in (("o3p","mu",-1), ("o3a_vb7","mu",-1), ("o3a","J",1), ("o3a_p","J",1)):
+    for method, u, sign in (("o3p","mu",-1), ("o3a_vb7","mu",-1), ("o3a","J",1), ("o3a_p","J",1),("o3","mu",1),("o3","J",1),("o3_p","J",1)):
         gdiff = sign*data_rel[f"gdc_{u}_{method}"]
         if u != "mu":
             method += "_" + u
@@ -1752,7 +1701,7 @@ def export_convergence_pgfplots(simplified_initial_conditions = False):
 
     suffix = "_bad" if simplified_initial_conditions else "_good"
     def export_img(name, array):
-        imwrite(f"figdata/{name}{suffix}.png", np.array(0xffff*plt.cm.viridis(array[::-1])[...,:3], dtype=np.uint16), format="PNG-FI")
+        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)
     iac_max = np.nanmax(iac_fit_a)
     i_max = max(idc_max, 2*iac_max)
@@ -1805,7 +1754,7 @@ def export_floquet_matrices_pgfplots():
 
     scale = lambda x: np.log10(np.abs(x))/9 + 1
     def export_img(name, array):
-        imwrite(f"figdata/{name}.png", np.array(0xffff*plt.cm.viridis(scale(array))[...,:3], dtype=np.uint16), format="PNG-FI")
+        imwrite(f"figdata/{name}.png", np.array(0xffff*viridis(scale(array))[...,:3], dtype=np.uint16), format="PNG-FI")
 
     z_ref_j = z_matrices[("J", True, True)]
     z_ref_mu = z_matrices[("mu", True, False)]
@@ -1921,7 +1870,7 @@ def export_Ga_max_pgfplots():
     grid_g2_max = griddata((vdc, vac), g2_max, (vdc_mesh, vac_mesh), method="cubic")
 
     def export_img(name, array):
-        imwrite(f"figdata/{name}.png", np.array(0xffff*plt.cm.viridis(array[::-1])[...,:3], dtype=np.uint16), format="PNG-FI")
+        imwrite(f"figdata/{name}.png", np.array(0xffff*viridis(array[::-1])[...,:3], dtype=np.uint16), format="PNG-FI")
 
     print(f"Extent: xmin={extent[0]/omega:.6g}, xmax={extent[1]/omega:,.6g}, ymin={extent[2]/omega:.6g}, ymax={extent[3]/omega:,.6g}")
     norm = np.nanmax(np.abs(grid_ga_max))
diff --git a/plot.py b/plot.py
index 10a0447..96a0c45 100644
--- a/plot.py
+++ b/plot.py
@@ -44,6 +44,7 @@ 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_O3 = 3.44428 # for rtol=1e-8, atol=1e-10, voltage_branches=4
 TK_VOLTAGE_ORDER2 = 10.1368086 # for rtol=1e-9, atol=1e-11, voltage_branches=10
 #TK_VOLTAGE_ORDER2 = 10.13754 # with correction
 
@@ -101,7 +102,7 @@ def main():
             help="Truncation order of RG equations")
     parser.add_argument("--integral_method", metavar="int", type=int, default=-15,
             help="method for solving frequency integral (-1 for exact solution)")
-    parser.add_argument("--include_Ga", metavar="bool", type=bool, default=False,
+    parser.add_argument("--include_Ga", metavar="bool", type=int, default=1,
             help="include Ga in RG equations")
     parser.add_argument("--good_flags", metavar="int", type=int, default=0,
             help = "Require these flags")
@@ -1366,8 +1367,10 @@ def convergence_overview(dm, **parameters):
     vdc_max = 82.686
     vac_max = 82.686
     exponent = 3
-    if parameters["good_flags"] & 8:
+    if parameters["good_flags"] & DataManager.SOLVER_FLAGS["simplified_initial_conditions"]:
         exponent = 2
+    elif parameters["good_flags"] & DataManager.SOLVER_FLAGS["improved_initial_conditions"]:
+        exponent = 4
     d_num = 5
     log10d_min = 5
     log10d_max = 9
@@ -1403,8 +1406,8 @@ def convergence_overview(dm, **parameters):
     ax3.plot(logd_inv3, -gdc_fit_b.reshape((1,-1))*logd_inv3.reshape((d_num,1)), color="gray", linewidth=0.5)
     ax3.plot(logd_inv3, gdc_fit_a.reshape((1,-1)) - gdc.reshape((d_num,-1)))
 
-    #fig, ((ax1, ax2, ax3), (ax4, ax5, ax6), (ax7, ax8, ax9)) = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True)
-    fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=True)
+    fig, ((ax1, ax2, ax3), (ax4, ax5, ax6), (ax7, ax8, ax9)) = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True)
+    #fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=True)
     extent = (-0.5*vdc_max/(vdc_num-1), vdc_max*(1+0.5/(vdc_num-1)), -0.5*vac_max/(vac_num-1), vac_max*(1+0.5/(vac_num-1)))
     ax1.set_title(r"$I_\mathrm{avg}(\Lambda_0=\infty)$")
     ax2.set_title(r"$\frac{|I_\mathrm{avg}(\Lambda_0=10^9 T_K') - I_\mathrm{avg}(\Lambda_0=\infty)|}{I_\mathrm{avg}(\Lambda_0=10^9T_K')}$")
@@ -1427,12 +1430,12 @@ def convergence_overview(dm, **parameters):
     fig.colorbar(img, ax=ax5)
     img = ax6.imshow(-iac_diff, extent=extent, origin="lower")
     fig.colorbar(img, ax=ax6)
-    #img = ax7.imshow(gdc_fit_a, extent=extent, origin="lower")
-    #fig.colorbar(img, ax=ax7)
-    #img = ax8.imshow(gdc_fit_b*logd_inv3[4]/gdc_fit_a, extent=extent, origin="lower")
-    #fig.colorbar(img, ax=ax8)
-    #img = ax9.imshow(gdc_diff, extent=extent, origin="lower")
-    #fig.colorbar(img, ax=ax9)
+    img = ax7.imshow(gdc_fit_a, extent=extent, origin="lower")
+    fig.colorbar(img, ax=ax7)
+    img = ax8.imshow(gdc_fit_b*logd_inv3[4]/gdc_fit_a, extent=extent, origin="lower")
+    fig.colorbar(img, ax=ax8)
+    img = ax9.imshow(gdc_diff, extent=extent, origin="lower")
+    fig.colorbar(img, ax=ax9)
 
     #cmap = plt.cm.seismic
     #gdc_diff = (merged.dc_conductance_1 - merged.dc_conductance_2)/merged.dc_conductance_2
@@ -1833,8 +1836,8 @@ def RGeq_comparison_contours2(
         **trashparams):
     omega = 16.5372
     kwargs = dict(linewidths=0.8)
-    colors = dict(o2="#ff0000", o3="#6000ff", o3a="#00f000", o3p="#000000", ref="#808080")
-    zorder = dict(o2=1, o3=2, o3a=3, o3p=4)
+    colors = dict(o2="#ff0000", o3="#6000ff", o3r="#ff00ff", o3a="#00f000", o3p="#000000", ref="#808080")
+    zorder = dict(o2=1, o3=2, o3r=5, o3a=3, o3p=4)
     approximations = list(zorder.keys())
     data = np.load(filename)
     vdc = data["vdc"] / omega
@@ -1842,8 +1845,11 @@ def RGeq_comparison_contours2(
     levels = (1/np.linspace(1, 12, 12)[::-1]) / prefactor
     fig, ax = plt.subplots(nrows=1, ncols=1, sharex=True, sharey=True)
     for approx in approximations:
-        contour = ax.contour(vdc, vac, data[f"{observable}_mu_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
-        save_contour_coordinates(contour, f"figdata/contour_{observable}_mu_{approx}.dat", "vdc", "vac", "gdc")
+        try:
+            contour = ax.contour(vdc, vac, data[f"{observable}_mu_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+            save_contour_coordinates(contour, f"figdata/contour_{observable}_mu_{approx}.dat", "vdc", "vac", "gdc")
+        except KeyError:
+            settings.logger.warning(f"No data available for {observable}_mu_{approx}")
         try:
             contour = ax.contour(vdc, vac, data[f"{observable}_J_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
             save_contour_coordinates(contour, f"figdata/contour_{observable}_J_{approx}.dat", "vdc", "vac", "gdc")
diff --git a/tikz/contour_wrapper.tex b/tikz/contour_wrapper.tex
deleted file mode 100644
index fd51447..0000000
--- a/tikz/contour_wrapper.tex
+++ /dev/null
@@ -1,59 +0,0 @@
-\documentclass{standalone}
-\newcommand\tikzsetnextfilename[1]{}
-\newcommand\tkv{T_\mathrm{K}}
-\newcommand\vdc{V_\mathrm{avg}}
-\newcommand\vac{V_\mathrm{osc}}
-\newcommand\gdc{G_\mathrm{avg}}
-\newcommand\gac{G_\mathrm{osc}}
-\usepackage{pgfplots}
-\pgfplotsset{compat=1.18}
-\pgfplotsset{surf shading/precision=pdf}
-\usepackage{fontspec}
-\defaultfontfeatures{Ligatures=TeX}
-\setsansfont[
-  UprightFont = *-Regular,
-  BoldFont = *-Bold,
-  ItalicFont = *-Italic,
-  %BoldItalicFont = LinBiolinum_K, % TODO: Check this font, use a better one
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSans}
-\setmainfont[
-  %UprightFont = *Display-Regular,
-  UprightFont = *-Regular,
-  BoldFont = *-Semibold,
-  ItalicFont = *-Italic,
-  BoldItalicFont = *-SemiboldItalic,
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  %SlantedFont = *-Italic,
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSerif}
-%\setmonofont[
-%  UprightFont = *-Regular
-%]{LibertinusMono}
-\setmonofont[
-  UprightFont = *-Regular,
-  BoldFont = *-Medium,
-  ItalicFont = *-It,
-  BoldItalicFont = *-MediumIt,
-]{SourceCodePro}
-
-\usepackage{mathtools}
-%\usepackage{amsmath}
-\usepackage{amssymb}
-%\usepackage{amsfonts}
-\usepackage[math-style=ISO, bold-style=TeX, partial=upright, nabla=upright]{unicode-math}
-\setmathfont{Libertinus Math}
-\setmathfont{STIXTwoMath-Regular.otf}[range=cal]
-\setmathfont{STIXTwoMath-Regular.otf}[range={\mathscr},StylisticSet=1]
-\setmathfont{STIXTwoMath-Regular.otf}[range={"02022, "02218, "02200, "0223C, "02248}] % Slightly larger \bullet and \circ; larger \forall; different \sim, \approx
-%\setmathfont{Latin Modern Math}[range={"027F8-"027FA}] % \iff, \implies, \impliedby
-%\setmathfont{STIX Two Math}[range={"027E8, "027E9}] % \langle, \rangle
-\setmathfont{STIXTwoMath-Regular.otf}[range={"0221A-"0221C, "023B7}] % root (sqrt, cbrt, \dots) symbols
-\begin{document}
-\input{contour}
-\end{document}
diff --git a/tikz/convergence-o2-rel.tex b/tikz/convergence-o2-rel.tex
new file mode 100644
index 0000000..0cee967
--- /dev/null
+++ b/tikz/convergence-o2-rel.tex
@@ -0,0 +1,158 @@
+\tikzsetnextfilename{convergence_o2}%
+\begin{tikzpicture}
+  \pgfplotsset{set layers = axis on top}
+  \begin{axis}[
+      name = img,
+      anchor = north west,
+      width = 4cm,
+      height = 4cm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0, ymax = 10,
+      xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\vac~(\Omega)$},
+      every axis plot post/.append style = {
+        line width=0.2pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \addplot graphics[xmin=\extentxmin, xmax=\extentxmax, ymin=\extentymin, ymax=\extentymax] {../figdata/convergence_o2_relative.png};
+    \node at (2,9.2) {(b)};
+    % vertical lines
+    \addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
+    \addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
+    \addplot[no marks, black] coordinates {(0,4.0) (10,4.0)};
+    \addplot[no marks, black] coordinates {(0,6.2) (10,6.2)};
+    \addplot[no marks, black] coordinates {(0,8.6) (10,8.6)};
+    \node[anchor=west] at (7.0,1.6) {\small △};
+    \node[anchor=west] at (7.5,2.6) {\small ▽};
+    \node[anchor=west] at (8.0,4.0) {\small □};
+    \node[anchor=west] at (8.5,6.2) {\small ○};
+    \node[anchor=west] at (9.0,8.6) {\small ◊};
+    % horizontal lines
+    \addplot[no marks, black] coordinates {(0,0) (0,10)};
+    \addplot[no marks, black] coordinates {(0.3,0) (0.3,10)};
+    \addplot[no marks, black] coordinates {(0.7,0) (0.7,10)};
+    \addplot[no marks, black] coordinates {(1,0) (1,10)};
+    \addplot[no marks, black] coordinates {(3,0) (3,10)};
+    \addplot[no marks, black] coordinates {(5.5,0) (5.5,10)};
+    \node[anchor=south] at (0,8.5) {\small ▶};
+    \node[anchor=south] at (0.3,7) {\small ◀};
+    \node[anchor=south] at (0.7,5) {\small ■};
+    \node[anchor=south] at (1,4.0) {\large •};
+    \node[anchor=south] at (3,0.2) {\small ♦};
+    \node[anchor=south] at (5.5,0.2) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      at = (img.north),
+      name = ac,
+      anchor = south,
+      width = 4cm,
+      height = 3cm,
+      yshift = 2.5mm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0.08, ymax = 0.88,
+      %xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\gdc~(2e^2/h)$},
+      %xticklabel pos = top,
+      xticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+      legend entries = {leading order, full RG equations},
+      legend cell align = left,
+      legend style = {
+        legend pos = outer north east,
+        draw = none,
+        fill = none,
+      },
+    ]
+    \node at (2,0.79) {(a)};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3p_ac.dat};
+    \node[anchor=west] at (7.0,0.61) {\small △};
+    \node[anchor=west] at (7.5,0.48) {\small ▽};
+    \node[anchor=west] at (8.0,0.365) {\small □};
+    \node[anchor=west] at (8.5,0.251) {\small ○};
+    \node[anchor=west] at (9.0,0.145) {\small ◊};
+  \end{axis}
+  \begin{axis}[
+      name = dc,
+      at = (img.east),
+      anchor = west,
+      width = 3cm,
+      height = 4cm,
+      xshift = 2.5mm,
+      scale only axis,
+      xmin = 0.05, xmax = 1.62,
+      ymin = 0, ymax = 10,
+      xlabel = {$\gdc~(2e^2/h)$},
+      %ylabel = {$\vac~(\Omega)$},
+      %yticklabel pos = right,
+      yticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \node at (1.3,9.2) {(c)};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \node[anchor=south] at (0.72,10) {\small ▶};
+    \node[anchor=south] at (0.59,10) {\small ◀};
+    \node[anchor=south] at (0.47,10) {\small ■};
+    \node[anchor=south] at (0.35,10) {\large •};
+    \node[anchor=south] at (0.23,10) {\small ♦};
+    \node[anchor=south] at (0.11,10) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      name = colorbar,
+      at = (dc.north west),
+      xmin = -\scaleolrel,
+      xmax = \scaleolrel,
+      xshift = 5mm,
+      yshift = 11mm,
+      anchor = south west,
+      enlargelimits = false,
+      axis on top,
+      height = .75em,
+      width = 2.5cm,
+      scale only axis,
+      tickpos = bottom,
+      tick align = outside,
+      ytick = \empty,
+      xtick = {-0.05, 0, 0.05},
+      xticklabels = {$-5\%$, $0$, $5\%$},
+      scaled ticks = false,
+      xlabel = {relative deviation},
+      xlabel shift = -14.6mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=-\scaleolrel, xmax=\scaleolrel] {../figdata/colorbar_seismic.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/convergence-o3-rel.tex b/tikz/convergence-o3-rel.tex
new file mode 100644
index 0000000..c8f5fdc
--- /dev/null
+++ b/tikz/convergence-o3-rel.tex
@@ -0,0 +1,158 @@
+\tikzsetnextfilename{convergence_o3}%
+\begin{tikzpicture}
+  \pgfplotsset{set layers = axis on top}
+  \begin{axis}[
+      name = img,
+      anchor = north west,
+      width = 4cm,
+      height = 4cm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0, ymax = 10,
+      xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\vac~(\Omega)$},
+      every axis plot post/.append style = {
+        line width=0.2pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \addplot graphics[xmin=\extentxmin, xmax=\extentxmax, ymin=\extentymin, ymax=\extentymax] {../figdata/convergence_o3_relative.png};
+    \node at (2,9.2) {(b)};
+    % vertical lines
+    \addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
+    \addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
+    \addplot[no marks, black] coordinates {(0,4.0) (10,4.0)};
+    \addplot[no marks, black] coordinates {(0,6.2) (10,6.2)};
+    \addplot[no marks, black] coordinates {(0,8.6) (10,8.6)};
+    \node[anchor=west] at (7.0,1.6) {\small △};
+    \node[anchor=west] at (7.5,2.6) {\small ▽};
+    \node[anchor=west] at (8.0,4.0) {\small □};
+    \node[anchor=west] at (8.5,6.2) {\small ○};
+    \node[anchor=west] at (9.0,8.6) {\small ◊};
+    % horizontal lines
+    \addplot[no marks, black] coordinates {(0,0) (0,10)};
+    \addplot[no marks, black] coordinates {(0.3,0) (0.3,10)};
+    \addplot[no marks, black] coordinates {(0.7,0) (0.7,10)};
+    \addplot[no marks, black] coordinates {(1,0) (1,10)};
+    \addplot[no marks, black] coordinates {(3,0) (3,10)};
+    \addplot[no marks, black] coordinates {(5.5,0) (5.5,10)};
+    \node[anchor=south] at (0,8.5) {\small ▶};
+    \node[anchor=south] at (0.3,7) {\small ◀};
+    \node[anchor=south] at (0.7,5) {\small ■};
+    \node[anchor=south] at (1,4.0) {\large •};
+    \node[anchor=south] at (3,0.2) {\small ♦};
+    \node[anchor=south] at (5.5,0.2) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      at = (img.north),
+      name = ac,
+      anchor = south,
+      width = 4cm,
+      height = 3cm,
+      yshift = 2.5mm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0.08, ymax = 0.88,
+      %xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\gdc~(2e^2/h)$},
+      %xticklabel pos = top,
+      xticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+      legend entries = {without $G^a$, full RG equations},
+      legend cell align = left,
+      legend style = {
+        legend pos = outer north east,
+        draw = none,
+        fill = none,
+      },
+    ]
+    \node at (2,0.79) {(a)};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3p_ac.dat};
+    \node[anchor=west] at (7.0,0.61) {\small △};
+    \node[anchor=west] at (7.5,0.48) {\small ▽};
+    \node[anchor=west] at (8.0,0.365) {\small □};
+    \node[anchor=west] at (8.5,0.251) {\small ○};
+    \node[anchor=west] at (9.0,0.145) {\small ◊};
+  \end{axis}
+  \begin{axis}[
+      name = dc,
+      at = (img.east),
+      anchor = west,
+      width = 3cm,
+      height = 4cm,
+      xshift = 2.5mm,
+      scale only axis,
+      xmin = 0.05, xmax = 1.62,
+      ymin = 0, ymax = 10,
+      xlabel = {$\gdc~(2e^2/h)$},
+      %ylabel = {$\vac~(\Omega)$},
+      %yticklabel pos = right,
+      yticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \node at (1.3,9.2) {(c)};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \node[anchor=south] at (0.72,10) {\small ▶};
+    \node[anchor=south] at (0.59,10) {\small ◀};
+    \node[anchor=south] at (0.47,10) {\small ■};
+    \node[anchor=south] at (0.35,10) {\large •};
+    \node[anchor=south] at (0.23,10) {\small ♦};
+    \node[anchor=south] at (0.11,10) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      name = colorbar,
+      at = (dc.north west),
+      xmin = -\scaleoxrel,
+      xmax = \scaleoxrel,
+      xshift = 5mm,
+      yshift = 11mm,
+      anchor = south west,
+      enlargelimits = false,
+      axis on top,
+      height = .75em,
+      width = 2.5cm,
+      scale only axis,
+      tickpos = bottom,
+      tick align = outside,
+      ytick = \empty,
+      xtick = {-0.05, 0, 0.05},
+      xticklabels = {$-5\%$, $0$, $5\%$},
+      scaled ticks = false,
+      xlabel = {relative deviation},
+      xlabel shift = -14.6mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=-\scaleoxrel, xmax=\scaleoxrel] {../figdata/colorbar_seismic.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/convergence-o3a-rel.tex b/tikz/convergence-o3a-rel.tex
new file mode 100644
index 0000000..f2217d4
--- /dev/null
+++ b/tikz/convergence-o3a-rel.tex
@@ -0,0 +1,158 @@
+\tikzsetnextfilename{convergence_o3a}%
+\begin{tikzpicture}
+  \pgfplotsset{set layers = axis on top}
+  \begin{axis}[
+      name = img,
+      anchor = north west,
+      width = 4cm,
+      height = 4cm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0, ymax = 10,
+      xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\vac~(\Omega)$},
+      every axis plot post/.append style = {
+        line width=0.2pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \addplot graphics[xmin=\extentxmin, xmax=\extentxmax, ymin=\extentymin, ymax=\extentymax] {../figdata/convergence_o3p_relative_fromdiff.png};
+    \node at (2,9.2) {(b)};
+    % vertical lines
+    \addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
+    \addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
+    \addplot[no marks, black] coordinates {(0,4.0) (10,4.0)};
+    \addplot[no marks, black] coordinates {(0,6.2) (10,6.2)};
+    \addplot[no marks, black] coordinates {(0,8.6) (10,8.6)};
+    \node[anchor=west] at (7.0,1.6) {\small △};
+    \node[anchor=west] at (7.5,2.6) {\small ▽};
+    \node[anchor=west] at (8.0,4.0) {\small □};
+    \node[anchor=west] at (8.5,6.2) {\small ○};
+    \node[anchor=west] at (9.0,8.6) {\small ◊};
+    % horizontal lines
+    \addplot[no marks, black] coordinates {(0,0) (0,10)};
+    \addplot[no marks, black] coordinates {(0.3,0) (0.3,10)};
+    \addplot[no marks, black] coordinates {(0.7,0) (0.7,10)};
+    \addplot[no marks, black] coordinates {(1,0) (1,10)};
+    \addplot[no marks, black] coordinates {(3,0) (3,10)};
+    \addplot[no marks, black] coordinates {(5.5,0) (5.5,10)};
+    \node[anchor=south] at (0,8.5) {\small ▶};
+    \node[anchor=south] at (0.3,7) {\small ◀};
+    \node[anchor=south] at (0.7,5) {\small ■};
+    \node[anchor=south] at (1,4.0) {\large •};
+    \node[anchor=south] at (3,0.2) {\small ♦};
+    \node[anchor=south] at (5.5,0.2) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      at = (img.north),
+      name = ac,
+      anchor = south,
+      width = 4cm,
+      height = 3cm,
+      yshift = 2.5mm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0.08, ymax = 0.88,
+      %xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\gdc~(2e^2/h)$},
+      %xticklabel pos = top,
+      xticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+      legend entries = {approx.~integral, full RG equations},
+      legend cell align = left,
+      legend style = {
+        legend pos = outer north east,
+        draw = none,
+        fill = none,
+      },
+    ]
+    \node at (2,0.79) {(a)};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3p_ac.dat};
+    \node[anchor=west] at (7.0,0.61) {\small △};
+    \node[anchor=west] at (7.5,0.48) {\small ▽};
+    \node[anchor=west] at (8.0,0.365) {\small □};
+    \node[anchor=west] at (8.5,0.251) {\small ○};
+    \node[anchor=west] at (9.0,0.145) {\small ◊};
+  \end{axis}
+  \begin{axis}[
+      name = dc,
+      at = (img.east),
+      anchor = west,
+      width = 3cm,
+      height = 4cm,
+      xshift = 2.5mm,
+      scale only axis,
+      xmin = 0.05, xmax = 1.62,
+      ymin = 0, ymax = 10,
+      xlabel = {$\gdc~(2e^2/h)$},
+      %ylabel = {$\vac~(\Omega)$},
+      %yticklabel pos = right,
+      yticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \node at (1.3,9.2) {(c)};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \node[anchor=south] at (0.72,10) {\small ▶};
+    \node[anchor=south] at (0.59,10) {\small ◀};
+    \node[anchor=south] at (0.47,10) {\small ■};
+    \node[anchor=south] at (0.35,10) {\large •};
+    \node[anchor=south] at (0.23,10) {\small ♦};
+    \node[anchor=south] at (0.11,10) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      name = colorbar,
+      at = (dc.north west),
+      xmin = -\scaleoxprelX,
+      xmax = \scaleoxprelX,
+      xshift = 5mm,
+      yshift = 11mm,
+      anchor = south west,
+      enlargelimits = false,
+      axis on top,
+      height = .75em,
+      width = 2.5cm,
+      scale only axis,
+      tickpos = bottom,
+      tick align = outside,
+      ytick = \empty,
+      xtick = {-0.005, 0, 0.005},
+      xticklabels = {$-0.5\%$, $0$, $0.5\%$},
+      scaled ticks = false,
+      xlabel = {relative deviation},
+      xlabel shift = -14.6mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=-\scaleoxprelX, xmax=\scaleoxprelX] {../figdata/colorbar_seismic.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/convergence-o3a_ispline.tex b/tikz/convergence-o3a_ispline.tex
index 53542c1..93e34fc 100644
--- a/tikz/convergence-o3a_ispline.tex
+++ b/tikz/convergence-o3a_ispline.tex
@@ -153,7 +153,7 @@
       xticklabels = {$-5\%$, $0$, $5\%$},
       scaled ticks = false,
       xlabel = {relative deviation},
-      xlabel shift = -14mm,
+      xlabel shift = -14.6mm,
     ]
     \addplot graphics[ymin=0, ymax=1, xmin=-\scaleoxaxisplinerel, xmax=\scaleoxaxisplinerel] {../figdata/colorbar_seismic.png};
   \end{axis}
diff --git a/tikz/g_overview_vdc_vac_wrapper.tex b/tikz/g_overview_vdc_vac_wrapper.tex
deleted file mode 100644
index c8580d4..0000000
--- a/tikz/g_overview_vdc_vac_wrapper.tex
+++ /dev/null
@@ -1,13 +0,0 @@
-\documentclass{standalone}
-\newcommand\tikzsetnextfilename[1]{}
-\newcommand\tkv{T_\mathrm{K}}
-\newcommand\vdc{V_\mathrm{dc}}
-\newcommand\vac{V_\mathrm{ac}}
-\newcommand\gdc{G_\mathrm{dc}}
-\newcommand\gac{G_\mathrm{ac}}
-\usepackage{pgfplots}
-\pgfplotsset{compat=1.18}
-\pgfplotsset{surf shading/precision=pdf}
-\begin{document}
-\input{g_overview_vdc_vac}
-\end{document}
diff --git a/tikz/omega5-3d_wrapper.tex b/tikz/omega5-3d_wrapper.tex
deleted file mode 100644
index 6eeb7c4..0000000
--- a/tikz/omega5-3d_wrapper.tex
+++ /dev/null
@@ -1,59 +0,0 @@
-\documentclass{standalone}
-\newcommand\tikzsetnextfilename[1]{}
-\newcommand\tkv{T_\mathrm{K}}
-\newcommand\vdc{V_\mathrm{avg}}
-\newcommand\vac{V_\mathrm{osc}}
-\newcommand\gdc{G_\mathrm{avg}}
-\newcommand\gac{G_\mathrm{osc}}
-\usepackage{pgfplots}
-\pgfplotsset{compat=1.18}
-\pgfplotsset{surf shading/precision=pdf}
-\usepackage{fontspec}
-\defaultfontfeatures{Ligatures=TeX}
-\setsansfont[
-  UprightFont = *-Regular,
-  BoldFont = *-Bold,
-  ItalicFont = *-Italic,
-  %BoldItalicFont = LinBiolinum_K, % TODO: Check this font, use a better one
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSans}
-\setmainfont[
-  %UprightFont = *Display-Regular,
-  UprightFont = *-Regular,
-  BoldFont = *-Semibold,
-  ItalicFont = *-Italic,
-  BoldItalicFont = *-SemiboldItalic,
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  %SlantedFont = *-Italic,
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSerif}
-%\setmonofont[
-%  UprightFont = *-Regular
-%]{LibertinusMono}
-\setmonofont[
-  UprightFont = *-Regular,
-  BoldFont = *-Medium,
-  ItalicFont = *-It,
-  BoldItalicFont = *-MediumIt,
-]{SourceCodePro}
-
-\usepackage{mathtools}
-%\usepackage{amsmath}
-\usepackage{amssymb}
-%\usepackage{amsfonts}
-\usepackage[math-style=ISO, bold-style=TeX, partial=upright, nabla=upright]{unicode-math}
-\setmathfont{Libertinus Math}
-\setmathfont{STIXTwoMath-Regular.otf}[range=cal]
-\setmathfont{STIXTwoMath-Regular.otf}[range={\mathscr},StylisticSet=1]
-\setmathfont{STIXTwoMath-Regular.otf}[range={"02022, "02218, "02200, "0223C, "02248}] % Slightly larger \bullet and \circ; larger \forall; different \sim, \approx
-%\setmathfont{Latin Modern Math}[range={"027F8-"027FA}] % \iff, \implies, \impliedby
-%\setmathfont{STIX Two Math}[range={"027E8, "027E9}] % \langle, \rangle
-\setmathfont{STIXTwoMath-Regular.otf}[range={"0221A-"0221C, "023B7}] % root (sqrt, cbrt, \dots) symbols
-\begin{document}
-\input{omega5-3d}
-\end{document}
diff --git a/tikz/overview_wrapper.tex b/tikz/overview_wrapper.tex
deleted file mode 100644
index 4efc8fc..0000000
--- a/tikz/overview_wrapper.tex
+++ /dev/null
@@ -1,59 +0,0 @@
-\documentclass{standalone}
-\newcommand\tikzsetnextfilename[1]{}
-\newcommand\tkv{T_\mathrm{K}}
-\newcommand\vdc{V_\mathrm{avg}}
-\newcommand\vac{V_\mathrm{osc}}
-\newcommand\gdc{G_\mathrm{avg}}
-\newcommand\gac{G_\mathrm{osc}}
-\usepackage{pgfplots}
-\pgfplotsset{compat=1.18}
-\pgfplotsset{surf shading/precision=pdf}
-\usepackage{fontspec}
-\defaultfontfeatures{Ligatures=TeX}
-\setsansfont[
-  UprightFont = *-Regular,
-  BoldFont = *-Bold,
-  ItalicFont = *-Italic,
-  %BoldItalicFont = LinBiolinum_K, % TODO: Check this font, use a better one
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSans}
-\setmainfont[
-  %UprightFont = *Display-Regular,
-  UprightFont = *-Regular,
-  BoldFont = *-Semibold,
-  ItalicFont = *-Italic,
-  BoldItalicFont = *-SemiboldItalic,
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  %SlantedFont = *-Italic,
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSerif}
-%\setmonofont[
-%  UprightFont = *-Regular
-%]{LibertinusMono}
-\setmonofont[
-  UprightFont = *-Regular,
-  BoldFont = *-Medium,
-  ItalicFont = *-It,
-  BoldItalicFont = *-MediumIt,
-]{SourceCodePro}
-
-\usepackage{mathtools}
-%\usepackage{amsmath}
-\usepackage{amssymb}
-%\usepackage{amsfonts}
-\usepackage[math-style=ISO, bold-style=TeX, partial=upright, nabla=upright]{unicode-math}
-\setmathfont{Libertinus Math}
-\setmathfont{STIXTwoMath-Regular.otf}[range=cal]
-\setmathfont{STIXTwoMath-Regular.otf}[range={\mathscr},StylisticSet=1]
-\setmathfont{STIXTwoMath-Regular.otf}[range={"02022, "02218, "02200, "0223C, "02248}] % Slightly larger \bullet and \circ; larger \forall; different \sim, \approx
-%\setmathfont{Latin Modern Math}[range={"027F8-"027FA}] % \iff, \implies, \impliedby
-%\setmathfont{STIX Two Math}[range={"027E8, "027E9}] % \langle, \rangle
-\setmathfont{STIXTwoMath-Regular.otf}[range={"0221A-"0221C, "023B7}] % root (sqrt, cbrt, \dots) symbols
-\begin{document}
-\input{overview}
-\end{document}
-- 
GitLab