From 1200a2822769e2c78c89d0260a0954d57606a6eb Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Sat, 31 Dec 2022 00:44:55 +0100
Subject: [PATCH] more tikz plots

---
 final_plots.py           | 268 ++++++++++++++++++++++++++++++++++++---
 plot.py                  |  31 +++--
 test/benchmark.py        |   7 +-
 tikz/commutator_RG2.tex  |  39 ++++++
 tikz/contour.tex         |  24 ++--
 tikz/convergence-o2.tex  | 164 ++++++++++++++++++++++++
 tikz/convergence-o3.tex  | 164 ++++++++++++++++++++++++
 tikz/convergence-o3a.tex | 163 ++++++++++++++++++++++++
 tikz/omega5-3d.tex       |  63 ++++++---
 tikz/overview.tex        |   6 +-
 10 files changed, 866 insertions(+), 63 deletions(-)
 create mode 100644 tikz/commutator_RG2.tex
 create mode 100644 tikz/convergence-o2.tex
 create mode 100644 tikz/convergence-o3.tex
 create mode 100644 tikz/convergence-o3a.tex

diff --git a/final_plots.py b/final_plots.py
index 8b05af2..e1d0c58 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -16,6 +16,7 @@ import numpy as np
 from scipy.interpolate import bisplrep, bisplev, splrep, BSpline, griddata, interp1d
 from scipy.special import jn
 from imageio import imwrite
+from scipy.integrate import quad, quad_vec
 import settings
 from data_management import DataManager, KondoImport
 from kondo import gen_init_matrix
@@ -101,13 +102,14 @@ def interp_vac0(dm, xL=0.5, include_Ga=True, truncation_order=3, integral_method
     return interp1d(vdc, gdc, kind="cubic")
 
 
-def photon_assisted_tunneling(dm, omega, vdc, vac, nmax=50, **parameters):
+def photon_assisted_tunneling(dm=None, omega=0, vdc=0, vac=0, nmax=50, interp=None, **parameters):
     """
     Return differential conductance G calculated from photon
     assisted tunneling at given parameters.
     """
     result = np.zeros(np.broadcast_shapes(np.shape(omega), np.shape(vdc), np.shape(vac)), dtype=np.float64)
-    interp = interp_vac0(dm, **parameters)
+    if interp is None:
+        interp = interp_vac0(dm, **parameters)
     for n in range(-nmax, nmax+1):
         try:
             result += jn(n, vac/omega)**2 * interp(vdc+n*omega)
@@ -176,7 +178,31 @@ def filter_grid_data(
     return *np.meshgrid(vdc_arr, vac_arr), gdc_arr, idc_arr, iac_arr, phase_arr
 
 
+def export_matrix_pgfplots(filename, *arrays, header="", fmt="%.6g"):
+    array = np.array(arrays).T
+    with open(filename, "w") as file:
+        if header:
+            file.write(header)
+        for sector in array:
+            file.write("\n")
+            np.savetxt(file, sector, fmt=fmt)
+
+
+def export_omega5_pgfplots(filename="figdata/omega5_interp.dat", dc_steps=3, ac_steps=2):
+    data = np.load("figdata/omega5_interp.npz")
+    omega = 16.5372
+    vdc = data["vdc"][::ac_steps,::dc_steps]/omega
+    vac = data["vac"][::ac_steps,::dc_steps]/omega
+    gdc = data["gdc_mu_o3a"][::ac_steps,::dc_steps]*np.pi
+    gac = data["gac_mu_o3a_ispline"][::ac_steps,::dc_steps]*np.pi
+    idc = data["idc_mu_o3a"][::ac_steps,::dc_steps]
+    iac = data["iac_mu_o3a"][::ac_steps,::dc_steps]
+    phase = data["phase_mu_o3a"][::ac_steps,::dc_steps]
+    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
@@ -191,12 +217,11 @@ def export_omega5():
             solver_tol_abs = 1e-10,
             voltage_branches = 4,
             truncation_order = 3,
-            include_Ga = False,
+            include_Ga = True,
             solve_integral_exactly = False,
             )
     bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
             | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
-            | DataManager.SOLVER_FLAGS["include_Ga"] \
             | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
             | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
             | DataManager.SOLVER_FLAGS["deleted"]
@@ -206,26 +231,26 @@ def export_omega5():
             omega = omega,
             vac_min = 0,
             vac_max = 165.372,
-            vac_num = 101,
+            vac_num = 51,
             vdc_min = 0,
             vdc_max = 165.372,
-            vdc_num = 101,
+            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(101):
+        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[10*i]/omega, vdc[10*i]/omega, np.pi*gdc[10*i], idc[10*i], iac[10*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[:,10*i]/omega, vdc[:,10*i]/omega, np.pi*gdc[:,10*i], idc[:,10*i], iac[:,10*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,
@@ -257,6 +282,7 @@ def export_interp(
         x_func = None,
         y_func = None,
         korder = 2,
+        spline_s = 2e-4,
         special_selection = None,
         o2_scale_x = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
         o2_scale_y = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
@@ -360,7 +386,7 @@ def export_interp(
             selected &= special_selection(thedata, order, method, padding, **kwargs)
         return thedata.loc[selected]
 
-    def gen_data(order, method, padding=False, s=5e-5, **kwargs):
+    def gen_data(order, method, padding=False, s=spline_s, **kwargs):
         suffix = "_p" if padding else ""
         reduced_data = selection(order, method, padding, **kwargs)
         settings.logger.info(f"Starting with {order}{suffix} {method}, found {reduced_data.shape[0]} data points")
@@ -409,9 +435,10 @@ def export_interp(
         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,
-                    reduced_data.ac_current_phase,
+                    (xy_data[0][sel], xy_data[1][sel]),
+                    reduced_data.ac_current_phase[sel],
                     mesh,
                     method="cubic")
         except:
@@ -571,6 +598,207 @@ def prepare_RGconvergence():
     print("AC voltage: " + " ".join("%.6g"%vac[i] for i in ac_indices))
 
 
+def adjust_argument_array(start, end, *funcs, initial_resolution=51, s1=0.01, maxit=5, scale_coef=-0.5, yscales=1):
+    x = np.linspace(start, end, initial_resolution)
+    ys = [func(x) for func in funcs]
+    try:
+        assert len(yscales) == len(funcs)
+    except:
+        yscales = len(funcs)*[1]
+    for _ in range(maxit):
+        y1s = [(y[1:] - y[:-1]) / (x[1:] - x[:-1]) for y in ys]
+        w1s = [(y1[1:] - y1[:-1]) / (1 + yscale*(y1[1:] + y1[:-1])) for y1,yscale in zip(y1s,yscales)]
+        if scale_coef != 0:
+            for i in range(len(funcs)):
+                w1s[i] *= (x[2:] - x[:-2])**scale_coef
+        w1max = max(w1.max() for w1 in w1s)
+        settings.logger.debug(f"adjust argument: {w1max} {s1}, iteration {_+1}")
+        if w1max <= s1:
+            return (x, *ys)
+        x_selection = np.zeros(x.shape, dtype=bool)
+        for w1 in w1s:
+            x_selection[1:-1] |= np.abs(w1) > s1
+        x_selection[:-2] |= x_selection[1:-1]
+        #x_selection[2:-1] |= np.abs(w2) > s2
+        new_x = (x[x_selection] + x[1:][x_selection[:-1]])/2
+        full_x = np.append(x, new_x)
+        sorting = np.argsort(full_x)
+        x = full_x[sorting]
+        for i, func in enumerate(funcs):
+            ys[i] = np.append(ys[i], func(new_x))[sorting]
+    return (x, *ys)
+
+
+
+def export_gapproximations_pgfplots(
+        omega = 16.5372,
+        vdc = 165.372,
+        vac1 = 16.5372,
+        vac2 = 165.372,
+        resolution = 500,
+        ):
+    parameters = dict(
+            d = 1e9,
+            include_Ga = True,
+            integral_method = -15,
+            solver_tol_rel = 1e-8,
+            solver_tol_abs = 1e-10,
+            xL = 0.5,
+            bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
+                    | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
+                    | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
+                    | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
+                    | DataManager.SOLVER_FLAGS["deleted"]
+                )
+
+    dm = DataManager()
+    data_omega = dm.list(vac=vac1, vdc=0.0, method="J", voltage_branches=0, **parameters).sort_values("omega")
+    data_omega.drop_duplicates("omega", inplace=True)
+    data_vac = dm.list(omega=omega, vdc=0.0, method="J", voltage_branches=0, **parameters).sort_values("vac")
+    data_vac.drop_duplicates("vac", inplace=True)
+    data_vdc = dm.list(vac=vac2, omega=omega, method="mu", voltage_branches=4, **parameters).sort_values("vdc")
+    data_vdc.drop_duplicates("vdc", inplace=True)
+    data_vac_vdc = dm.list(vdc=vdc, omega=omega, method="J", voltage_branches=0, resonant_dc_shift=10, **parameters)
+    # messy selection of data points for data_vac_vdc
+    sel = (data_vac_vdc.padding>10) & (data_vac_vdc.nmax>=40) & ((data_vac_vdc.version_minor==15) | (data_vac_vdc.vac>=165.372))
+    data_vac_vdc = data_vac_vdc[sel].sort_values("vac")
+    data_vac_vdc.drop_duplicates("vac", keep="last", inplace=True)
+    interp = interp_vac0(dm, voltage_branches=4, method="mu", **parameters)
+
+    def adiabatic(vac, vdc):
+        if np.all(vdc == 0):
+            result, error = quad_vec((lambda x: interp(vac*np.sin(x))), 0, np.pi/2, epsrel=1e-7, epsabs=1e-14)
+            result *= 2
+            error *= 2
+            assert error < 1e-6
+            return result
+        result, error = quad_vec((lambda x: interp(vdc+vac*np.sin(x))), -np.pi/2, np.pi/2, epsrel=1e-7, epsabs=1e-14)
+        assert error < 1e-6
+        return result
+
+    def kaminski_high_vac(vac):
+        """
+        Comparison to Kaminski et al PRB 62.8154 (2000)
+        Low frequency (adiabatic) limit
+        Eq. (77) in Kaminski et al 2000
+        """
+        g = 3*np.pi**2/(16*np.log(vac/TK_VOLTAGE)**2)
+        if isinstance(g, Number):
+            return g
+        g[vac <= TK_VOLTAGE] = np.nan
+        g[g > 1.2] = np.nan
+        return g
+
+    def kaminski_high_omega(vac, omega):
+        """
+        Comparison to Kaminski et al PRB 62.8154 (2000)
+        High frequency limit
+        Eqs. (72) - (74) in Kaminski et al 2000
+        """
+        # Eq. (72) in Kaminski et al 2000
+        decoherence_rate = (vac/TK_VOLTAGE)**2 / (np.pi*omega/TK_VOLTAGE*np.log(omega/TK_VOLTAGE)**2)
+        # Eq. (74) in Kaminski et al 2000
+        g = 1 - decoherence_rate
+        # Eq. (73) in Kaminski et al 2000
+        g[decoherence_rate>1] = 3*np.pi**2/(16*np.log(decoherence_rate[decoherence_rate>1])**2)
+        return g
+
+    def kaminski_low_energy(vac):
+        """
+        Comparison to Kaminski et al PRB 62.8154 (2000)
+        low energies, near equilibrium
+        Eq. (76) in Kaminski et al 2000
+        """
+        return 1 - 3/16 * (vac/TK_VOLTAGE)**2
+
+    common_adjust_kw = dict(s1=1e-4, initial_resolution=101, maxit=4, scale_coef=0.9)
+
+    # first plot: swipe omega at constant vac and vdc=0
+    g_pat_func = lambda log_omega: np.pi*photon_assisted_tunneling(omega=np.exp(log_omega), vdc=0, vac=vac1, interp=interp)
+    gdc_func = interp1d(np.log(data_omega.omega), np.pi*data_omega.dc_conductance, kind="cubic")
+    log_omega_arr, g_pat_arr, gdc_arr = adjust_argument_array(
+            np.log(data_omega.omega.iloc[0]), np.log(data_omega.omega.iloc[-1]),
+            g_pat_func, gdc_func,
+            **common_adjust_kw,
+            yscales=(10,10))
+    omega_arr = np.exp(log_omega_arr)
+    g_adiabatic_arr = np.ones_like(omega_arr) * adiabatic(vac1, 0)
+    g_adiabatic_arr[omega_arr > 2.5*TK_VOLTAGE] = np.nan
+    g_kaminski1_arr = np.ones_like(omega_arr) * kaminski_high_vac(vac1)
+    g_kaminski1_arr[omega_arr > 2.5*TK_VOLTAGE] = np.nan
+    g_kaminski2_arr = kaminski_high_omega(vac1, omega_arr)
+    g_kaminski2_arr[np.exp(vac1/(np.pi*omega_arr*TK_VOLTAGE)**0.5)*TK_VOLTAGE > omega_arr] = np.nan
+    np.savetxt("figdata/omega_vac5_vdc0.dat",
+               np.array([omega_arr/TK_VOLTAGE,
+                         gdc_arr,
+                         g_adiabatic_arr,
+                         g_pat_arr,
+                         g_kaminski1_arr,
+                         g_kaminski2_arr
+                         ]).T,
+                header="omega gdc gdc_adiabatic gdc_pat gdc_kaminski1 gdc_kaminski2",
+                comments = "")
+
+    # second plot: swipe vac at constant omega and vdc=0
+    g_pat_func = lambda log_vac: np.pi*photon_assisted_tunneling(omega=omega, vac=np.exp(log_vac), vdc=0, interp=interp)
+    gdc_func = interp1d(np.log(data_vac.vac), np.pi*data_vac.dc_conductance, kind="cubic")
+    log_vac_arr, g_pat_arr, gdc_arr = adjust_argument_array(
+            np.log(data_vac.vac.iloc[0]), np.log(data_vac.vac.iloc[-1]),
+            g_pat_func, gdc_func,
+            **common_adjust_kw,
+            yscales=(10,10))
+    vac_arr = np.exp(log_vac_arr)
+    g_adiabatic_arr = adiabatic(vac_arr, 0)
+    g_kaminski1_arr = kaminski_high_vac(vac_arr)
+    g_kaminski1_arr[g_kaminski1_arr>=1] = np.nan
+    g_kaminski2_arr = kaminski_high_omega(vac_arr, omega)
+    g_kaminski2_arr[np.exp(vac_arr/(np.pi*omega*TK_VOLTAGE)**0.5)*TK_VOLTAGE > omega] = np.nan
+    np.savetxt("figdata/vac_omega5_vdc0.dat",
+               np.array([vac_arr/TK_VOLTAGE,
+                         gdc_arr,
+                         g_adiabatic_arr,
+                         g_pat_arr,
+                         g_kaminski1_arr,
+                         g_kaminski2_arr
+                         ]).T,
+                header="vac gdc gdc_adiabatic gdc_pat gdc_kaminski1 gdc_kaminski2",
+                comments = "")
+
+    # third plot: swipe vdc at constant omega and vac
+    g_pat_func = lambda vdc: np.pi*photon_assisted_tunneling(omega=omega, vac=vac2, vdc=vdc, interp=interp)
+    gdc_func = interp1d(data_vdc.vdc, np.pi*data_vdc.dc_conductance, kind="cubic")
+    g_adiabatic_func = lambda vdc: adiabatic(vac2, vdc)
+    vdc_arr, g_pat_arr, gdc_arr, g_adiabatic_arr = adjust_argument_array(
+            data_vdc.vdc.iloc[0], data_vdc.vdc.iloc[-1],
+            g_pat_func, gdc_func, g_adiabatic_func,
+            **common_adjust_kw)
+    np.savetxt("figdata/vdc_omega5_vac10.dat",
+               np.array([vdc_arr/omega,
+                         gdc_arr,
+                         g_adiabatic_arr,
+                         g_pat_arr,
+                         ]).T,
+                header="vdc gdc gdc_adiabatic gdc_pat",
+                comments = "")
+
+    # fourth plot: swipe vac at constant omega and vdc
+    g_pat_func = lambda vac: np.pi*photon_assisted_tunneling(omega=omega, vac=vac, vdc=vdc, interp=interp)
+    gdc_func = interp1d(data_vac_vdc.vac, np.pi*data_vac_vdc.dc_conductance, kind="cubic")
+    g_adiabatic_func = lambda vac: adiabatic(vac, vdc)
+    vac_arr, g_pat_arr, gdc_arr, g_adiabatic_arr = adjust_argument_array(
+            data_vac_vdc.vac.iloc[0], data_vac_vdc.vac.iloc[-1],
+            g_pat_func, gdc_func, g_adiabatic_func,
+            **common_adjust_kw)
+    np.savetxt("figdata/vac_omega5_vdc10.dat",
+               np.array([vac_arr/omega,
+                         gdc_arr,
+                         g_adiabatic_arr,
+                         g_pat_arr,
+                         ]).T,
+                header="vac gdc gdc_adiabatic gdc_pat",
+                comments = "")
+
+
 if __name__ == "__main__":
     from sys import argv
     globals()[argv[1]]()
diff --git a/plot.py b/plot.py
index 012003d..b9c4879 100644
--- a/plot.py
+++ b/plot.py
@@ -18,6 +18,7 @@ from final_plots import filter_grid_data, photon_assisted_tunneling, gen_photon_
 from scipy.interpolate import bisplrep, bisplev, splrep, BSpline, griddata
 from scipy.optimize import curve_fit, leastsq
 from PIL import Image
+from imageio import imwrite
 import settings
 from logging import _levelToName
 from data_management import DataManager, KondoImport
@@ -1681,18 +1682,6 @@ def compare_RGeq_interp(
     #    contour_kwargs["alpha"] = 1
     #    ax.contour(vdc, vac, ref_data, missing_levels, **contour_kwargs)
 
-
-class GNorm(plt.Normalize):
-    def __init__(self, vmin=0.1, vmax=1, clip=False):
-        assert vmin > 0
-        assert vmax > vmin
-        self._vmin = vmin
-        self._vmax = vmax
-        self._clip = clip
-
-    def __call__(self, value, clip=None):
-        return (1 - self._vmin/value) / (1 - self._vmin/self._vmax)
-
 def GNorm(vmin, vmax):
     assert vmax > vmin > 0
     scale = 1 - vmin/vmax
@@ -2023,14 +2012,27 @@ def plot_RG2_commutator(
     s2 = ax2.scatter(vdc, vac, c=commutator01_max)
     s3 = ax3.scatter(vdc, vac, c=commutator00_max/g2_max)
     s4 = ax4.scatter(vdc, vac, c=commutator01_max/g2_max)
+
     fig.colorbar(s1, ax=ax1)
     fig.colorbar(s2, ax=ax2)
     fig.colorbar(s3, ax=ax3)
     fig.colorbar(s4, ax=ax4)
 
     fig, ax = plt.subplots()
-    s = ax.scatter(vdc/omega, vac/omega, c=commutator01_max/g2_max)
-    fig.colorbar(s, ax=ax)
+    vdc_mesh = np.linspace(0, 165.372, 201)
+    vac_mesh = np.linspace(0, 165.372, 201)
+    extent = (1.5*vdc_mesh[0] - 0.5*vdc_mesh[1], 1.5*vdc_mesh[-1] - 0.5*vdc_mesh[-2], 1.5*vac_mesh[0] - 0.5*vac_mesh[1], 1.5*vac_mesh[-1] - 0.5*vac_mesh[-2])
+    extent = tuple(v/omega for v in extent)
+    vdc_mesh, vac_mesh = np.meshgrid(vdc_mesh, vac_mesh)
+    interpolated = griddata((vdc, vac), commutator01_max/g2_max, (vdc_mesh, vac_mesh), method="cubic")
+    interp_max = interpolated.max()
+    print("Maximum of color scale:", interp_max)
+    norm = plt.Normalize(0, interp_max)
+    img = ax.imshow(interpolated, extent=extent, norm=norm, origin="lower")
+    imwrite(f"figdata/colorbar_viridis.png", np.array(0xffff*plt.cm.viridis(np.linspace(0, 1, 0x1000))[::-1,:3], dtype=np.uint16).reshape((0x1000,1,3)), format="PNG-FI")
+    imwrite(f"figdata/commutator_RG2_relative.png", np.array(0xffff*plt.cm.viridis(interpolated[::-1]/interp_max)[...,:3], dtype=np.uint16), format="PNG-FI")
+    s = ax.scatter(vdc/omega, vac/omega, c=commutator01_max/g2_max, norm=norm)
+    fig.colorbar(img, ax=ax)
     ax.set_xlabel(r"$V_\mathrm{avg}~(\Omega)$")
     ax.set_ylabel(r"$V_\mathrm{osc}~(\Omega)$")
     plt.show()
@@ -2178,6 +2180,7 @@ def compare_limits(
     data_vac = dm.list(omega=omega, vdc=0.0, method="J", voltage_branches=0, **parameters).sort_values("vac")
     data_vdc = dm.list(vac=vac_vdc, omega=omega_vdc, method="mu", voltage_branches=4, **parameters).sort_values("vdc")
     data_vac_vdc = dm.list(vdc=vdc, omega=omega_vac, **parameters).sort_values("vac")
+    #print(data_omega.shape, data_vac.shape, data_vdc.shape, data_vac_vdc.shape)
     interp = interp_vac0(dm, voltage_branches=4, method="mu", **parameters)
 
     def adiabatic(vac, vdc):
diff --git a/test/benchmark.py b/test/benchmark.py
index 982e6c0..9d7ebf0 100644
--- a/test/benchmark.py
+++ b/test/benchmark.py
@@ -5,7 +5,12 @@ Benchmarks for Kondo FRTRG solver
 import numpy as np
 import cProfile
 import pstats
-from pstats import SortKey
+try:
+       from pstats import SortKey
+except ImportError:
+       SortKey = lambda:None
+       SortKey.CUMULATIVE = None
+
 import settings
 from kondo import Kondo
 
diff --git a/tikz/commutator_RG2.tex b/tikz/commutator_RG2.tex
new file mode 100644
index 0000000..5491144
--- /dev/null
+++ b/tikz/commutator_RG2.tex
@@ -0,0 +1,39 @@
+\begin{tikzpicture}
+  \pgfplotsset{set layers = axis on top}
+  \begin{axis}[
+      name = img,
+      width = 5cm,
+      height = 5cm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0, ymax = 10,
+      xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\vac~(\Omega)$},
+      clip = true,
+    ]
+    \addplot graphics[xmin=-.05, xmax=10.05, ymin=-.05, ymax=10.05] {../figdata/commutator_RG2_relative.png};
+  \end{axis}
+  \begin{axis}[
+      name = colorbar,
+      at = (img.east),
+      anchor = west,
+      ymin = 0,
+      ymax = 0.0288735933,
+      xshift = 5mm,
+      enlargelimits = false,
+      axis on top,
+      clip = true,
+      width = .75em,
+      height = 48mm,
+      scale only axis,
+      tickpos = right,
+      tick align = outside,
+      xtick = \empty,
+      ytick = {0, 0.01, 0.02},
+      yticklabels = {$0$, $0.01$, $0.02$},
+      scaled ticks = false,
+      ylabel = {$\|[\chi^{-1},ZG^2_{LR}]\|\,/\,J$},
+    ]
+    \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=0.0288735933] {../figdata/colorbar_viridis.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/contour.tex b/tikz/contour.tex
index 0a6ca5e..a8ecaeb 100644
--- a/tikz/contour.tex
+++ b/tikz/contour.tex
@@ -1,24 +1,26 @@
+\definecolor{darkgreen}{HTML}{00c000}
 \begin{tikzpicture}
   \begin{axis}[
       width = 12cm,
       height = 10cm,
       xmin = 0, xmax = 10,
       ymin = 0, ymax = 10,
-      xlabel = {$\vac~(\Omega)$},
-      ylabel = {$\vdc~(\Omega)$},
+      xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\vac~(\Omega)$},
       every axis plot post/.append style = {
         line width=0.2pt,
         line cap=round,
         line join=round
       },
+      legend cell align = left,
       legend style = {
         anchor = south east,
         at = {(0.99,0.01)}
       }
     ]
-    \addplot[green] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_mu_o2.dat};
-    \addplot[green, forget plot] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_J_o2.dat};
-    \addplot[green, forget plot] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_J_o2_p.dat};
+    \addplot[darkgreen] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_mu_o2.dat};
+    \addplot[darkgreen, forget plot] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_J_o2.dat};
+    \addplot[darkgreen, forget plot] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_J_o2_p.dat};
     \addplot[blue] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_mu_o3.dat};
     \addplot[blue, forget plot] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_J_o3.dat};
     \addplot[blue, forget plot] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_J_o3_p.dat};
@@ -30,10 +32,14 @@
     \addplot[black, forget plot] table[x=vdc, y=vac, z=gdc]{../figdata/contour_gdc_J_o3p_p.dat};
 
     \legend{
-      $O(J^2)$,
-      {$O(J^3)$, $G^a=0$},
-      {$O(J^3)$, approx.~$I[X]$},
-      $O(J^3)$,
+      %$O(J^2)$,
+      %{$O(J^3)$, $G^a=0$},
+      %{$O(J^3)$, approx.~$I[X]$},
+      %$O(J^3)$,
+      leading order,
+      without $G^a$,
+      approx.~integral,
+      full RG equations
     }
   \end{axis}
 \end{tikzpicture}
diff --git a/tikz/convergence-o2.tex b/tikz/convergence-o2.tex
new file mode 100644
index 0000000..43f417f
--- /dev/null
+++ b/tikz/convergence-o2.tex
@@ -0,0 +1,164 @@
+\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=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o2_diff.png};
+    % vertical lines
+    \addplot[no marks, black] coordinates {(0,1.5) (10,1.5)};
+    \addplot[no marks, black] coordinates {(0,2.5) (10,2.5)};
+    \addplot[no marks, black] coordinates {(0,4) (10,4)};
+    \addplot[no marks, black] coordinates {(0,6.25) (10,6.25)};
+    \addplot[no marks, black] coordinates {(0,8.5) (10,8.5)};
+    \node[anchor=west] at (7.0,1.5) {\small △};
+    \node[anchor=west] at (7.5,2.5) {\small ▽};
+    \node[anchor=west] at (8.0,4) {\small □};
+    \node[anchor=west] at (8.5,6.25) {\small ○};
+    \node[anchor=west] at (9.0,8.5) {\small ◊};
+    % horizontal lines
+    \addplot[no marks, black] coordinates {(0,0) (0,10)};
+    \addplot[no marks, black] coordinates {(0.26666667,0) (0.26666667,10)};
+    \addplot[no marks, black] coordinates {(0.66666667,0) (0.66666667,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.26666667,7) {\small ◀};
+    \node[anchor=south] at (0.66666667,5) {\small ■};
+    \node[anchor=south] at (1,4) {\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,
+      },
+    ]
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.50}+0.48}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.50}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.50}+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.25}+0.12}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc8.50}}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc2.50}+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.25}+0.12}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc8.50}}] {../figdata/convergence_mu_o3p_ac.dat};
+    %\node[anchor=west] at (10.5,0.58) {\small$\vac=1.5\Omega$};
+    %\node[anchor=west] at (10.5,0.46) {\small$\vac=2.5\Omega$};
+    %\node[anchor=west] at (10.5,0.34) {\small$\vac=4\Omega$};
+    %\node[anchor=west] at (10.5,0.22) {\small$\vac=6.25\Omega$};
+    %\node[anchor=west] at (10.5,0.10) {\small$\vac=8.5\Omega$};
+    %\node[anchor=west] at (10,0.57) {\small 1};
+    %\node[anchor=west] at (10,0.45) {\small 2};
+    %\node[anchor=west] at (10,0.33) {\small 3};
+    %\node[anchor=west] at (10,0.22) {\small 4};
+    %\node[anchor=west] at (10,0.11) {\small 5};
+    \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.255) {\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,
+    ]
+    \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.27}+0.48}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.67}+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.27}+0.48}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.67}+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 = -0.0313315776,
+      xmax = 0.0313315776,
+      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.02, 0, 0.02},
+      xticklabels = {$-0.02$, $0$, $0.02$},
+      scaled ticks = false,
+      xlabel = {deviation $(2e^2/h)$},
+      xlabel shift = -15mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=-0.0313315776, xmax=0.0313315776] {../figdata/colorbar_seismic.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/convergence-o3.tex b/tikz/convergence-o3.tex
new file mode 100644
index 0000000..fab270c
--- /dev/null
+++ b/tikz/convergence-o3.tex
@@ -0,0 +1,164 @@
+\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=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o3_diff.png};
+    % vertical lines
+    \addplot[no marks, black] coordinates {(0,1.5) (10,1.5)};
+    \addplot[no marks, black] coordinates {(0,2.5) (10,2.5)};
+    \addplot[no marks, black] coordinates {(0,4) (10,4)};
+    \addplot[no marks, black] coordinates {(0,6.25) (10,6.25)};
+    \addplot[no marks, black] coordinates {(0,8.5) (10,8.5)};
+    \node[anchor=west] at (7.0,1.5) {\small △};
+    \node[anchor=west] at (7.5,2.5) {\small ▽};
+    \node[anchor=west] at (8.0,4) {\small □};
+    \node[anchor=west] at (8.5,6.25) {\small ○};
+    \node[anchor=west] at (9.0,8.5) {\small ◊};
+    % horizontal lines
+    \addplot[no marks, black] coordinates {(0,0) (0,10)};
+    \addplot[no marks, black] coordinates {(0.26666667,0) (0.26666667,10)};
+    \addplot[no marks, black] coordinates {(0.66666667,0) (0.66666667,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.26666667,7) {\small ◀};
+    \node[anchor=south] at (0.66666667,5) {\small ■};
+    \node[anchor=south] at (1,4) {\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,
+      },
+    ]
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.50}+0.48}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.50}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.50}+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.25}+0.12}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc8.50}}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc2.50}+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.25}+0.12}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc8.50}}] {../figdata/convergence_mu_o3p_ac.dat};
+    %\node[anchor=west] at (10.5,0.58) {\small$\vac=1.5\Omega$};
+    %\node[anchor=west] at (10.5,0.46) {\small$\vac=2.5\Omega$};
+    %\node[anchor=west] at (10.5,0.34) {\small$\vac=4\Omega$};
+    %\node[anchor=west] at (10.5,0.22) {\small$\vac=6.25\Omega$};
+    %\node[anchor=west] at (10.5,0.10) {\small$\vac=8.5\Omega$};
+    %\node[anchor=west] at (10,0.57) {\small 1};
+    %\node[anchor=west] at (10,0.45) {\small 2};
+    %\node[anchor=west] at (10,0.33) {\small 3};
+    %\node[anchor=west] at (10,0.22) {\small 4};
+    %\node[anchor=west] at (10,0.11) {\small 5};
+    \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.255) {\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,
+    ]
+    \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.27}+0.48}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.67}+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.27}+0.48}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.67}+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 = -0.0282119544,
+      xmax = 0.0282119544,
+      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.02, 0, 0.02},
+      xticklabels = {$-0.02$, $0$, $0.02$},
+      scaled ticks = false,
+      xlabel = {deviation $(2e^2/h)$},
+      xlabel shift = -15mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=-0.0282119544, xmax=0.0282119544] {../figdata/colorbar_seismic.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/convergence-o3a.tex b/tikz/convergence-o3a.tex
new file mode 100644
index 0000000..26d4835
--- /dev/null
+++ b/tikz/convergence-o3a.tex
@@ -0,0 +1,163 @@
+\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=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o3a_diff.png};
+    % vertical lines
+    \addplot[no marks, black] coordinates {(0,1.5) (10,1.5)};
+    \addplot[no marks, black] coordinates {(0,2.5) (10,2.5)};
+    \addplot[no marks, black] coordinates {(0,4) (10,4)};
+    \addplot[no marks, black] coordinates {(0,6.25) (10,6.25)};
+    \addplot[no marks, black] coordinates {(0,8.5) (10,8.5)};
+    \node[anchor=west] at (7.0,1.5) {\small △};
+    \node[anchor=west] at (7.5,2.5) {\small ▽};
+    \node[anchor=west] at (8.0,4) {\small □};
+    \node[anchor=west] at (8.5,6.25) {\small ○};
+    \node[anchor=west] at (9.0,8.5) {\small ◊};
+    % horizontal lines
+    \addplot[no marks, black] coordinates {(0,0) (0,10)};
+    \addplot[no marks, black] coordinates {(0.26666667,0) (0.26666667,10)};
+    \addplot[no marks, black] coordinates {(0.66666667,0) (0.66666667,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.26666667,7) {\small ◀};
+    \node[anchor=south] at (0.66666667,5) {\small ■};
+    \node[anchor=south] at (1,4) {\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,
+      },
+    ]
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.50}+0.48}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.50}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.50}+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.25}+0.12}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc8.50}}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc2.50}+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.25}+0.12}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc8.50}}] {../figdata/convergence_mu_o3p_ac.dat};
+    %\node[anchor=west] at (10.5,0.58) {\small$\vac=1.5\Omega$};
+    %\node[anchor=west] at (10.5,0.46) {\small$\vac=2.5\Omega$};
+    %\node[anchor=west] at (10.5,0.34) {\small$\vac=4\Omega$};
+    %\node[anchor=west] at (10.5,0.22) {\small$\vac=6.25\Omega$};
+    %\node[anchor=west] at (10.5,0.10) {\small$\vac=8.5\Omega$};
+    %\node[anchor=west] at (10,0.57) {\small 1};
+    %\node[anchor=west] at (10,0.45) {\small 2};
+    %\node[anchor=west] at (10,0.33) {\small 3};
+    %\node[anchor=west] at (10,0.22) {\small 4};
+    %\node[anchor=west] at (10,0.11) {\small 5};
+    \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.255) {\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,
+    ]
+    \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.27}+0.48}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.67}+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.27}+0.48}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.67}+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 = -0.003390334477,
+      xmax = 0.003390334477,
+      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.002, 0.002},
+      scaled ticks = false,
+      xlabel = {deviation $(2e^2/h)$},
+      xlabel shift = -15.85mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=-0.003390334477, xmax=0.003390334477] {../figdata/colorbar_seismic.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/omega5-3d.tex b/tikz/omega5-3d.tex
index 7492d2c..f53c84c 100644
--- a/tikz/omega5-3d.tex
+++ b/tikz/omega5-3d.tex
@@ -2,10 +2,10 @@
 \begin{tikzpicture}
   \begin{axis}[
       name = G3d,
-      view={110}{22},
-      width = 7cm,
-      height = 7cm,
-      zmin = 0.075, zmax = 80,
+      view={110}{15},
+      width = 5.5cm,
+      height = 6.5cm,
+      zmin = 0.075, zmax = 30,
       ymin = 0, ymax = 10,
       xmin = 0, xmax = 10,
       zmode = log,
@@ -18,9 +18,9 @@
       every tick label/.append style = {font=\small},
       label shift = {-.8ex},
       zlabel shift = {-1.5ex},
-      ztick = {0.1, 1, 8, 80},
+      ztick = {0.1, 1, 3, 30},
       zticklabels = {$0.1$, $1$, $0.1$, $1$},
-      minor ztick = {0.08,0.09,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 5.6,6.4,7.2,16,32,40,48,56,64,72},
+      minor ztick = {0.08,0.09,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 2.1,2.4,2.7,3,6,9,12,15,18,21,24,27},
       %zlabel style = {at = (ticklabel cs:0.65)},
       xtick = {0,5,10},
       xticklabels = {$\hspace*{-.5em}0$, $5$, $10$},
@@ -30,23 +30,52 @@
       clip = false,
       scale only axis,
     ]
-    \addplot3[gray, no marks, line width=0.2pt] coordinates {(10,0,0.3) (0,0,0.3) (0,10,0.3)};
+    \addplot3[gray, no marks, line width=0.2pt] coordinates {(10,0,0.1) (0,0,0.1) (0,10,0.1)};
     \addplot3[
       line width = 0.18pt,
       surf,
-      point meta min = -1.15,
-      point meta max = -0.4,
-      point meta = {log10(\thisrow{gac})},
-    ] table [x expr={\thisrow{vac}/5.000007155}, y expr={\thisrow{vdc}/5.000007155}, z expr={\thisrow{gac}}]{../figdata/vdc_vac_omega16.5372_interp.dat};
-    \addplot3[gray, no marks, line width=0.2pt] coordinates {(10,0,8) (0,0,8) (0,10,8)};
+      point meta min = 0,
+      point meta max = 1,
+      point meta = {(1-0.18/(\thisrow{gac}+0.1))/0.72},
+    ] table [x=vac, y=vdc, z expr={\thisrow{gac}}]{../figdata/omega5_interp.dat};
+    \addplot3[gray, no marks, line width=0.2pt] coordinates {(10,0,0.1) (10,10,0.1) (0,10,0.1)};
+    \addplot3[gray, no marks, line width=0.2pt] coordinates {(10,0,3) (0,0,3) (0,10,3)};
     \addplot3[
       line width = 0.18pt,
       surf,
-      point meta min = -1.15,
-      point meta max = -0.4,
-      point meta = {log10(\thisrow{gdc_g})},
-    ] table [x expr={\thisrow{vac}/5.000007155}, y expr={\thisrow{vdc}/5.000007155}, z expr={80*\thisrow{gdc_g}}]{../figdata/vdc_vac_omega16.5372_interp.dat};
-    \addplot3[gray, no marks, line width=0.2pt] coordinates {(10,0,8) (10,10,8) (0,10,8)};
+      point meta min = 0,
+      point meta max = 1,
+      point meta = {(1-0.18/(\thisrow{gdc}+0.1))/0.72},
+    ] table [x=vac, y=vdc, z expr={30*\thisrow{gdc}}]{../figdata/omega5_interp.dat};
+    \addplot3[gray, no marks, line width=0.2pt] coordinates {(10,0,3) (10,10,3) (0,10,3)};
     \node[black] at (axis description cs:0.75,0.82) {$\Omega=4.8\tkv$};
   \end{axis}
+  \begin{axis}[
+      at = (G3d.north east),
+      view={0}{90},
+      anchor = north west,
+      width = 0.75em,
+      height = 65mm,
+      yshift = -1.5ex,
+      colormap/viridis,
+      scale only axis,
+      xmin = 0, xmax = 1,
+      xshift = 1.2em,
+      ymin = 0.08, ymax = 1,
+      axis on top,
+      xtick = \empty,
+      tickpos = right,
+      tick align = outside,
+    ]
+    \addplot3[
+      surf,
+      domain = 0:1,
+      domain y = 0.08:1,
+      samples = 2,
+      samples y = 100,
+      shader = interp,
+      mesh/color input = colormap,
+      point meta = {(1-0.18/(y+0.1))/0.72}
+    ] {(1-0.18/(y+0.1))/0.72};
+  \end{axis}
 \end{tikzpicture}
diff --git a/tikz/overview.tex b/tikz/overview.tex
index cf5a448..bf44b66 100644
--- a/tikz/overview.tex
+++ b/tikz/overview.tex
@@ -1,7 +1,7 @@
 \begin{tikzpicture}
   \begin{axis}[
-      width=11cm,
-      height=11cm,
+      width=9.5cm,
+      height=9.5cm,
       scale only axis,
       xmin=0, xmax=10,
       ymin=0, ymax=10,
@@ -11,6 +11,8 @@
       xlabel = {$\vdc/\Omega$},
       ylabel = {$\vac/\Omega$},
       zlabel = {$\Omega~(\tkv)$},
+      axis lines* = left,
+      axis on top,
     ]
     \addplot3[surf] graphics[
       points={
-- 
GitLab