diff --git a/final_plots.py b/final_plots.py
index 4a5359615832f69a67d1ccc34428ee8f8132750c..d770e7c47829d6ae69faca085102fc9408f08149 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -107,6 +107,7 @@ def photon_assisted_tunneling(dm, omega, vdc, vac, nmax=50, **parameters):
             pass
     return result
 
+
 def gen_photon_assisted_tunneling(dm, fourier_coef, omega, vdc, nmax=50, **parameters):
     """
     Return differential conductance G calculated from generalized
@@ -272,7 +273,7 @@ def export_interp(
     The parameters vdc and vac are included as 2d arrays.
     Results have the naming convention
 
-        {observable}_{order}{padding}_{method}_{suffix}
+        {observable}_{method}_{order}{padding}{suffix}
 
     observable:
         gdc   DC differential conductance, may have suffix "_ispline"
diff --git a/gen_pulse_data.py b/gen_pulse_data.py
index a639cc9f6d868a273a5378bd20643f812af79bd4..bf188e7f39900143bfec72cecb0705cccf2746e4 100644
--- a/gen_pulse_data.py
+++ b/gen_pulse_data.py
@@ -35,16 +35,18 @@ def main():
 
     # Physical parameters
     phys_group = parser.add_argument_group(title="Physical parameters")
-    phys_group.add_argument("pulse_shape", type=str, choices=["gauss", "gauss_symmetric"], default="gauss",
+    phys_group.add_argument("pulse_shape", type=str, choices=["gauss", "gauss_symmetric", "updown_gauss"], default="gauss",
             help="name for preset parameters")
     phys_group.add_argument("--omega", metavar='float', type=float, default=0.,
             help = "Frequency, units of Tkrg")
     phys_group.add_argument("--pulse_duration", metavar='float', type=float,
-            help = "duration of voltage pulses")
+            help = "duration of voltage pulses, units of 1/Tkrg")
     phys_group.add_argument("--pulse_phase", metavar='float', type=float,
             help = "phase/2π accumulated during of voltage pulses")
     phys_group.add_argument("--pulse_height", metavar='float', type=float,
             help = "height of voltage pulses")
+    phys_group.add_argument("--pulse_separation", metavar='float', type=float,
+            help = "separation of two voltage pulses, units of 1/Tkrg")
     phys_group.add_argument("--baseline_voltage", metavar='float', type=float, default=0.,
             help="baseline voltage, units of Tkrg")
     phys_group.add_argument("--xL", metavar='float', type=float, default=0.5,
@@ -162,26 +164,39 @@ def main():
             )
 
     pulse_shape = options.pop("pulse_shape")
-    if pulse_shape == "gauss":
-        vdc, fourier_coef = fourier_coef_gauss(
-                omega = options["omega"],
-                baseline = options.pop("baseline_voltage"),
-                duration = options.pop("pulse_duration"),
-                height = options.pop("pulse_height"),
-                phase = options.pop("pulse_phase"),
-                nmax = options["nmax"],
-                )
-    elif pulse_shape == "gauss_symmetric":
-        vdc, fourier_coef = fourier_coef_gauss_symmetric(
-                omega = options["omega"],
-                baseline = options.pop("baseline_voltage"),
-                duration = options.pop("pulse_duration"),
-                height = options.pop("pulse_height"),
-                phase = options.pop("pulse_phase"),
-                nmax = options["nmax"],
-                )
-    else:
-        raise ValueError("invalid pulse shape")
+    match pulse_shape:
+        case "gauss":
+            vdc, fourier_coef = fourier_coef_gauss(
+                    omega = options["omega"],
+                    baseline = options.pop("baseline_voltage"),
+                    duration = options.pop("pulse_duration"),
+                    height = options.pop("pulse_height"),
+                    phase = options.pop("pulse_phase"),
+                    nmax = options["nmax"],
+                    )
+            options.pop("pulse_separation")
+        case "gauss_symmetric":
+            vdc, fourier_coef = fourier_coef_gauss_symmetric(
+                    omega = options["omega"],
+                    baseline = options.pop("baseline_voltage"),
+                    duration = options.pop("pulse_duration"),
+                    height = options.pop("pulse_height"),
+                    phase = options.pop("pulse_phase"),
+                    nmax = options["nmax"],
+                    )
+            options.pop("pulse_separation")
+        case "updown_gauss":
+            vdc, fourier_coef = fourier_coef_updown_gauss(
+                    omega = options["omega"],
+                    baseline = options.pop("baseline_voltage"),
+                    duration = options.pop("pulse_duration"),
+                    height = options.pop("pulse_height"),
+                    phase = options.pop("pulse_phase"),
+                    separation = options.pop("pulse_separation"),
+                    nmax = options["nmax"],
+                    )
+        case _:
+            raise ValueError("invalid pulse shape")
 
     # Generate data
     dm = DataManager()
@@ -260,5 +275,21 @@ def fourier_coef_gauss(
     return fourier_coef[0].real, fourier_coef[1:]
 
 
+def fourier_coef_updown_gauss(
+        nmax,
+        omega,
+        separation,
+        duration = None,
+        height = None,
+        phase = None,
+        baseline = 0.,
+        ):
+    assert baseline == 0.
+    narr = np.arange(1, nmax+1)
+    trash, coef_gauss = fourier_coef_gauss(nmax, omega, duration, height, phase, 0.)
+    coef = 2j*np.sin(0.5*narr*separation*omega)*coef_gauss
+    return 0., coef
+
+
 if __name__ == "__main__":
     main()
diff --git a/plot.py b/plot.py
index 6862d67ab587943de8f47ad7a3fc430700a22277..e791d52565ccab6f3956d600b0a8d93b9b463b07 100644
--- a/plot.py
+++ b/plot.py
@@ -2175,5 +2175,49 @@ def compare_limits(
     plt.show()
 
 
+def overview_3d_contours(dm, show_contours=False, **trash):
+    import mpl_toolkits.mplot3d
+    data_vdc0 = np.load("figdata/vdc0_interp.npz")
+    data_omega5 = np.load("figdata/omega5_interp.npz")
+    interp = interp_vac0(dm, method="mu", d=1e9, include_Ga=True, voltage_branches=4, solver_tol_rel=1e-8, solver_tol_abs=1e-10, integral_method=-15)
+    fig = plt.figure()
+    ax = fig.add_subplot(projection="3d")
+    fill_levels = 1.1/np.linspace(1, 6.12, 256)[::-1] - 0.1
+    colors = plt.cm.viridis(np.linspace(0, 1, 256))
+    levels = 1/np.arange(1, 50, 0.5)[::-1]
+    vdc_resolution = 400
+    omega_resolution = 200
+    omega_max = 16.5372
+
+    vdc0_selection = data_vdc0["omega"][0] <= omega_max
+    vdc0_y = data_vdc0["vac_omega"][:,vdc0_selection]
+    vdc0_x = np.zeros_like(vdc0_y)
+    vdc0_z = data_vdc0["omega"][:,vdc0_selection]
+    vdc0_m = np.pi*data_vdc0["gdc_J_o3a_p"][:,vdc0_selection]
+
+    ax.contourf(vdc0_m, vdc0_y, vdc0_z, zdir="x", levels=fill_levels, offset=0, zorder=-10, colors=colors)
+    if show_contours:
+        ax.contour(vdc0_m, vdc0_y, vdc0_z, zdir="x", levels=levels, offset=-0.1, zorder=10, colors="black")
+
+    omega5 = data_omega5["omega"]
+    omega5_x = data_omega5["vdc"]/omega5
+    omega5_y = data_omega5["vac"]/omega5
+    omega5_m = np.pi*data_omega5["gdc_mu_o3a"]
+    ax.contourf(omega5_x, omega5_y, omega5_m, zdir="z", levels=fill_levels, offset=omega5, colors=colors)
+    if show_contours:
+        ax.contour(omega5_x, omega5_y, omega5_m, zdir="z", levels=levels, offset=omega5+0.1, colors="black")
+
+    vac0_x = np.linspace(0, 10, vdc_resolution)
+    vac0_z = np.linspace(1e-3, omega_max, omega_resolution)
+    vac0_vdc = vac0_x.reshape((-1,1)) * vac0_z.reshape((1,-1))
+    vac0_m = np.pi*interp(vac0_vdc)
+    ax.contourf(vac0_x, vac0_m, vac0_z, zdir="y", levels=fill_levels, offset=0, zorder=-10, colors=colors)
+    if show_contours:
+        ax.contour(vac0_x, vac0_m, vac0_z, zdir="y", levels=levels, offset=-0.1, zorder=10, colors="black")
+
+    #ax.plot_surface(vdc0_x, vdc0_y, vdc0_z, rstride=1, cstride=1, facecolors=plt.cm.viridis(vdc0_m), shade=False, zorder=-10)
+    plt.show()
+
+
 if __name__ == '__main__':
     main()
diff --git a/plot_pulses.py b/plot_pulses.py
index 06c48abbd00725a19aeadf710a710ea1a67042f9..df8e730980ebf91cc56765f65c796dc93e7294eb 100644
--- a/plot_pulses.py
+++ b/plot_pulses.py
@@ -5,10 +5,12 @@
 """
 Kondo FRTRG, plot data for pulsed voltage
 """
+import os
 import argparse
 import numpy as np
 import settings
-from data_management import DataManager
+import tables as tb
+from data_management import DataManager, KondoImport
 from scipy.optimize import leastsq
 from gen_pulse_data import fourier_coef_gauss, fourier_coef_gauss_symmetric
 import matplotlib.pyplot as plt
@@ -25,7 +27,7 @@ def plot_gammaL(dm, omega=0., pulse_duration=0.1, pulse_height=4., pulse_phase=N
         else:
             settings.logger.warn(f"unknown pulse shape: {shape}")
             continue
-        for kondo in dm.load_all(vdc=vdc, has_fourier_coef=True, fourier_coef=fourier_coef, **parameters):
+        for kondo in dm.load_all(vdc=vdc, omega=omega, has_fourier_coef=True, fourier_coef=fourier_coef, **parameters):
             nmax = kondo.nmax
             gammaL = kondo.gammaL
             img = plt.imshow(np.abs(gammaL), norm=LogNorm(1e-6))
@@ -33,16 +35,42 @@ def plot_gammaL(dm, omega=0., pulse_duration=0.1, pulse_height=4., pulse_phase=N
             plt.show()
 
 
-def current_time(dm, omega=0., pulse_duration=0.1, pulse_height=4., pulse_phase=None, baseline=0., pulse_shape=["gauss", "gauss_symmetric"], resolution=500, **parameters):
+def current_time_all(dm, pulse_duration=0.1, pulse_phase=0.5, pulse_height=None, baseline=0., resolution=500, **parameters):
+    x = np.linspace(-np.pi, np.pi, resolution)
+    fig, (ax_voltage, ax_current) = plt.subplots(nrows=2, ncols=1, sharex=True)
+    for (row, kondo) in load_all_pulses_full(dm, pulse_duration=pulse_duration, pulse_phase=pulse_phase, pulse_height=pulse_height):
+            nmax = kondo.nmax
+            gammaL = kondo.gammaL
+            if gammaL.ndim == 2:
+                current_coef = gammaL[:,nmax]
+            elif gammaL.ndim == 1:
+                current_coef = gammaL
+            assert current_coef.size == 2*nmax+1
+            current = np.zeros(resolution, dtype=np.complex128)
+            for n, c in enumerate(current_coef, -nmax):
+                current += c*np.exp(-1j*n*x)
+            #print(kondo.omega, (0.5*kondo.dc_current + integrate_ft(-0.1, 0.4, current_coef[nmax-1::-1]))/kondo.omega*2*np.pi, kondo.dc_current/kondo.omega*2*np.pi)
+            ax_current.plot(x/kondo.omega, current.real)
+            ax_current.plot(x/kondo.omega, current.imag)
+            voltage = np.ones(resolution, dtype=np.float64) * kondo.vdc
+            for n, c in enumerate(kondo.fourier_coef, 1):
+                voltage += 2*(c*np.exp(1j*n*x)).real
+            ax_voltage.plot(x/kondo.omega, voltage)
+    plt.show()
+
+
+
+def current_time(dm, omega=0., pulse_duration=0.1, pulse_phase=None, pulse_height=4., baseline=0., pulse_shape=["gauss", "gauss_symmetric"], resolution=500, **parameters):
     x = np.linspace(-np.pi, np.pi, resolution)
     for shape in pulse_shape:
-        if shape == "gauss":
-            vdc, fourier_coef = fourier_coef_gauss(10, omega, pulse_duration, pulse_height, pulse_phase, baseline)
-        elif shape == "gauss_symmetric":
-            vdc, fourier_coef = fourier_coef_gauss_symmetric(10, omega, pulse_duration, pulse_height, pulse_phase, baseline)
-        else:
-            settings.logger.warn(f"unknown pulse shape: {shape}")
-            continue
+        match shape:
+            case "gauss":
+                vdc, fourier_coef = fourier_coef_gauss(10, omega, pulse_duration, pulse_height, pulse_phase, baseline)
+            case "gauss_symmetric":
+                vdc, fourier_coef = fourier_coef_gauss_symmetric(10, omega, pulse_duration, pulse_height, pulse_phase, baseline)
+            case _:
+                settings.logger.warn(f"unknown pulse shape: {shape}")
+                continue
         for kondo in dm.load_all(vdc=vdc, has_fourier_coef=True, fourier_coef=fourier_coef, **parameters):
             nmax = kondo.nmax
             gammaL = kondo.gammaL
@@ -63,44 +91,128 @@ def current_time(dm, omega=0., pulse_duration=0.1, pulse_height=4., pulse_phase=
     plt.show()
 
 
-def fit_fourier_coef(fourier_coef):
+def integrate_ft(t0, t1, fourier_coef):
+    """
+    t1                             i2πnt
+    ∫ f(t) dt  where f(t) =  Σ fn e
+    t0                       n
+
+    t0 and t1 should be given in units of T=2π/Ω.
+    fourier_coef should be given in units of Ω and start with the first
+    (not zeroth) coefficient.
+    """
+    fourier_coef = np.asarray(fourier_coef)
+    phase_arr = 2*np.pi*np.arange(1, fourier_coef.size+1)
+    return ((fourier_coef * (np.exp(1j*t1*phase_arr) - np.exp(1j*t0*phase_arr))).imag / phase_arr).sum() * 2
+
+
+def plot_transported_charge(dm, **parameters):
+    table = load_all_pulses(dm)
+    table = table[table.pulse_type == "gauss"]
+    fig, ax = plt.subplots()
+    ax.scatter(table.pulse_duration, table.dc_current/table.omega*2*np.pi, c=table.pulse_phase, s=100/table.omega**2, marker="o")
+
+    sym_pulse_omega = []
+    sym_pulse_duration = []
+    sym_pulse_phase = []
+    sym_pulse_charge = []
+    _h5files = set()
+    for (row, kondo) in load_all_pulses_full(dm, pulse_type="gauss_symmetric"):
+        ref_time = min(3*row.pulse_duration*row.omega/(2*np.pi), 0.15)
+        nmax = kondo.nmax
+        sym_pulse_charge.append(integrate_ft(-ref_time, 0.5-ref_time, kondo.gammaL[nmax-1::-1,nmax])/row.omega*2*np.pi)
+        _h5files.add(kondo._h5file)
+        sym_pulse_omega.append(row.omega)
+        sym_pulse_duration.append(row.pulse_duration)
+        sym_pulse_phase.append(row.pulse_phase)
+    for file in _h5files:
+        file.close()
+
+    ax.scatter(sym_pulse_duration, sym_pulse_charge, c=sym_pulse_phase, marker="^", s=80/np.array(sym_pulse_omega)**2)
+    plt.show()
+
+
+def fit_fourier_coef(fourier_coef, tol=1e-2):
     # gauss
     narr = np.arange(1, fourier_coef.size+1)
-    func = lambda param: (param[0]*np.exp(-param[1]*narr**2) - fourier_coef).view(dtype=np.float64)
-    fit, cov = leastsq(func, (0.5, 0.1))
-    residuals = (func(fit)**2).sum()
-    settings.logger.debug(f"fit_fourier_coef: fit {fit}, residuals {residuals:.3g}")
-    if residuals < 1e-3:
-        return dict(type="gauss", prefactor=fit[0], smoothen=fit[1])
+    fit_functions = dict(
+            gauss = lambda param: (param[0]*np.exp(-param[1]*narr**2) - fourier_coef).view(dtype=np.float64),
+            gauss_symmetric = lambda param: (param[0]*np.exp(-param[1]*((narr+1)/2)**2)*(narr%2) - fourier_coef).view(dtype=np.float64),
+            )
+    for name, func in fit_functions.items():
+        fit, cov = leastsq(func, (0.2, 0.1))
+        residuals = (func(fit)**2).sum()
+        settings.logger.debug(f"fit_fourier_coef: fit {name} yields {fit}, residuals {residuals:.3g}")
+        if residuals < tol:
+            return dict(type=name, prefactor=fit[0], smoothen=fit[1])
     return dict(type="unknown")
 
 
-def parameter_overview(dm, **trash):
-    table = dm.list_fourier_coef()
+def load_all_pulses(dm, **parameters):
+    table = dm.list_fourier_coef(**parameters)
     table["pulse_type"] = "unknown"
-    table["pulse_param1"] = 0.
-    table["pulse_param2"] = 0.
+    table["pulse_duration"] = 0.
+    table["pulse_phase"] = 0.
+    table["pulse_height"] = 0.
     fourier_coef = np.array([table[f"fcr{i}"]+1j*table[f"fci{i}"] for i in range(10)]).T
     for coef, idx in zip(fourier_coef, table.index):
         parameters = fit_fourier_coef(coef)
         match parameters["type"]:
             case "gauss":
-                table.pulse_type[idx] = "gauss"
-                table.pulse_param1[idx] = parameters["prefactor"]
-                table.pulse_param2[idx] = parameters["smoothen"]
+                table.loc[idx, "pulse_type"] = "gauss"
+                table.loc[idx, "pulse_duration"] = 4*np.sqrt(parameters["smoothen"]*np.log(2)) / table.omega[idx]
+                table.loc[idx, "pulse_height"] = parameters["prefactor"] / np.sqrt(parameters["smoothen"]) * np.pi**0.5
+                table.loc[idx, "pulse_phase"] = parameters["prefactor"] / table.omega[idx]
+            case "gauss_symmetric":
+                table.loc[idx, "pulse_type"] = "gauss_symmetric"
+                table.loc[idx, "pulse_duration"] = 2*np.sqrt(parameters["smoothen"]*np.log(2)) / table.omega[idx]
+                table.loc[idx, "pulse_height"] = parameters["prefactor"] / np.sqrt(parameters["smoothen"]) * np.pi**0.5
+                table.loc[idx, "pulse_phase"] = parameters["prefactor"] / (2*table.omega[idx])
             case _:
-                print(f"unknown type: {parameters['type']}")
+                settings.logger.debug(f"unknown type: {parameters['type']}")
+    return table
+
+
+def load_all_pulses_full(dm, pulse_type="all", pulse_duration=None, pulse_phase=None, pulse_height=None, **parameters):
+    table = load_all_pulses(dm, **parameters)
+    if pulse_type == "all":
+        sel = table.pulse_type != "unknown"
+    else:
+        sel = table.pulse_type == pulse_type
+    if pulse_duration is not None:
+        sel &= np.abs(table.pulse_duration - pulse_duration) < 1e-3 * pulse_phase
+    if pulse_phase is not None:
+        sel &= np.abs(table.pulse_phase - pulse_phase) < 1e-3
+    if pulse_height is not None:
+        sel &= np.abs(table.pulse_height - pulse_height) < 1e-3 * pulse_height
+    for ((dirname, basename), subtable) in table[sel].groupby(["dirname", "basename"]):
+        try:
+            h5file = tb.open_file(os.path.join(dirname, basename))
+        except:
+            settings.logger.exception("Error while loading HDF5 file")
+            continue
+        #metadatatable = h5file.get_node('/metadata/mdtable')
+        for index, row in subtable.iterrows():
+            try:
+                datanode = h5file.get_node("/data/" + row.hash)
+                #metadatarow = metadatatable.where("hash == '%s'"%(row.hash))
+                yield row, KondoImport(row, datanode, h5file, owns_h5file=subtable.shape[0]==1)
+            except:
+                settings.logger.exception("Error while loading data")
+
+
+def parameter_overview(dm, **trash):
+    table = load_all_pulses(dm)
     fig, ax_gauss = plt.subplots()
 
-    sel = table.pulse_type == "gauss"
-    gauss_frequency = table.omega[sel]
-    gauss_prefactor = table.pulse_param1[sel]
-    gauss_smoothen = table.pulse_param2[sel]
-    gauss_duration = 4*np.sqrt(gauss_smoothen*np.log(2)) / gauss_frequency
-    gauss_height = gauss_prefactor / np.sqrt(gauss_smoothen) * np.pi**0.5
-    gauss_phase = gauss_height * gauss_duration / (4*(np.log(2)*np.pi)**0.5)
-    gauss_phase = gauss_prefactor / gauss_frequency
-    ax_gauss.plot(gauss_duration, gauss_phase, "o")
+    sel = (table.pulse_type == "gauss") | (table.pulse_type == "gauss_symmetric")
+    gauss_frequency = table.omega[sel] * (1 + (table.pulse_type[sel] == "gauss_symmetric"))
+    gauss_duration = table.pulse_duration[sel]
+    gauss_height = table.pulse_height[sel]
+    gauss_phase = table.pulse_phase[sel]
+    gauss_good = table.solver_flags[sel] & (DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["second_order_rg_equations"]) == DataManager.SOLVER_FLAGS["include_Ga"]
+    ax_gauss.plot(gauss_duration[~gauss_good], gauss_phase[~gauss_good], "x")
+    ax_gauss.plot(gauss_duration[gauss_good], gauss_phase[gauss_good], "o")
     plt.show()
 
 
diff --git a/plot_pyqtgraph.py b/plot_pyqtgraph.py
index e81fde6d469b8849616ed6427003fc051932e374..0e6b92104d623f9d930b4c118c7bcd260cee442d 100644
--- a/plot_pyqtgraph.py
+++ b/plot_pyqtgraph.py
@@ -33,6 +33,7 @@ def full_overview(dm,
         plot_value = "G",
         cmap = "viridis",
         show_gpat = False,
+        units = "tkrg",
         **parameters):
     """
     Show all data for conductance.
@@ -76,10 +77,18 @@ def full_overview(dm,
         values = np.broadcast_to(vac0_data.gamma.to_numpy()/20, (2, vac0_data.shape[0])).T
     z = np.zeros((vac0_data.shape[0], 2))
     z[:,1] = 25
-    su = gl.GLSurfacePlotItem(x=vac0_data.vdc.to_numpy(), y=np.array([0,0]), z=z, colors=cmap(values))
-    w.addItem(su)
+
     # overview
-    pos = np.array([data.vdc, data.vac, data.omega]).T
+    match units:
+        case "omega":
+            pos = np.array([data.vdc/data.omega, data.vac/data.omega, data.omega]).T
+        case "tkrg":
+            su = gl.GLSurfacePlotItem(x=vac0_data.vdc.to_numpy(), y=np.array([0,0]), z=z, colors=cmap(values))
+            w.addItem(su)
+            pos = np.array([data.vdc, data.vac, data.omega]).T
+        case _:
+            raise ValueError(f"Invalid value for units: {units}. Possible values are omega and tkrg.")
+
     if xyscale == "log":
         pos[:,0] = np.log10(pos[:,0])
         pos[:,1] = np.log10(pos[:,1])
@@ -192,6 +201,27 @@ def fixed_parameter(dm,
     app.exec_()
 
 
+def stacked_frequencies(dm, omegas=(7.1271, 16.5372), shift=2, scale=1, size=0.1, xyscale=None, zscale=None, grid=None, cmap="viridis", gl_preset=None, **parameters):
+    parameters.pop("omega", None)
+    app = QtGui.QApplication([])
+    w = gl.GLViewWidget()
+    w.show()
+    w.setWindowTitle("Kondo model: stacked frequencies")
+    plots = []
+    cmap = cm.get_cmap(cmap)
+    for n, omega in enumerate(omegas):
+        data = dm.list(omega=omega, **parameters)
+        pos = np.array([
+                data.vdc/omega,
+                data.vac/omega,
+                np.log(data.dc_conductance*np.pi)*scale + n*shift
+            ]).T
+        sp = gl.GLScatterPlotItem(pos=pos, size=size, color=cmap(data.dc_conductance*np.pi), pxMode=False)
+        w.addItem(sp)
+        plots.append(sp)
+    app.exec_()
+
+
 def vdc0_scaled(dm,
         scale = 80,
         size = None,
@@ -221,7 +251,7 @@ def vdc0_scaled(dm,
         pos[:,0] = np.log10(pos[:,0])
         pos[:,1] = np.log10(pos[:,1])
     if zscale == "log":
-        pos[:,2] = np.log10(pos[:,2])
+        pos[:,2] = scale*np.log10(pos[:,2]/scale)
     if size is None:
         size = scale**0.5/20
     sp = gl.GLScatterPlotItem(
@@ -681,6 +711,7 @@ PRESETS = dict(
         iac = dict(d=1e9, truncation_order=3, plot_value="ac_current_abs", scale=20, size=1),
         phase = dict(d=1e9, truncation_order=3, plot_value="ac_current_phase", scale=100, size=1),
         compare_RGeq = dict(function=compare_RGeq),
+        stacked = dict(function=stacked_frequencies),
         )
 
 def main(dm, preset=None, **parameters):