Skip to content
Snippets Groups Projects
Select Git revision
  • ffa6a4b906c57ba3e720681a1ab0d580aba51188
  • 5.4 default protected
  • 5.5
  • dev/5.5
  • dev/5.4
  • dev/5.3_downgrade
  • feature/experimenttime_hack
  • 5.3 protected
  • _IntenSelect5.3
  • IntenSelect5.3
  • 4.27 protected
  • 4.26 protected
  • 5.0 protected
  • 4.22 protected
  • 4.21 protected
  • UE5.4-2024.1
  • UE5.4-2024.1-rc1
  • UE5.3-2023.1-rc3
  • UE5.3-2023.1-rc2
  • UE5.3-2023.1-rc
20 results

RWTHVRUtilities.h

Blame
  • final_plots.py 23.30 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 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, BSpline, griddata, interp1d
    from scipy.special import jn
    import settings
    from data_management import DataManager, KondoImport
    from kondo import gen_init_matrix
    from numbers import Number
    
    # 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, omega, vdc, vac, nmax=50, **parameters):
        """
        Return differential conductance G calculated from photon
        assisted tunneling at given parameters.
        """
        result = np.zeros(np.broadcast_shapes(np.shape(omega), np.shape(vdc), np.shape(vac)), dtype=np.float64)
        interp = interp_vac0(dm, **parameters)
        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_omega5():
        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 = False,
                solve_integral_exactly = False,
                )
        bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
                | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
                | DataManager.SOLVER_FLAGS["include_Ga"] \
                | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
                | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
                | DataManager.SOLVER_FLAGS["deleted"]
        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 = 101,
                vdc_min = 0,
                vdc_max = 165.372,
                vdc_num = 101,
                )
        with open("figdata/vdc_vac_omega16.5372.dat", "w") as file:
            file.write("vac vdc gdc idc iac")
            for i in range(101):
                file.write("\n")
                np.savetxt(file, np.array([vac[i]/omega, vdc[i]/omega, np.pi*gdc[i], idc[i], iac[i]]).T)
        # Lines at constant Vac
        for i in range(11):
            with open(f"figdata/vdc_vac{i}_omega16.5372.dat", "w") as file:
                file.write("vac vdc gdc idc iac\n")
                np.savetxt(file, np.array([vac[10*i]/omega, vdc[10*i]/omega, np.pi*gdc[10*i], idc[10*i], iac[10*i]]).T)
        # Lines at constant Vdc
        for i in range(11):
            with open(f"figdata/vac_vdc{i}_omega16.5372.dat", "w") as file:
                file.write("vac vdc gdc idc iac\n")
                np.savetxt(file, np.array([vac[:,10*i]/omega, vdc[:,10*i]/omega, np.pi*gdc[:,10*i], idc[:,10*i], iac[:,10*i]]).T)
        # 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,
            korder = 2,
            special_selection = None,
            o2_scale_x = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
            o2_scale_y = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
            o2_scale_fixed = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
            **parameters
            ):
        """
        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
            korder: order of spline interpolation
            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=5e-5, **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")
            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=korder, ky=korder)
                #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=korder, ky=korder)
                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=korder, ky=korder)
                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:
                results[f"phase_{method}_{order}{suffix}"] = griddata(
                        xy_data,
                        reduced_data.ac_current_phase,
                        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_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,
            korder = 2,
            vdc_max_J = 82.686,
            vac_max_J = 82.686,
            voltage_branches = 4,
            ):
        """
        Export interpolated data for fixed frequency omega.
        """
        parameters = dict(d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, voltage_branches=voltage_branches, xL=0.5)
        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
        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 = np.linspace(vac_min, vac_max, ac_res),
                      korder = korder,
                      special_selection = special_selection,
                      **parameters)
    
    
    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
        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())
    
    
    if __name__ == "__main__":
        from sys import argv
        globals()[argv[1]]()