Skip to content
Snippets Groups Projects
Select Git revision
  • 08af1bfa6b8bbb08dca246fbc45655c75f132707
  • master default protected
  • b24_tutorial
  • 2.4.1
  • 1.3.0
  • 1.2.0
6 results

setup.py

Blame
  • final_plots.py 85.74 KiB
    #!/usr/bin/env python3
    
    # Copyright 2022 Valentin Bruch <valentin.bruch@rwth-aachen.de>
    # License: MIT
    """
    Kondo FRTRG, generate high-quality plots for publication
    """
    
    import os
    import itertools
    import scipy.constants as sc
    import pandas as pd
    import matplotlib.pyplot as plt
    import matplotlib.colors as mplcolors
    from matplotlib.widgets import Slider
    import argparse
    import numpy as np
    from scipy.interpolate import bisplrep, bisplev, splrep, splev, BSpline, griddata, interp1d
    from scipy.special import jn
    from imageio import imwrite
    from scipy.integrate import quad, quad_vec
    from scipy.optimize import newton, curve_fit, leastsq
    import settings
    from data_management import DataManager, KondoImport
    from kondo import gen_init_matrix
    from numbers import Number
    from gen_pulse_data import fourier_coef_gauss_symmetric
    from plot_pulses import load_all_pulses_full, integrate_ft
    from kondo import solveTV0_scalar
    
    # 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
    # G(V=Tk)=G(V=0)/2=e²/h. The ratio Tk/Tkrg is:
    
    #TK_VOLTAGE = 3.44249 # for rtol=1e-10, atol=1e-12, voltage_branches=10
    TK_VOLTAGE = 3.4425351 # for rtol=1e-9, atol=1e-11, voltage_branches=10
    #TK_VOLTAGE = 3.44334 # with correction
    TK_VOLTAGE_O3P = 3.458524 # for rtol=1e-9, atol=1e-11, voltage_branches=10
    #TK_VOLTAGE_O3P = 3.4593 # with correction
    TK_VOLTAGE_ORDER2 = 10.1368086 # for rtol=1e-9, atol=1e-11, voltage_branches=10
    #TK_VOLTAGE_ORDER2 = 10.13754 # with correction
    
    def save_overview(
            omega = 16.5372,
            vdc_res = 501,
            vac_res = 501,
            vdc_max = 165.372,
            vac_max = 165.372,
            method = "mu",
            d = 1e9,
            xL = 0.5,
            solver_tol_rel = 1e-8,
            solver_tol_abs = 1e-10,
            lazy_inverse_factor = 0.25,
            voltage_branches = 4,
            s_g = 1e-4,
            s_idc = 5e-6,
            s_iac = 5e-6,
            filename = "figdata/vdc_vac_omega16.5372_interp.npz",
            **kwargs
            ):
        """
        3d plot of ac and dc differential conductance and current as function of Vac and Vdc.
        """
        dm = DataManager()
        data = dm.list(omega=omega, vdc=None, vac=None, method=method, d=d, xL=xL, solver_tol_abs=solver_tol_abs, solver_tol_rel=solver_tol_rel, voltage_branches=voltage_branches, lazy_inverse_factor=lazy_inverse_factor, **kwargs)
        bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
                | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
                | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
                | DataManager.SOLVER_FLAGS["deleted"]
        filters = (data.vac < 1.2*vac_max) \
                & (data.vdc < 1.2*vdc_max) \
                & np.isfinite(data.dc_conductance) \
                & (data.solver_flags & bad_flags == 0)
        data = data.loc[filters]
    
        print("Interpolate Gdc", flush=True)
        gdc_tck = bisplrep(data.vac, data.vdc, data.dc_conductance, s=s_g, kx=3, ky=3)
        print("Interpolate Idc", flush=True)
        idc_tck = bisplrep(data.vac, data.vdc, data.dc_current, s=s_idc, kx=3, ky=3)
        print("Interpolate Iac", flush=True)
        iac_tck = bisplrep(data.vac, data.vdc, data.ac_current_abs, s=s_iac, kx=3, ky=3)
        print("done.")
        vac_arr = np.linspace(0, min(vac_max, data.vac.max()), vac_res)
        vdc_arr = np.linspace(0, min(vdc_max, data.vdc.max()), vdc_res)
        gdc_g_interp = bisplev(vac_arr, vdc_arr, gdc_tck).T
        idc_interp   = bisplev(vac_arr, vdc_arr, idc_tck).T
        iac_interp   = bisplev(vac_arr, vdc_arr, iac_tck).T
        gdc_i_interp = bisplev(vac_arr, vdc_arr, idc_tck, dy=1).T
        gac_interp   = 2*bisplev(vac_arr, vdc_arr, iac_tck, dx=1).T
        np.savez(filename, vac=vac_arr/TK_VOLTAGE, vdc=vdc_arr/TK_VOLTAGE, gdc_g=np.pi*gdc_g_interp, gdc_i=np.pi*gdc_i_interp, gac=np.pi*gac_interp, idc=idc_interp/TK_VOLTAGE, iac=iac_interp/TK_VOLTAGE)
        #vac_data, vdc_data = np.meshgrid(vac_arr, vdc_arr)
        #with open(filename, 'w') as file:
        #    file.write("vdc,vac,gdc_g,gdc_i,gac,idc,iac\n")
        #    np.savetxt(file, np.array([vdc_data/TK_VOLTAGE, vac_data/TK_VOLTAGE, np.pi*gdc_g_interp, np.pi*gdc_i_interp, np.pi*gac_interp, idc_interp/TK_VOLTAGE, iac_interp/TK_VOLTAGE]).reshape((7,-1)).T, fmt="%.9e", delimiter=",")
        #    #for i in range(vdc_res):
        #    #    file.write('\n')
        #    #    np.savetxt(file, np.array([vdc_data[i]/omega, vac_data[i]/omega, np.pi*gdc_g_interp[i], np.pi*gdc_i_interp[i], np.pi*gac_interp[i], idc_interp[i], iac_interp[i]]).T)
    
    
    def interp_vac0(dm, xL=0.5, include_Ga=True, truncation_order=3, integral_method=-15, **trashparams):
        """Return differential conductance G as function of Vdc at Vac=0"""
        data = dm.list(vac=0, omega=0, include_Ga=include_Ga, d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, nmax=0, voltage_branches=4, xL=xL, truncation_order=truncation_order, integral_method=integral_method)
        data = data.sort_values("vdc")
        vdc = np.concatenate((-data.vdc[:0:-1], data.vdc))
        gdc = np.concatenate((data.dc_conductance[:0:-1], data.dc_conductance))
        return interp1d(vdc, gdc, kind="cubic")
    
    
    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)
        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)
            except ValueError:
                pass
        return result
    
    
    def gen_photon_assisted_tunneling(dm, fourier_coef, omega, vdc, nmax=50, **parameters):
        """
        Return differential conductance G calculated from generalized
        photon assisted tunneling at given parameters.
    
        check:
        >>> gen_photon_assisted_tunneling(dm, (vac/2,0), omega, vdc, ...) \
        >>>         == photon_assisted_tunneling(dm, omega, vdc, vac, ...)
        """
        result = np.zeros(np.broadcast_shapes(np.shape(omega), np.shape(vdc)), dtype=np.complex128)
        interp = interp_vac0(dm, **parameters)
        f0 = gen_init_matrix(nmax, *(f/omega for f in fourier_coef))
        for n in range(-nmax, nmax+1):
            try:
                result += (abs(f0[nmax+n,nmax])**2 + abs(f0[nmax,nmax+n])**2)/2 * interp(vdc+n*omega)
            except ValueError:
                pass
        assert (np.abs(result.imag) < 1e-9).all()
        return result.real
    
    
    def filter_grid_data(
            data : pd.DataFrame,
            omega = 16.5372,
            vac_min = 0,
            vac_max = 165.372,
            vac_num = 101,
            vdc_min = 0,
            vdc_max = 165.372,
            vdc_num = 101,
            v_tol = 1e-6,
            ):
        vac_step = (vac_max - vac_min) / (vac_num - 1)
        vdc_step = (vdc_max - vdc_min) / (vdc_num - 1)
        data = data.loc[(data.vac >= vac_min - v_tol) & (data.vdc >= vdc_min - v_tol) & (data.vac <= vac_max + v_tol) & (data.vdc <= vdc_max + v_tol) & np.isfinite(data.dc_conductance)]
        grid_data = data.loc[(np.abs(((data.vac - vac_min + v_tol) % vac_step) - v_tol)  < v_tol) & (np.abs(((data.vdc - vdc_min + v_tol) % vdc_step) - v_tol)  < v_tol)]
        data = grid_data.sort_values(["vac", "vdc"])
        vac_arr = np.linspace(vac_min, vac_max, vac_num)
        vdc_arr = np.linspace(vdc_min, vdc_max, vdc_num)
        gdc_arr = np.empty((vac_num, vdc_num), dtype=np.float64)
        idc_arr = np.empty((vac_num, vdc_num), dtype=np.float64)
        iac_arr = np.empty((vac_num, vdc_num), dtype=np.float64)
        phase_arr = np.empty((vac_num, vdc_num), dtype=np.float64)
        gdc_arr.fill(np.nan)
        idc_arr.fill(np.nan)
        iac_arr.fill(np.nan)
        phase_arr.fill(np.nan)
        lower_index = 0
        for i, vac in enumerate(vac_arr):
            upper_index = data.vac.searchsorted(vac + v_tol)
            indices = vdc_arr.searchsorted(data.vdc[lower_index:upper_index] - v_tol)
            indices = indices[np.abs(vdc_arr[indices] - data.vdc[lower_index:upper_index]) < v_tol]
            gdc_arr[i, indices] = data.dc_conductance[lower_index:upper_index]
            idc_arr[i, indices] = data.dc_current[lower_index:upper_index]
            iac_arr[i, indices] = data.ac_current_abs[lower_index:upper_index]
            phase_arr[i, indices] = data.ac_current_phase[lower_index:upper_index]
            lower_index = upper_index
        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
        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,
            fixed_value,
            x_parameter,
            y_parameter,
            x_arr,
            y_arr,
            x_sort_parameter = None,
            y_sort_parameter = None,
            x_func = None,
            y_func = None,
            kxorder = 2,
            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
            ):
        """
        Export interpolated data. Parameters:
    
            filename: string, file path
            fixed_parameter: name of fixed parameter (e.g. "omega")
            fixed_value: value of fixed parameter (e.g. 16.5372)
            x_parameter: name of x axis parameter (e.g. "vdc")
            y_parameter: name of y axis parameter (e.g. "vac")
            x_sort_parameter: name of x axis parameter for sorting
            y_sort_parameter: name of y axis parameter for sorting
            x_func: function to obtain x value from data
            y_func: function to obtain y value from data
            kxorder: order of spline interpolation in x direction
            kyorder: order of spline interpolation in y direction
            special_selection: function for filtering input data. Arguments:
                data, order, method, padding, **kwargs
    
        Data are stored with units e=hbar=kB=Tkrg=1.
        All results are exported as 2d arrays of the same shape.
        The parameters vdc and vac are included as 2d arrays.
        Results have the naming convention
    
            {observable}_{method}_{order}{padding}{suffix}
    
        observable:
            gdc   DC differential conductance, may have suffix "_ispline"
            gac   AC differential conductance, must have suffix "_ispline"
            idc   average current, may have suffix "_spline"
            iac   oscillating current, may have suffix "_spline"
            phase AC phase
            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
    
        padding:
            ""    no Floquet matrix extrapolation
            "_p"  Floquet matrix extrapolation to mitigate truncation effects
    
        method:
            mu    no unitary transformation, oscillating chemical potentials
            J     include oscillating voltage in coupling by unitary transformation
    
        suffix
            ""    interpolation is done using griddata
            "_spline"   interpolation is done using spline
            "_ispline"  data is derived from spline interpolation for the
                        current (dc or ac)
        """
        assert isinstance(filename, str)
        if x_func is None:
            x_func = lambda data: data[x_parameter]
        if y_func is None:
            y_func = lambda data: data[y_parameter]
        x_mesh, y_mesh = np.meshgrid(x_arr, y_arr)
        mesh = (x_mesh, y_mesh)
        dm = DataManager()
        data = dm.list(**{fixed_parameter:fixed_value}, **parameters)
        global_bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
                | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
                | DataManager.SOLVER_FLAGS["deleted"]
        data = data.loc[data.solver_flags & global_bad_flags == 0]
        data = data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
        required_flags = dict(
                o2 = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
                o3 = 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"],
                o3a = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
                o3p = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
                )
    
        o2_data = dm.list(**{fixed_parameter:fixed_value*o2_scale_fixed}, **parameters, min_version=(14,15,-1,-1))
        o2_data = o2_data.loc[o2_data.solver_flags & (global_bad_flags | bad_flags["o2"] | required_flags["o2"]) == required_flags["o2"]]
        o2_data = o2_data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
    
        def selection(order, method, padding=False, **kwargs):
            thedata = o2_data if order == "o2" else data
            selected = (thedata.solver_flags & required_flags[order] == required_flags[order]) \
                    & (thedata.solver_flags & bad_flags[order] == 0) \
                    & (thedata.method == method)
            if padding:
                selected &= (thedata["padding"] > 0) | (thedata["nmax"] == 0)
            else:
                selected &= thedata["padding"] == 0
            for key, value in kwargs.items():
                selected &= (thedata[key] == value)
            if special_selection is not None:
                selected &= special_selection(thedata, order, method, padding, **kwargs)
            return thedata.loc[selected]
    
        def gen_data(order, method, padding=False, s=spline_s, extend_vdc=10, extend_vac=10, **kwargs):
            suffix = "_p" if padding else ""
            reduced_data = selection(order, method, padding, **kwargs)
            settings.logger.info(f"Starting with {order}{suffix} {method}, found {reduced_data.shape[0]} data points")
            if extend_vdc or extend_vac:
                extended_data = [reduced_data]
                if extend_vdc:
                    dc_mirrored = reduced_data[reduced_data.vdc<extend_vdc].copy()
                    dc_mirrored.vdc *= -1
                    dc_mirrored.dc_current *= -1
                    extended_data.append(dc_mirrored)
                    del dc_mirrored
                if extend_vac:
                    ac_mirrored = reduced_data[reduced_data.vac<extend_vac].copy()
                    ac_mirrored.vac *= -1
                    ac_mirrored.ac_current_phase += np.pi
                    extended_data.append(ac_mirrored)
                    if extend_vdc:
                        acdc_mirrored = ac_mirrored[ac_mirrored.vdc<extend_vdc].copy()
                        acdc_mirrored.vdc *= -1
                        acdc_mirrored.dc_current *= -1
                        extended_data.append(acdc_mirrored)
                        del acdc_mirrored
                    del ac_mirrored
                reduced_data = pd.concat(extended_data)
            results = {}
            if reduced_data.shape[0] < 100:
                return results
            xy_data = (x_func(reduced_data), y_func(reduced_data))
            if order == "o2":
                xy_data = (xy_data[0]/o2_scale_x, xy_data[1]/o2_scale_y)
            try:
                results[f"gdc_{method}_{order}{suffix}"] = griddata(
                        xy_data,
                        reduced_data.dc_conductance,
                        mesh,
                        method="cubic")
                #gdc_tck = bisplrep(reduced_data[y], reduced_data[x], reduced_data.dc_conductance, s=s, kx=kxorder, ky=kyorder)
                #results[f"gdc_{method}_{order}{suffix}"] = bisplev(y_arr, x_arr, gdc_tck)
            except:
                settings.logger.exception(f"in griddata interpolation of Gdc for {order}{suffix} {method}")
            try:
                results[f"idc_{method}_{order}{suffix}"] = griddata(
                        xy_data,
                        reduced_data.dc_current,
                        mesh,
                        method="cubic")
            except:
                settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method}")
            try:
                idc_tck = bisplrep(*xy_data[::-1], reduced_data.dc_current, s=s, kx=kxorder, ky=kyorder)
                results[f"idc_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, idc_tck)
                results[f"gdc_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, idc_tck, dy=1)
            except:
                settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method}")
            try:
                results[f"iac_{method}_{order}{suffix}"] = griddata(
                        xy_data,
                        2*reduced_data.ac_current_abs,
                        mesh,
                        method="cubic")
            except:
                settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
            try:
                iac_tck = bisplrep(*xy_data[::-1], 2*reduced_data.ac_current_abs, s=s, kx=kxorder, ky=kyorder)
                results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, iac_tck)
                results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
            except:
                settings.logger.exception(f"in spline interpolation of Iac for {order}{suffix} {method}")
            try:
                sel = reduced_data.vac > 0
                results[f"phase_{method}_{order}{suffix}"] = griddata(
                        (xy_data[0][sel], xy_data[1][sel]),
                        reduced_data.ac_current_phase[sel],
                        mesh,
                        method="cubic")
            except:
                settings.logger.exception(f"in griddata interpolation of phase for {order}{suffix} {method}")
            try:
                results[f"gamma_{method}_{order}{suffix}"] = griddata(
                        xy_data,
                        reduced_data.gamma,
                        mesh,
                        method="cubic")
            except:
                settings.logger.exception(f"in griddata interpolation of gamma for {order}{suffix} {method}")
            return results
    
        np.savez(
                filename,
                **{fixed_parameter:fixed_value, x_parameter:x_mesh, y_parameter:y_mesh},
                **gen_data("o2", "mu"),
                **gen_data("o3", "mu"),
                **gen_data("o3a", "mu"),
                **gen_data("o3p", "mu", integral_method=-1),
                **gen_data("o2", "J", False),
                **gen_data("o3", "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("o3a", "J", True),
                **gen_data("o3p", "J", True, integral_method=-1),
                )
    
    
    def export_interp_relative_o3a(
            filename,
            fixed_parameter,
            fixed_value,
            x_parameter,
            y_parameter,
            x_arr,
            y_arr,
            x_sort_parameter = None,
            y_sort_parameter = None,
            x_func = None,
            y_func = None,
            kxorder = 2,
            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
            ):
        """
        Export interpolated data. Parameters:
    
            filename: string, file path
            fixed_parameter: name of fixed parameter (e.g. "omega")
            fixed_value: value of fixed parameter (e.g. 16.5372)
            x_parameter: name of x axis parameter (e.g. "vdc")
            y_parameter: name of y axis parameter (e.g. "vac")
            x_sort_parameter: name of x axis parameter for sorting
            y_sort_parameter: name of y axis parameter for sorting
            x_func: function to obtain x value from data
            y_func: function to obtain y value from data
            kxorder: order of spline interpolation in x direction
            kyorder: order of spline interpolation in y direction
            special_selection: function for filtering input data. Arguments:
                data, order, method, padding, **kwargs
    
        Data are stored with units e=hbar=kB=Tkrg=1.
        All results are exported as 2d arrays of the same shape.
        The parameters vdc and vac are included as 2d arrays.
        Results have the naming convention
    
            {observable}_{method}_{order}{padding}{suffix}
    
        observable:
            gdc   DC differential conductance, may have suffix "_ispline"
            gac   AC differential conductance, must have suffix "_ispline"
            idc   average current, may have suffix "_spline"
            iac   oscillating current, may have suffix "_spline"
            phase AC phase
            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
    
        padding:
            ""    no Floquet matrix extrapolation
            "_p"  Floquet matrix extrapolation to mitigate truncation effects
    
        method:
            mu    no unitary transformation, oscillating chemical potentials
            J     include oscillating voltage in coupling by unitary transformation
    
        suffix
            ""    interpolation is done using griddata
            "_spline"   interpolation is done using spline
            "_ispline"  data is derived from spline interpolation for the
                        current (dc or ac)
        """
        assert isinstance(filename, str)
        if x_func is None:
            x_func = lambda data: data[x_parameter]
        if y_func is None:
            y_func = lambda data: data[y_parameter]
        x_mesh, y_mesh = np.meshgrid(x_arr, y_arr)
        mesh = (x_mesh, y_mesh)
        dm = DataManager()
        data = dm.list(**{fixed_parameter:fixed_value}, **parameters)
        global_bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
                | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
                | DataManager.SOLVER_FLAGS["deleted"]
        data = data.loc[data.solver_flags & global_bad_flags == 0]
        data = data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
        data.vdc = np.round(data.vdc, 9)
        data.vac = np.round(data.vac, 9)
        required_flags = dict(
                o3 = 0,
                o3a = DataManager.SOLVER_FLAGS["include_Ga"],
                o3p = DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
                o3a_vb7 = DataManager.SOLVER_FLAGS["include_Ga"],
                )
        bad_flags = dict(
                o3 = 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"],
                o3a_vb7 = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
                )
        voltage_branches = dict(
                o3 = 4,
                o3a = 4,
                o3p = 4,
                o3a_vb7 = 7,
                )
        ref_data = data.loc[(data.solver_flags & (required_flags["o3a"] | bad_flags["o3a"]) == required_flags["o3a"]) & (data.method == "mu")]
        merged = pd.merge(
                data,
                ref_data,
                how="inner",
                on=["vdc","vac","d","omega","energy_im","energy_re","solver_tol_rel","solver_tol_abs","xL","lazy_inverse_factor","resonant_dc_shift"],
                suffixes=("", "_ref"))
    
        def selection(order, method, padding=False, voltage_branches=4, **kwargs):
            selected = (merged.solver_flags & required_flags[order] == required_flags[order]) \
                    & (merged.solver_flags & bad_flags[order] == 0) \
                    & (merged.method == method) \
                    & (merged.voltage_branches == voltage_branches)
            if padding:
                selected &= (merged.padding > 0) | (merged.nmax == 0)
            else:
                selected &= merged.padding == 0
            for key, value in kwargs.items():
                selected &= (merged[key] == value)
            if special_selection is not None:
                selected &= special_selection(merged, order, method, padding, **kwargs)
            return merged.loc[selected]
    
        def gen_data(order, method, padding=False, s=spline_s, extend_vdc=10, extend_vac=10, **kwargs):
            suffix = "_p" if padding else ""
            reduced_data = selection(order, method, padding, voltage_branches[order], **kwargs)
            settings.logger.info(f"Starting with {order}{suffix} {method}, found {reduced_data.shape[0]} data points")
            results = {}
            if reduced_data.shape[0] < 100:
                return results
            if extend_vdc or extend_vac:
                extended_data = [reduced_data]
                if extend_vdc:
                    dc_mirrored = reduced_data[reduced_data.vdc<extend_vdc].copy()
                    dc_mirrored.vdc *= -1
                    dc_mirrored.dc_current *= -1
                    dc_mirrored.dc_current_ref *= -1
                    extended_data.append(dc_mirrored)
                    del dc_mirrored
                if extend_vac:
                    ac_mirrored = reduced_data[reduced_data.vac<extend_vac].copy()
                    ac_mirrored.vac *= -1
                    ac_mirrored.ac_current_phase += np.pi
                    ac_mirrored.ac_current_phase_ref += np.pi
                    extended_data.append(ac_mirrored)
                    if extend_vdc:
                        acdc_mirrored = ac_mirrored[ac_mirrored.vdc<extend_vdc].copy()
                        acdc_mirrored.vdc *= -1
                        acdc_mirrored.dc_current *= -1
                        acdc_mirrored.dc_current_ref *= -1
                        extended_data.append(acdc_mirrored)
                        del acdc_mirrored
                    del ac_mirrored
                reduced_data = pd.concat(extended_data)
            xy_data = (x_func(reduced_data), y_func(reduced_data))
            try:
                results[f"gdc_{method}_{order}{suffix}"] = griddata(
                        xy_data,
                        reduced_data.dc_conductance - reduced_data.dc_conductance_ref,
                        mesh,
                        method="cubic")
            except:
                settings.logger.exception(f"in griddata interpolation of Gdc for {order}{suffix} {method}")
            try:
                results[f"idc_{method}_{order}{suffix}"] = griddata(
                        xy_data,
                        reduced_data.dc_current - reduced_data.dc_current_ref,
                        mesh,
                        method="cubic")
            except:
                settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method}")
            try:
                idc_tck = bisplrep(*xy_data[::-1], reduced_data.dc_current - reduced_data.dc_current_ref, s=s, kx=kxorder, ky=kyorder)
                results[f"idc_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, idc_tck)
                results[f"gdc_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, idc_tck, dy=1)
            except:
                settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method}")
            try:
                results[f"iac_{method}_{order}{suffix}"] = griddata(
                        xy_data,
                        2*(reduced_data.ac_current_abs - reduced_data.ac_current_abs_ref),
                        mesh,
                        method="cubic")
            except:
                settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
            try:
                iac_tck = bisplrep(*xy_data[::-1], 2*(reduced_data.ac_current_abs - reduced_data.ac_current_abs_ref), s=s, kx=kxorder, ky=kyorder)
                results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, iac_tck)
                results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
            except:
                settings.logger.exception(f"in spline interpolation of Iac for {order}{suffix} {method}")
            try:
                sel = reduced_data.vac > 0
                results[f"phase_{method}_{order}{suffix}"] = griddata(
                        (xy_data[0][sel], xy_data[1][sel]),
                        reduced_data.ac_current_phase[sel] - reduced_data.ac_current_phase_ref[sel],
                        mesh,
                        method="cubic")
            except:
                settings.logger.exception(f"in griddata interpolation of phase for {order}{suffix} {method}")
            try:
                results[f"gamma_{method}_{order}{suffix}"] = griddata(
                        xy_data,
                        reduced_data.gamma - reduced_data.gamma_ref,
                        mesh,
                        method="cubic")
            except:
                settings.logger.exception(f"in griddata interpolation of gamma for {order}{suffix} {method}")
            return results
    
        np.savez(
                filename,
                **{fixed_parameter:fixed_value, x_parameter:x_mesh, y_parameter:y_mesh},
                **gen_data("o3", "mu"),
                **gen_data("o3a_vb7", "mu"),
                **gen_data("o3p", "mu", integral_method=-1),
                **gen_data("o3", "J", False),
                **gen_data("o3a", "J", False),
                **gen_data("o3p", "J", False, integral_method=-1),
                **gen_data("o3", "J", True),
                **gen_data("o3a", "J", True),
                **gen_data("o3p", "J", True, integral_method=-1),
                )
    
    
    def export_omega5_interp(
            filename = "figdata/omega5_interp.npz",
            omega = 16.5372,
            vdc_min = 0,
            vdc_max = 165.372,
            vac_min = 0,
            vac_max = 165.372,
            dc_res = 301,
            ac_res = 201,
            vdc_max_J = 82.686,
            vac_max_J = 82.686,
            voltage_branches = 4,
            ):
        """
        Export interpolated data for fixed frequency omega.
        """
        def special_selection(data, order, method, padding, **kwargs):
            if method == "J":
                return (data["vdc"] < vdc_max_J + 1e-3) & (data["vac"] < vac_max_J + 1e-3)
            return True
        yarr = np.linspace(vac_min, vac_max, ac_res)
        yarr[-1] -= 1e-10
        yarr[ac_res//2] -= 1e-10
        export_interp(filename,
                      fixed_parameter = "omega",
                      fixed_value = omega,
                      x_parameter = "vdc",
                      y_parameter = "vac",
                      x_arr = np.linspace(vdc_min, vdc_max, dc_res),
                      y_arr = yarr,
                      kyorder = 3,
                      special_selection = special_selection,
                      spline_s = 2e-5,
                      d = 1e9,
                      solver_tol_rel = 1e-8,
                      solver_tol_abs = 1e-10,
                      voltage_branches = voltage_branches,
                      xL = 0.5)
    
    
    def export_omega5_interp_deviation(
            filename = "figdata/omega5_interp_deviation.npz",
            omega = 16.5372,
            vdc_min = 0,
            vdc_max = 165.372,
            vac_min = 0,
            vac_max = 165.372,
            dc_res = 301,
            ac_res = 201,
            vdc_max_J = 82.686,
            vac_max_J = 82.686,
            ):
        """
        Export interpolated data for fixed frequency omega.
        """
        def special_selection(data, order, method, padding, **kwargs):
            if method == "J":
                return (data["vdc"] < vdc_max_J + 1e-3) & (data["vac"] < vac_max_J + 1e-3)
            return True
        yarr = np.linspace(vac_min, vac_max, ac_res)
        yarr[-1] -= 1e-10
        yarr[ac_res//2] -= 1e-10
        export_interp_relative_o3a(filename,
                      fixed_parameter = "omega",
                      fixed_value = omega,
                      x_parameter = "vdc",
                      y_parameter = "vac",
                      x_arr = np.linspace(vdc_min, vdc_max, dc_res),
                      y_arr = yarr,
                      kyorder = 3,
                      special_selection = special_selection,
                      spline_s = 2e-5,
                      d = 1e9,
                      solver_tol_rel = 1e-8,
                      solver_tol_abs = 1e-10,
                      xL = 0.5)
    
    
    def export_vdc0_interp(
            filename = "figdata/vdc0_interp.npz",
            omega_min = 0.1,
            omega_max = 40,
            vac_omega_min = 0.01,
            vac_omega_max = 10,
            omega_res = 201,
            vac_res = 201,
            korder = 2,
            ):
        """
        Export interpolated data for fixed frequency omega.
        """
        parameters = dict(d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, voltage_branches=0, xL=0.5)
        export_interp(filename,
                      fixed_parameter = "vdc",
                      fixed_value = 0,
                      x_parameter = "omega",
                      y_parameter = "vac_omega",
                      y_sort_parameter = "vac",
                      x_arr = np.linspace(omega_min, omega_max, omega_res),
                      y_arr = np.linspace(vac_omega_min, vac_omega_max, vac_res),
                      y_func = (lambda data: data["vac"]/data["omega"]),
                      korder = korder,
                      special_selection = (lambda data, *args, **kwargs: data["omega"] > 0),
                      **parameters)
    
    
    def prepare_plotly_csv():
        dm = DataManager()
        # General overview
        reduction_dict = dict(omega="omega", vdc="vdc", vac="vac", dc_conductance="g", dc_current="idc", ac_current_abs="iac", ac_current_phase="ac_phase")
        data = dm.list(d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, truncation_order=3)
        bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
                | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
                | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
                | DataManager.SOLVER_FLAGS["include_Ga"] \
                | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
                | DataManager.SOLVER_FLAGS["deleted"]
        data = data.loc[data.solver_flags & bad_flags == 0]
        data = data.rename(columns=reduction_dict)
        data.g *= np.pi
        data.omega /= TK_VOLTAGE
        data.vdc /= TK_VOLTAGE
        data.vac /= TK_VOLTAGE
        data.idc /= TK_VOLTAGE
        data.iac /= TK_VOLTAGE
        data = data.sort_values(["vac","vdc","omega"])
        #data.to_csv("html/full.csv", columns=reduction_dict.values())
        # omega = 5Tk
        # TODO: this is outdated!
        for omega, name in zip((16.5372, 9.2159791, 5.8206184, 7.1271), ("omega5", "compare_bruhat18a", "compare_bruhat18b", "compare_kogan04")):
            reduced = data.loc[np.isclose(data.omega*TK_VOLTAGE, omega) & (data.method == "mu") & (data.voltage_branches == 4)]
            reduced.to_csv(f"html/{name}.csv", columns=reduction_dict.values())
        # vdc = 0
        reduced = data.loc[(data.vdc == 0) & (data.method != "mu")]
        reduced.to_csv("html/vdc0.csv", columns=reduction_dict.values())
    
    
    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")
        def export_img(name, array):
            imwrite(f"figdata/{name}.png", np.array(0xffff*real_cmap(array[::-1]), dtype=np.uint16), format="PNG-FI")
    
        data = np.load("figdata/omega5_interp.npz")
        data_vb7 = np.load("figdata/omega5_interp_vb7.npz")
        data_rel = np.load("figdata/omega5_interp_deviation.npz")
    
        # O2
        omega = 16.5372
        omega_o2 = 48.695054
        vdc = data["vdc"][0] / omega
        vac = data["vac"][:,0] / omega
    
        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")):
            gdiff = data[f"gdc_{u}_{method}"] - data[f"gdc_mu_{ref}"]
            if u != "mu":
                method += "_" + u
            normalize = np.nanmax(np.abs(gdiff))
            export_img(f"convergence_{method}_diff", gdiff/normalize)
            print(f"Color scale for {method} difference: {np.pi*normalize}")
            grel = gdiff/data[f"gdc_mu_{ref}"]
            normalize = np.nanmax(np.abs(grel))
            export_img(f"convergence_{method}_relative", grel/normalize)
            print(f"Color scale for {method} relative: {normalize}")
    
        for method, u, sign in (("o3p","mu",-1), ("o3a_vb7","mu",-1), ("o3a","J",1), ("o3a_p","J",1)):
            gdiff = sign*data_rel[f"gdc_{u}_{method}"]
            if u != "mu":
                method += "_" + u
            normalize = np.nanmax(np.abs(gdiff))
            export_img(f"convergence_{method}_diff_fromdiff", gdiff/normalize)
            print(f"Color scale for {method} difference (interpolate in the end): {np.pi*normalize}")
            grel = gdiff/data[f"gdc_mu_o3a"]
            normalize = np.nanmax(np.abs(grel))
            export_img(f"convergence_{method}_relative_fromdiff", grel/normalize)
            print(f"Color scale for {method} relative (interpolate in the end): {normalize}")
    
        export_img(f"convergence_o3a_J_relative_fromdiff_pnorm", data_rel[f"gdc_J_o3a"]/normalize)
    
        g1rel = data_rel["gdc_J_o3a"] / data["gdc_mu_o3a"]
        g2rel = data_rel["gdc_J_o3a_p"] / data["gdc_mu_o3a"]
        g2rel[:16,:8] = np.nan
        normalize = max(np.nanmax(np.abs(g1rel)), np.nanmax(np.abs(g2rel)))
        g1img = np.array(0xffff*real_cmap(g1rel[::-1]/normalize), dtype=np.uint16)
        g2img = np.array(0xffff*real_cmap(g2rel[::-1]/normalize), dtype=np.uint16)
        g2img[-16:,:8] = 0x8000
        imwrite(f"figdata/convergence_o3a_J_relative_fromdiff_maskednorm.png", g1img, format="PNG-FI")
        imwrite(f"figdata/convergence_o3a_J_p_relative_fromdiff_masked.png", g2img, format="PNG-FI")
        print(f"Color scale for o3a J relative (interpolate in the end, masked): {normalize}")
    
        grel = -data_rel["gdc_mu_o3a_vb7"]/data["gdc_mu_o3a"]
        normalize = max(np.nanmax(np.abs(grel[16:])), np.nanmax(np.abs(grel[:16,8:])))
        export_img(f"convergence_o3a_vb7_relative_fromdiff_masked", grel/normalize)
        print(f"Color scale for vb7 masked relative (interpolate in the end): {normalize}")
    
        gdiff = data["gdc_mu_o3a"] - data_vb7["gdc_mu_o3a"]
        normalize = np.nanmax(np.abs(gdiff))
        export_img("convergence_o3a_vb7_diff", gdiff/normalize)
        print(f"Color scale for o3a_vb7 difference: {np.pi*normalize}")
        grel = gdiff/data_vb7["gdc_mu_o3a"]
        normalize = np.nanmax(np.abs(grel))
        export_img("convergence_o3a_vb7_relative", grel/normalize)
        print(f"Color scale for o3a_vb7 relative: {normalize}")
    
        for method, u  in (("o2","mu"), ("o3","mu"), ("o3a","mu"), ("o3p","mu"), ("o3a_ispline","mu")):
            if u != "mu":
                method += "_" + u
            g = np.pi*data[f"gdc_{u}_{method}"]
            with open(f"figdata/convergence_{u}_{method}_dc.dat", "w") as file:
                file.write("vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices) + "\n")
                np.savetxt(file, np.array([vac, *(g[:,i] for i in dc_indices)]).T)
            with open(f"figdata/convergence_{u}_{method}_ac.dat", "w") as file:
                file.write("vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices) + "\n")
                np.savetxt(file, np.array([vdc, *(g[i] for i in ac_indices)]).T)
        g = np.pi*data_vb7[f"gdc_mu_o3a"]
        with open(f"figdata/convergence_mu_o3a_vb7_dc.dat", "w") as file:
            file.write("vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices) + "\n")
            np.savetxt(file, np.array([vac, *(g[:,i] for i in dc_indices)]).T)
        with open(f"figdata/convergence_mu_o3a_vb7_ac.dat", "w") as file:
            file.write("vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices) + "\n")
            np.savetxt(file, np.array([vdc, *(g[i] for i in ac_indices)]).T)
    
        print("DC voltage: " + " ".join("%.6g"%vdc[i] for i in dc_indices))
        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]
        if isinstance(yscales, Number):
            yscales = len(funcs)*[yscales]
        else:
            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,
            ):
        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 = "",
                    fmt = "%.6g")
    
        # 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 = "",
                    fmt = "%.6g")
    
    
    def export_asymmetry_pgfplots(
            omega = 16.5372,
            vac = 26.45952,
            xL = (0.5, 0.4, 0.3, 0.2, 0.1, 0.001),
            vdc_max = 82.686,
            ):
        dm = DataManager()
        data = dm.list(
                omega = omega,
                vac = vac,
                d = 1e9,
                method = "mu",
                voltage_branches = 4,
                include_Ga = True,
                integral_method = -15,
                solver_tol_rel = 1e-8,
                solver_tol_abs = 1e-10,
                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"]
                ).sort_values("vdc")
        data.drop_duplicates(("vdc", "xL"), inplace=True)
        gdc_funcs = []
        for x in xL:
            sel_data = data[np.abs(data.xL - x) < 1e-6]
            gdc_funcs.append(interp1d(sel_data.vdc, np.pi/(4*x*(1-x))*sel_data.dc_conductance, kind="cubic"))
        vdc_arr, *gdc_arrs = adjust_argument_array(0, vdc_max, *gdc_funcs, s1=1e-3, yscales=[0.1 for x in xL], initial_resolution=101, maxit=3)
        np.savetxt("figdata/asymmetry.dat",
                   np.array([vdc_arr/omega,
                             *gdc_arrs
                             ]).T,
                   header="vdc " + " ".join(f"xL{x:.3g}" for x in xL),
                   comments = "",
                   fmt = "%.6g")
    
    
    def export_kogan04_pgfplots():
        """
        In the paper the numbers Tk≃300mK and Ω≃2Tk are given.
        These correspond to 7.4182 and 6.8851, respectively.
        """
        #omega = 3.1
        #xL = 0.05
        omega = 3.5
        xL = 0.166667
        #omega = 4
        #xL = 0.2
        dm = DataManager()
        # adjusted by factor 1.4:
        vac_omega = np.array([0.72881, 1.13091, 1.50788, 1.6838,  3.61891])
        parameters = dict(
                d = 1e9,
                method = "mu",
                voltage_branches = 4,
                include_Ga = True,
                integral_method = -15,
                solver_tol_rel = 1e-8,
                solver_tol_abs = 1e-10,
                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"]
                )
        data = dm.list(omega=omega, **parameters).sort_values(["xL", "vac", "vdc"])
        gdc_funcs = []
        gdc_funcs_asym = []
        #prefactor_sym = 0.113*np.pi
        #prefactor_asym = 0.62*np.pi
        #g_shift_sym = 0.023
        #g_shift_asym = 0.022
        #prefactor_sym = 0.1*np.pi
        #g_shift_sym = 0.015
        # for omega=3.5, xL=0.166667:
        prefactor_sym = 0.096*np.pi
        g_shift_sym = 0.02
        # for omega=4, xL=0.2:
        #prefactor_sym = 0.093*np.pi
        #g_shift_sym = 0.024
        #prefactor_sym = 0.0915*np.pi
        #g_shift_sym = 0.0245
        prefactor_asym = prefactor_sym / (4*xL*(1-xL))
        g_shift_asym = g_shift_sym
        for v in vac_omega:
            data_sel = data[(np.abs(data.xL-0.5) < 1e-6) & (np.abs(data.vac-v*omega)<1e-6)]
            try:
                gdc_funcs.append(interp1d(data_sel.vdc, g_shift_sym + prefactor_sym*data_sel.dc_conductance, kind="cubic", bounds_error=False))
            except ValueError:
                gdc_funcs.append(lambda x: np.nan*np.empty_like(x))
            data_sel = data[(np.abs(data.xL-xL) < 1e-6) & (np.abs(data.vac-v*omega)<1e-6)]
            try:
                gdc_funcs_asym.append(interp1d(data_sel.vdc, g_shift_asym + prefactor_asym*data_sel.dc_conductance, kind="cubic", bounds_error=False))
            except ValueError:
                gdc_funcs_asym.append(lambda x: np.nan*np.empty_like(x))
        vdc, *arrays = adjust_argument_array(0, 25.5, *gdc_funcs, *gdc_funcs_asym, s1=5e-4, maxit=4, initial_resolution=51)
        vdc = np.append(-vdc[:0:-1], vdc)
    
        interp = interp_vac0(dm, **parameters)
        g_pat = np.array([g_shift_sym + prefactor_sym*photon_assisted_tunneling(omega=omega, vdc=vdc, vac=vac, interp=interp) for vac in omega*vac_omega])
    
        np.savetxt(
                "figdata/kogan04.dat",
                np.array([vdc/omega, *(np.append(a[:0:-1], a) for a in arrays), *g_pat]).T,
                header = "vdc " + " ".join(f"g{v:.3g}" for v in vac_omega) + " " + " ".join(f"g{v:.3g}asym" for v in vac_omega) + " " + " ".join(f"g{v:.3g}pat" for v in vac_omega),
                fmt = "%.6g",
                comments = "")
    
    
    def export_bruhat18_pgfplots(
            plot_s = 2e-3,
            init_res = 41,
            ):
        dm = DataManager()
        omega1 = 4.14
        omega2 = 2.61
        prefactor = 0.282*np.pi
        g_shift = 0.062
        xL = 0.333333
        prefactor_asym = prefactor/(4*xL*(1-xL))
        # including correction factors:
        vac_omega1 = np.array([0.90611, 1.132637, 1.359164, 1.585692, 1.812219, 2.038747, 2.265274, 2.491802, 2.718329, 3.397911])
        vac_omega2 = np.array([0.431208, 0.862416, 1.293624, 1.724832, 2.15604, 2.587248, 3.018457, 3.449665, 3.880873, 4.743289, 6.468121])
        parameters = dict(
                d = 1e9,
                method = "mu",
                voltage_branches = 4,
                include_Ga = True,
                integral_method = -15,
                solver_tol_rel = 1e-8,
                solver_tol_abs = 1e-10,
                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"],
                )
        gdc_funcs1 = []
        gdc_funcs2 = []
        gdc_funcs1_asym = []
        # Collect data
        data1 = dm.list(omega=omega1, **parameters).sort_values(["xL", "vac", "vdc"])
        for v in vac_omega1:
            data_sel = data1[(np.abs(data1.xL-0.5) < 1e-6) & (np.abs(data1.vac-v*omega1)<1e-6)]
            try:
                gdc_funcs1.append(interp1d(data_sel.vdc, g_shift + prefactor*data_sel.dc_conductance, kind="cubic", bounds_error=False))
            except ValueError:
                gdc_funcs1.append(lambda x: np.nan*np.empty_like(x))
            data_sel = data1[(np.abs(data1.xL-0.333333) < 1e-6) & (np.abs(data1.vac-v*omega1)<1e-6)]
            try:
                gdc_funcs1_asym.append(interp1d(data_sel.vdc, g_shift + prefactor_asym*data_sel.dc_conductance, kind="cubic", bounds_error=False))
            except ValueError:
                gdc_funcs1_asym.append(lambda x: np.nan*np.empty_like(x))
        data2 = dm.list(omega=omega2, xL=0.5, **parameters).sort_values(["vac", "vdc"])
        for v in vac_omega2:
            data_sel = data2[np.abs(data2.vac-v*omega2)<1e-6]
            try:
                gdc_funcs2.append(interp1d(data_sel.vdc, g_shift + prefactor*data_sel.dc_conductance, kind="cubic", bounds_error=False))
            except ValueError:
                gdc_funcs2.append(lambda x: np.nan*np.empty_like(x))
    
        interp = interp_vac0(dm, **parameters)
        get_pat = lambda omega, vac: (lambda vdc: g_shift + prefactor*photon_assisted_tunneling(omega=omega, vdc=vdc, vac=vac, interp=interp))
        g1_pat_funcs = [get_pat(omega1, omega1*vac) for vac in vac_omega1]
        g2_pat_funcs = [get_pat(omega2, omega2*vac) for vac in vac_omega2]
    
        vdc1, *arrays1 = adjust_argument_array(0, 21.5, *gdc_funcs1, *g1_pat_funcs, s1=plot_s, maxit=4, initial_resolution=init_res)
        vdc1 = np.append(-vdc1[:0:-1], vdc1) / omega1
        n1 = len(vac_omega1)
        g1 = np.array([np.append(a[:0:-1], a) for a in arrays1[:n1]])
        g1_pat = np.array([np.append(a[:0:-1], a) for a in arrays1[n1:]])
        vdc2, *arrays2 = adjust_argument_array(0, 21.5, *gdc_funcs2, *g2_pat_funcs, s1=plot_s, maxit=4, initial_resolution=init_res)
        vdc2 = np.append(-vdc2[:0:-1], vdc2) / omega2
        n2 = len(vac_omega2)
        g2 = np.array([np.append(a[:0:-1], a) for a in arrays2[:n2]])
        g2_pat = np.array([np.append(a[:0:-1], a) for a in arrays2[n2:]])
    
        vdc1_asym, *g1_asym = adjust_argument_array(0, 21.5, *gdc_funcs1_asym, s1=plot_s, maxit=4, initial_resolution=init_res)
        vdc1_asym = np.append(-vdc1_asym[:0:-1], vdc1_asym) / omega1
        g1_asym = np.array([np.append(a[:0:-1], a) for a in g1_asym])
    
        # collect experimental data
        exp_data2 = np.genfromtxt("exp_data/KondoAC_Freq12GHz.dat", names=True)
        vdc_omega_exp2 = exp_data2["Vsd_mV"] * sc.eV/(12e12*sc.h)
        g_exp2_20 = exp_data2["Vac20"]
        # get background
        g_exp2_bg = g_exp2_20 - gdc_funcs2[0](np.abs(vdc_omega_exp2)*omega2)
        sel = ~np.isnan(g_exp2_bg)
        g_exp2_bg_spl = splrep(vdc_omega_exp2[sel], g_exp2_bg[sel], s=1.1e-2)
        bg2 = np.array([quad_vec((lambda t: splev(vdc2+vac*np.cos(t), g_exp2_bg_spl, ext=3)), 0, np.pi)[0]/np.pi for vac in vac_omega2])
        bg1 = np.array([quad_vec((lambda t: splev((vdc1+vac*np.cos(t))*omega1/omega2, g_exp2_bg_spl, ext=3)), 0, np.pi)[0]/np.pi for vac in vac_omega1])
    
        np.savetxt("figdata/bruhat18_19GHz.dat",
                   np.array([vdc1, *g1, *bg1, *(g1 + bg1), *g1_pat, *(g1_pat + bg1)]).T,
                   header = " ".join((
                       "vdc",
                       " ".join(f"g{v:.3g}" for v in vac_omega1),
                       " ".join(f"bg{v:.3g}" for v in vac_omega1),
                       " ".join(f"gbg{v:.3g}" for v in vac_omega1),
                       " ".join(f"gpat{v:.3g}" for v in vac_omega1),
                       " ".join(f"gpatbg{v:.3g}" for v in vac_omega1),
                       )),
                   fmt = "%.6g",
                   comments = "")
        np.savetxt("figdata/bruhat18_12GHz.dat",
                   np.array([vdc2, *g2, *bg2, *(g2 + bg2), *g2_pat, *(g2_pat + bg2)]).T,
                   header = " ".join((
                       "vdc",
                       " ".join(f"g{v:.3g}" for v in vac_omega2),
                       " ".join(f"bg{v:.3g}" for v in vac_omega2),
                       " ".join(f"gbg{v:.3g}" for v in vac_omega2),
                       " ".join(f"gpat{v:.3g}" for v in vac_omega2),
                       " ".join(f"gpatbg{v:.3g}" for v in vac_omega2),
                       )),
                   fmt = "%.6g",
                   comments = "")
        np.savetxt("figdata/bruhat18_19GHz_asym.dat",
                   np.array([vdc1_asym, *g1_asym]).T,
                   header = "vdc " + " ".join(f"g{v:.3g}" for v in vac_omega1),
                   fmt = "%.6g",
                   comments = "")
    
    
    def export_pulse_current_pgfplots(omega=1.5, pulse_duration=0.01):
        dm = DataManager()
        i_coef = []
        g_coef = []
        for phase in (0.25, 0.5, 0.75, 1):
            vdc, fourier_coef = fourier_coef_gauss_symmetric(10, omega, pulse_duration, None, phase)
            kondo, = dm.load_all(
                    omega = omega,
                    d = 1e9,
                    method = "J",
                    voltage_branches = 0,
                    include_Ga = True,
                    integral_method = -15,
                    solver_tol_rel = 1e-8,
                    solver_tol_abs = 1e-10,
                    has_fourier_coef = True,
                    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"],
                    fourier_coef = fourier_coef
                    )
            nmax = kondo.nmax
            i_coef.append(kondo.gammaL[nmax:,nmax]/TK_VOLTAGE)
            g_coef.append(np.pi*kondo.deltaGammaL[nmax:,nmax])
            del kondo
        t_prefactor = 2*np.pi/omega * TK_VOLTAGE
        t = np.concatenate((np.linspace(-0.01, -0.05/t_prefactor, 10, endpoint=False),
                            np.linspace(-0.05/t_prefactor, 0.08/t_prefactor, 130, endpoint=False),
                            np.linspace(0.08/t_prefactor, 0.5/t_prefactor, 42, endpoint=False),
                            np.linspace(0.5/t_prefactor, 3.33/t_prefactor, 40, endpoint=False),
                            np.linspace(3.33/t_prefactor, 0.49, 41)))
        vdc, fourier_coef = fourier_coef_gauss_symmetric(1000, omega, pulse_duration, None, 1)
        u = vdc + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(fourier_coef, 1))
        u /= u.max()
        i_arrs = [-coef[0].real + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(coef)) for coef in i_coef]
        g_arrs = [-coef[0].real + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(coef)) for coef in g_coef]
        t *= t_prefactor
    
        array = np.array([t, u, *i_arrs, *g_arrs]).T
        np.savetxt("figdata/pulse_current_full.dat",
                   array,
                   header = "t u " + " ".join(f"i{i}" for i in range(1,5)) + " " + " ".join(f"g{i}" for i in range(1,5)),
                   fmt = "%.6g",
                   comments = "")
    
        np.savetxt("figdata/pulse_current_zoom.dat",
                   array[10:141],
                   header = "t u " + " ".join(f"i{i}" for i in range(1,5)) + " " + " ".join(f"g{i}" for i in range(1,5)),
                   fmt = "%.6g",
                   comments = "")
    
        np.savetxt("figdata/pulse_current_tail.dat",
                   array[140:223],
                   header = "t u " + " ".join(f"i{i}" for i in range(1,5)) + " " + " ".join(f"g{i}" for i in range(1,5)),
                   fmt = "%.6g",
                   comments = "")
    
    
    def export_pulse_charge_pgfplots(omega=1.5):
        dm = DataManager()
        duration_data = [0.04]
        phase_data = [0]
        charge_data = [0]
        trec_data = [np.nan]
        _h5files = set()
        parameters = dict(
                    d = 1e9,
                    omega = omega,
                    method = "J",
                    voltage_branches = 0,
                    include_Ga = True,
                    integral_method = -15,
                    solver_tol_rel = 1e-8,
                    solver_tol_abs = 1e-10,
                    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"],
                    pulse_type = "gauss_symmetric",
                )
        parameters2 = parameters.copy()
        parameters2.update(omega=2, solver_tol_rel=2e-8, solver_tol_abs=2e-10)
        for (row, kondo) in itertools.chain(
                load_all_pulses_full(dm, **parameters),
                load_all_pulses_full(dm, **parameters2),
                ):
            ref_time = min(4*row.pulse_duration*row.omega/(2*np.pi), 0.15)
            nmax = kondo.nmax
            charge_data.append(integrate_ft(-ref_time, 0.5-ref_time, kondo.gammaL[nmax-1::-1,nmax])/row.omega*2*np.pi)
            _h5files.add(kondo._h5file)
            duration_data.append(row.pulse_duration)
            phase_data.append(row.pulse_phase)
            nmax = kondo.nmax
            g_coef = np.pi*kondo.deltaGammaL[nmax::-1,nmax]
            g_coef[1:] *= 2
            aux = 2j*np.pi*np.arange(g_coef.size)
            func = lambda t: (g_coef*np.exp(t*aux)).real.sum() - 0.8
            func_d1 = lambda t: (aux*g_coef*np.exp(t*aux)).real.sum()
            func_d2 = lambda t: (aux**2*g_coef*np.exp(t*aux)).real.sum()
            #func = lambda t: sum((c*np.exp(t*2j*np.pi*n)).real for n,c in enumerate(g_coef)) - 1
            trec = np.nan
            for init in np.array([1.2, 2, 4, 8, 16]) * row.pulse_duration:
                try:
                    trec = newton(func, init, fprime=func_d1, fprime2=func_d2)
                    assert trec > 0.8*row.pulse_duration and trec < 0.45
                    break
                except (RuntimeError, AssertionError):
                    pass
                trec = np.nan
            trec_data.append(TK_VOLTAGE*2*np.pi/row.omega * trec)
            #print(trec, func(trec), charge_data[-1], row.pulse_phase, row.pulse_duration)
        for file in _h5files:
            file.close()
    
        phase_data = np.array(phase_data)
        duration_data = np.array(duration_data)
        charge_data = np.array(charge_data)
        trec_data = np.array(trec_data)
    
        order = np.argsort(duration_data)
        duration_data = duration_data[order]
        charge_data = charge_data[order]
        phase_data = phase_data[order]
        trec_data = trec_data[order]
    
        for phase in (0.25, 0.5, 0.75, 1):
            sel = np.abs(phase_data - phase) < 1e-5
            np.savetxt(f"figdata/pulse_charge_duration_p{phase:.2g}.dat",
                       np.array([duration_data[sel]*TK_VOLTAGE, charge_data[sel], trec_data[sel]]).T,
                       header = "t q trec",
                       fmt = "%.6g",
                       comments = "")
    
        sel = np.abs(duration_data - 0.04) < 1e-6
        np.savetxt("figdata/pulse_charge_phase.dat",
                   np.array([phase_data[sel], charge_data[sel], trec_data[sel]]).T,
                   header = "phase q trec",
                   fmt = "%.6g",
                   comments = "")
    
        q_duration_funcs = []
        trec_duration_funcs = []
        for phase in (0.25, 0.5, 0.75, 1):
            sel = np.abs(phase_data - phase) < 1e-5
            q_duration_funcs.append(interp1d(duration_data[sel], charge_data[sel], kind="cubic", bounds_error=False, fill_value=np.nan))
            sel *= np.isfinite(trec_data)
            trec_duration_funcs.append(interp1d(duration_data[sel], trec_data[sel], kind="cubic", bounds_error=False, fill_value=np.nan))
    
        t, *arrs = adjust_argument_array(0.01, 0.2, *q_duration_funcs, *trec_duration_funcs, s1=1e-3, maxit=3)
        np.savetxt("figdata/pulse_charge_duration_interp.dat",
                   np.array([t*TK_VOLTAGE, *arrs]).T,
                   header = "t " + " ".join(f"q{i}" for i in range(1,5)) + " " + " ".join(f"trec{i}" for i in range(1,5)),
                   fmt = "%.6g",
                   comments = "")
    
        sel = np.abs(duration_data - 0.04) < 1e-5
        q_duration_func = interp1d(phase_data[sel], charge_data[sel], kind="cubic", bounds_error=False, fill_value=np.nan)
        sel *= np.isfinite(trec_data)
        trec_duration_func = interp1d(phase_data[sel], trec_data[sel], kind="cubic", bounds_error=False, fill_value=np.nan)
    
        phases, *arrs = adjust_argument_array(0.0, 2.8, q_duration_func, trec_duration_func, s1=1e-3, maxit=3)
        np.savetxt("figdata/pulse_charge_phase_interp.dat",
                   np.array([phases, *arrs]).T,
                   header = "phase q trec",
                   fmt = "%.6g",
                   comments = "")
    
    
    def export_harmonic_modes(
            #omega = 3.4425351,
            omega = 68.850702,
            #omega = 16.5372,
            n_phase_crossings = 10,
            ):
        label = "-omega20"
        vac_omega = np.array([5, 10, 20, 40, 80])
        dm = DataManager()
        data = dm.list(
                    omega = omega,
                    d = 1e9,
                    vdc = 0,
                    method = "J",
                    voltage_branches = 0,
                    include_Ga = True,
                    integral_method = -15,
                    solver_tol_rel = 1e-8,
                    solver_tol_abs = 1e-10,
                    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"],
                    )
        current = []
        t = np.linspace(-0.5, 0.5, 201)
        current = []
        phase = []
        valid_vac_omega = []
        phase_crossing_times = []
        phase_crossing_currents = []
        for vac in vac_omega:
            sel = np.abs(data.vac-vac*omega) < 1e-6
            if sel.sum() != 1:
                continue
            row = data[sel].iloc[-1]
            filename = os.path.join(row.dirname, row.basename)
            kondo, = KondoImport.read_from_h5(filename, row.hash)
            nmax = row.nmax
            coef = kondo.gammaL[nmax::-1,nmax]
            iarr = 2*sum((c*np.exp(2j*np.pi*t*n)).real for n,c in enumerate(coef)) - coef.real[0]
            iarr /= TK_VOLTAGE
            current.append(iarr)
            pc_times = np.arcsin(2*np.pi/vac*np.arange(n_phase_crossings)-1)/(2*np.pi)
            pc_currents = 2*sum((c*np.exp(2j*np.pi*pc_times*n)).real for n,c in enumerate(coef)) - coef.real[0]
            pc_currents /= TK_VOLTAGE
            phase_crossing_times.append(pc_times)
            phase_crossing_currents.append(pc_currents)
            phase.append(vac/(2*np.pi)*(1+np.sin(2*np.pi*t)))
            valid_vac_omega.append(vac)
            coef /= coef[1]
            np.savetxt(
                    f"figdata/harmonic_modes_vac{vac:.3g}{label}.dat",
                    np.array([np.arange(1, coef.size, 2), np.abs(coef[1::2])]).T,
                    header = "n i",
                    comments = "",
                    fmt = "%.6g")
        np.savetxt(
                f"figdata/harmonic_current{label}.dat",
                np.array([t, *current, *phase]).T,
                header = "t " + " ".join(f"i{vac:.3g}" for vac in valid_vac_omega) + " " + " ".join(f"p{vac:.3g}" for vac in valid_vac_omega),
                comments = "",
                fmt = "%.6g")
        np.savetxt(
                f"figdata/phase_crossings{label}.dat",
                np.array([np.arange(n_phase_crossings), *phase_crossing_times, *phase_crossing_currents]).T,
                header = "phase " + " ".join(f"t{vac:.3g}" for vac in valid_vac_omega) + " " + " ".join(f"i{vac:.3g}" for vac in valid_vac_omega),
                comments = "",
                fmt = "%.6g")
    
        # adiabatic result
        n = np.arange(1, 101, 2)
        adiabatic_modes = []
    
        data = dm.list(vac=0, omega=0, include_Ga=True, d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, nmax=0, voltage_branches=4, xL=0.5, truncation_order=3, integral_method=-15)
        data = data.sort_values("vdc")
        vdc = np.concatenate((-data.vdc[:0:-1], data.vdc))
        idc = np.concatenate((data.dc_current[:0:-1], data.dc_current))
        interp = interp1d(vdc, idc, kind="cubic")
        for i, vac in enumerate(omega*vac_omega):
            ni, nierr = quad_vec((lambda t: np.sin(n*t) * interp(vac*np.sin(t))), 0, np.pi, limit=200)
            adiabatic_modes.append(np.abs(ni/ni[0]))
        np.savetxt(
                f"figdata/harmonic_modes_adiabatic{label}.dat",
                np.array([n, *adiabatic_modes]).T,
                header = "n " + " ".join(f"mode{vac:.0f}" for vac in omega*vac_omega),
                comments = "",
                fmt = "%.6g")
    
        fig, (ax1, ax2) = plt.subplots(nrows=2)
        for i, (c, p, pt, pc) in enumerate(zip(current, phase, phase_crossing_times, phase_crossing_currents)):
            color = f"C{i+1}"
            ax1.plot(t, c, color=color)
            ax1.plot(pt, pc, "x", color=color)
            ax2.plot(p, c, color=color)
            ax2.plot(np.arange(n_phase_crossings), pc, "x", color=color)
        plt.show()
    
    
    def export_convergence_pgfplots(simplified_initial_conditions = False):
        """
        Plot convergence of current as function of D (Λ₀) for a fine
        grid of dc and ac voltage at fixed omega.
        """
        dm = DataManager()
        parameters = dict(
                omega = 16.5372,
                method = "mu",
                include_Ga = True,
                voltage_branches = 4,
                integral_method = -15,
                solver_tol_rel = 1e-9,
                solver_tol_abs = 1e-11,
                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"],
                good_flags = DataManager.SOLVER_FLAGS["include_Ga"],
                    )
        if simplified_initial_conditions:
            parameters["good_flags"] ^= DataManager.SOLVER_FLAGS["simplified_initial_conditions"]
            parameters["bad_flags"] ^= DataManager.SOLVER_FLAGS["simplified_initial_conditions"]
    
        vac_num = 26
        vdc_num = 26
        vdc_max = 82.686
        vac_max = 82.686
        exponent = 2 if parameters["good_flags"] & DataManager.SOLVER_FLAGS["simplified_initial_conditions"] else 3
        d_num = 5
        log10d_min = 5
        log10d_max = 9
        d_arr = np.logspace(log10d_min, log10d_max, d_num)
        gdc = np.empty((d_num, vac_num, vdc_num), dtype=np.float64)
        idc = np.empty((d_num, vac_num, vdc_num), dtype=np.float64)
        iac = np.empty((d_num, vac_num, vdc_num), dtype=np.float64)
        phase = np.empty((d_num, vac_num, vdc_num), dtype=np.float64)
        for i, d in enumerate(d_arr):
            data = dm.list(d=d, **parameters)
            vdc, vac, gdc[i], idc[i], iac[i], phase[i] = filter_grid_data(data, vac_max=vac_max, vdc_max=vdc_max, vac_num=vac_num, vdc_num=vdc_num)
    
        gdc_fit_a = np.ndarray((vac_num, vdc_num), dtype=np.float64)
        gdc_fit_b = np.ndarray((vac_num, vdc_num), dtype=np.float64)
        idc_fit_a = np.ndarray((vac_num, vdc_num), dtype=np.float64)
        idc_fit_b = np.ndarray((vac_num, vdc_num), dtype=np.float64)
        iac_fit_a = np.ndarray((vac_num, vdc_num), dtype=np.float64)
        iac_fit_b = np.ndarray((vac_num, vdc_num), dtype=np.float64)
        j = np.vectorize(lambda d: solveTV0_scalar(d, rtol=1e-9, atol=1e-11)[2])
        logd_inv3 = j(d_arr)**exponent
        #logd_inv3 = np.log(d_arr)**-exponent
        fit_func = lambda x, y: x[0] + x[1] * logd_inv3 - y
        for i in range(vac_num):
            for j in range(vdc_num):
                (gdc_fit_a[i,j], gdc_fit_b[i,j]), trash = leastsq(fit_func, (gdc[-1,i,j], 0.), args=(gdc[:,i,j]))
                (idc_fit_a[i,j], idc_fit_b[i,j]), trash = leastsq(fit_func, (idc[-1,i,j], 0.), args=(idc[:,i,j]))
                (iac_fit_a[i,j], iac_fit_b[i,j]), trash = leastsq(fit_func, (iac[-1,i,j], 0.), args=(iac[:,i,j]))
    
        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)))
        #extent_dc = (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)))
        #extent_ac = (-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)))
        print(f"Extent: xmin={extent[0]/TK_VOLTAGE:.6g}, xmax={extent[1]/TK_VOLTAGE:,.6g}, ymin={extent[2]/TK_VOLTAGE:.6g}, ymax={extent[3]/TK_VOLTAGE:,.6g}")
        #print(f"Extent: xmin={extent_dc[0]/TK_VOLTAGE:.6g}, xmax={extent_dc[1]/TK_VOLTAGE:,.6g}, ymin={extent_dc[2]/TK_VOLTAGE:.6g}, ymax={extent_dc[3]/TK_VOLTAGE:,.6g}")
        #print(f"Extent: xmin={extent_ac[0]/TK_VOLTAGE:.6g}, xmax={extent_ac[1]/TK_VOLTAGE:,.6g}, ymin={extent_ac[2]/TK_VOLTAGE:.6g}, ymax={extent_ac[3]/TK_VOLTAGE:,.6g}")
    
        suffix = "_bad" if simplified_initial_conditions else "_good"
        def export_img(name, array):
            imwrite(f"figdata/{name}{suffix}.png", np.array(0xffff*plt.cm.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)
        export_img("idc_converged", idc_fit_a/i_max)
        export_img("iac_converged", iac_fit_a*(2/i_max))
        print(f"I range: 0, {i_max/TK_VOLTAGE:.6g}")
    
        idc_diff = (idc_fit_b*(-logd_inv3[4])/idc_fit_a)[:,1:]
        iac_diff = (iac_fit_b*(-logd_inv3[4])/iac_fit_a)[1:]
        idc_diff_max = np.nanmax(idc_diff)
        iac_diff_max = np.nanmax(iac_diff)
        i_diff_max = max(idc_diff_max, iac_diff_max)
        export_img("idc_diff", idc_diff/i_diff_max)
        export_img("iac_diff", iac_diff/i_diff_max)
        print(f"I diff range: 0, {i_diff_max:.6g}")
    
    
    def export_floquet_matrices_pgfplots():
        dm = DataManager()
        omega = 16.5371763
        vdc = 19.84461156
        vac = 29.76691734
        vb = 4
        params = dict(d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, voltage_branches=vb, omega=omega, vdc=vdc, vac=vac)
        kondos = {}
        g200_matrices = {}
        g201_matrices = {}
        z_matrices = {}
        gamma_matrices = {}
        deltaGammaL_matrices = {}
        nmax_ref = 19
        for kondo in dm.load_all(**params):
            nmax = kondo.nmax
            idx = (kondo.method, kondo.nmax > 20, kondo.padding>0)
            if nmax > nmax_ref:
                z_matrices[idx] = kondo.z[vb, nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
                gamma_matrices[idx] = kondo.gamma[vb, nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
                deltaGammaL_matrices[idx] = kondo.deltaGammaL[nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
                g200_matrices[idx] = kondo.g2[0,0,vb, nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
                g201_matrices[idx] = kondo.g2[0,1,vb, nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
            elif nmax == nmax_ref:
                z_matrices[idx] = kondo.z[vb]
                gamma_matrices[idx] = kondo.gamma[vb]
                deltaGammaL_matrices[idx] = kondo.deltaGammaL
                g200_matrices[idx] = kondo.g2[0,0,vb]
                g201_matrices[idx] = kondo.g2[0,1,vb]
            else:
                raise ValueError(f"nmax={nmax} should not be smaller than nmax_ref={nmax_ref}")
            kondos[idx] = kondo
    
        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")
    
        z_ref_j = z_matrices[("J", True, True)]
        z_ref_mu = z_matrices[("mu", True, False)]
    
        g201_ref_j = g201_matrices[("J", True, True)]
        g201_ref_mu = g201_matrices[("mu", True, False)]
    
        export_img("floquet_matrix_z_mu_small_p", z_matrices[("mu",False,True)])
        export_img("floquet_matrix_z_mu_small", z_matrices[("mu",False,False)])
        export_img("floquet_matrix_z_mu_large_p", z_matrices[("mu",True,True)])
        export_img("floquet_matrix_z_mu_large", z_matrices[("mu",True,False)])
        export_img("floquet_matrix_z_J_small_p", z_matrices[("J",False,True)])
        export_img("floquet_matrix_z_J_small", z_matrices[("J",False,False)])
        export_img("floquet_matrix_z_J_large_p", z_matrices[("J",True,True)])
        export_img("floquet_matrix_z_J_large", z_matrices[("J",True,False)])
        export_img("floquet_matrix_z_mu_diff_p", z_matrices[("mu",False,True)] - z_ref_mu)
        export_img("floquet_matrix_z_mu_diff", z_matrices[("mu",False,False)] - z_ref_mu)
        export_img("floquet_matrix_z_J_diff_p", z_matrices[("J",False,True)] - z_ref_j)
        export_img("floquet_matrix_z_J_diff", z_matrices[("J",False,False)] - z_ref_j)
    
        export_img("floquet_matrix_g201_mu_small_p", g201_matrices[("mu",False,True)])
        export_img("floquet_matrix_g201_mu_small", g201_matrices[("mu",False,False)])
        export_img("floquet_matrix_g201_mu_large_p", g201_matrices[("mu",True,True)])
        export_img("floquet_matrix_g201_mu_large", g201_matrices[("mu",True,False)])
        export_img("floquet_matrix_g201_J_small_p", g201_matrices[("J",False,True)])
        export_img("floquet_matrix_g201_J_small", g201_matrices[("J",False,False)])
        export_img("floquet_matrix_g201_J_large_p", g201_matrices[("J",True,True)])
        export_img("floquet_matrix_g201_J_large", g201_matrices[("J",True,False)])
        export_img("floquet_matrix_g201_mu_diff_p", g201_matrices[("mu",False,True)] - g201_ref_mu)
        export_img("floquet_matrix_g201_mu_diff", g201_matrices[("mu",False,False)] - g201_ref_mu)
        export_img("floquet_matrix_g201_J_diff_p", g201_matrices[("J",False,True)] - g201_ref_j)
        export_img("floquet_matrix_g201_J_diff", g201_matrices[("J",False,False)] - g201_ref_j)
    
    
    def export_Ga_max_pgfplots():
        """
        Show relevance of Ga: Plot maximum matrix element of Ga as function
        of Vdc and Vac for fixed Ω.
        """
        dm = DataManager()
        omega = 16.5372
        parameters = dict(
                omega = omega,
                d = 1e9,
                include_Ga = True,
                integral_method = -1,
                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["extrapolate_voltage"] \
                        | DataManager.SOLVER_FLAGS["deleted"],
                good_flags = DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
                        | DataManager.SOLVER_FLAGS["include_Ga"]
                    )
        vdc = []
        vac = []
        ga_max = []
        g2_max = []
        g200c = []
        ga00_max_c = []
        nmax_arr = []
        argmax = []
        for kondo in dm.load_all(**parameters):
            try:
                try:
                    nmax = kondo.nmax
                    ga = kondo.ga
                    g2 = kondo.g2
                except AttributeError:
                    continue
                vdc.append(kondo.vdc)
                vac.append(kondo.vac)
                nmax_arr.append(nmax)
                if ga.ndim == 4:
                    idx = np.abs(ga).argmax()
                    ga_max.append(np.abs(ga.flat[idx]))
                    g2_max.append(np.abs(g2).max())
                    argmax.append(idx)
                    ga00_max_c.append(np.abs(ga[0,0,nmax-4:nmax+5,nmax-4:nmax+5]).max())
                    g200c.append(np.abs(g2[0,0,nmax,nmax]))
                elif ga.ndim == 5:
                    vb = kondo.voltage_branches
                    idx = np.abs(ga[:,:,vb]).argmax()
                    ga_max.append(np.abs(ga[:,:,vb].flat[idx]))
                    g2_max.append(np.abs(g2[:,:,vb]).max())
                    argmax.append(idx)
                    ga00_max_c.append(np.abs(ga[0,0,vb,nmax-4:nmax+5,nmax-4:nmax+5]).max())
                    g200c.append(np.abs(g2[0,0,vb,nmax,nmax]))
                else:
                    raise ValueError("Invalid shape: %s"%ga.shape)
            except:
                settings.logger.exception("Error while reading data:")
    
        vdc = np.array(vdc)
        vac = np.array(vac)
        nmax = np.array(nmax_arr)
        argmax = np.array(argmax)
        ga_max = np.array(ga_max)
        g2_max = np.array(g2_max)
        g200c = np.array(g200c)
    
        vdc_max = 165.372
        vac_max = 165.372
        vdc_num = 301
        vac_num = 201
        vdc_arr = np.linspace(0, vdc_max, vdc_num)
        vac_arr = np.linspace(0, vac_max, vac_num)
        vdc_mesh, vac_mesh = np.meshgrid(vdc_arr, vac_arr)
        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)))
        grid_ga_max = griddata((vdc, vac), ga_max, (vdc_mesh, vac_mesh), method="cubic")
        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")
    
        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))
        export_img("Ga_max_abs", grid_ga_max/norm)
        print("Ga_max_abs norm:", norm)
        grid_ga_max_rel = grid_ga_max/grid_g2_max**2
        norm = np.nanmax(np.abs(grid_ga_max_rel))
        export_img("Ga_max_rel", grid_ga_max_rel/norm)
        print("Ga_max_rel norm:", norm)
    
    
    
    
    if __name__ == "__main__":
        from sys import argv
        globals()[argv[1]]()