diff --git a/experiments.py b/experiments.py
index 0ac3fd0ec76cd6193c37fd372ee67531cf74cf60..9df517b0ddcf7ba58dba5558b6140fb5843fcc26 100644
--- a/experiments.py
+++ b/experiments.py
@@ -274,31 +274,33 @@ def kogan04_calibrate(dm, **kwargs):
     #omega = round(omega, 4)
     colors = ("red", "green", "blue", "orange", "violet")
     #omega = 7.1271
-    omega = 3.5
-    fig, ax = adjust_background(
+    #omega = 3.5
+    omega = 3.62
+    fig1, ax1 = adjust_background(
             dm,
             omega = omega,
-            vac_arr = np.array([0.72881, 1.13091, 1.50788, 1.6838, 3.61891])*omega,
+            #vac_arr = np.array([0.72881, 1.13091, 1.50788, 1.6838, 3.61891])*omega,
+            vac_arr = np.array([0.73922, 1.14706, 1.52942, 1.70785, 3.6706])*omega,
             xscale = 1/(omega * sc.eV/(13.47e12*sc.h)),
             colors = colors,
             include_Ga = True,
             integral_method = -15,
             **kwargs,
             )
-    #fig, ax = plot_calibration_lines(
-    #        dm,
-    #        omega = omega,
-    #        #omega = 7.418,
-    #        vac = np.array([29, 45, 60, 67, 144]) * omega * sc.eV/(13.47e15*sc.h),
-    #        colors = colors,
-    #        xscale = 1/(omega * sc.eV/(13.47e12*sc.h)),
-    #        yscale = np.pi,
-    #        calibrate_min = 1,
-    #        calibrate_max = 1.6,
-    #        include_Ga = True,
-    #        integral_method = -15,
-    #        **kwargs,
-    #        )
+    fig2, ax2 = plot_calibration_lines(
+            dm,
+            omega = omega,
+            #omega = 7.418,
+            vac = np.array([29, 45, 60, 67, 144]) * omega * sc.eV/(13.47e15*sc.h),
+            colors = colors,
+            xscale = 1/(omega * sc.eV/(13.47e12*sc.h)),
+            yscale = np.pi,
+            calibrate_min = 1,
+            calibrate_max = 1.6,
+            include_Ga = True,
+            integral_method = -15,
+            **kwargs,
+            )
 
     shift_Vdc = 0.017
     for i, trace in enumerate((0, 10, 21, 26, 40)):
@@ -306,7 +308,8 @@ def kogan04_calibrate(dm, **kwargs):
         vdc_exp = data[:,1] # in mV
         g_exp = data[:,0] / 2 # in 2e²/h = 1/π
         # Plot original data
-        ax.plot(vdc_exp, g_exp, '.', markersize=1, color=colors[i])
+        ax1.plot(vdc_exp, g_exp, '.', markersize=1, color=colors[i])
+        ax2.plot(vdc_exp, g_exp, '.', markersize=1, color=colors[i])
         # Plot symmetrized data
         sorting = np.argsort(vdc_exp)
         vdc_exp = vdc_exp[sorting]
@@ -316,9 +319,10 @@ def kogan04_calibrate(dm, **kwargs):
         assert truncate > 0
         g_mean = (g_exp[truncate:] + g_exp[:truncate-1:-1]) / 2
         spline = BSpline(*splrep(vdc_exp[truncate:], g_mean, k=3, s=2e-4 if trace==40 else 6e-5), extrapolate=False)
-        ax.plot(vdc_exp[truncate:], spline(vdc_exp[truncate:]), color=colors[i], alpha=0.5, zorder=10)
+        ax1.plot(vdc_exp[truncate:], spline(vdc_exp[truncate:]), color=colors[i], alpha=0.5, zorder=10)
+        ax2.plot(vdc_exp[truncate:], spline(vdc_exp[truncate:]), color=colors[i], alpha=0.5, zorder=10)
 
-    return fig, ax
+    return (fig1, fig2), (ax1, ax2)
 
 def bruhat18_calibrate(dm, frequency_GHz=12, **kwargs):
     """
@@ -640,7 +644,8 @@ def kogan04(dm, **parameters):
     #omega = round(omega, 5)
     #omega = 7.1271
     #omega = 7
-    omega = 3.5
+    omega = 3.62
+    correction = 1.42
     vac_mueV_arr = np.array([29, 45, 60, 67, 144])
     vac_omega_arr = vac_mueV_arr / omega_mueV
     vac_arr = vac_mueV_arr / tkrg_mueV
@@ -648,8 +653,8 @@ def kogan04(dm, **parameters):
     print(f"omega = {omega}")
     print(f"Vac = {vac_arr}")
     print(f"Vac/omega = {vac_omega_arr}")
-    print(f"corrected Vac = {1.4*vac_arr}")
-    print(f"corrected Vac = {(1.4*vac_omega_arr).round(5)} Ω")
+    print(f"corrected Vac = {correction*vac_arr}")
+    print(f"corrected Vac = {(correction*vac_omega_arr).round(5)} Ω")
     vdc_mueV_max = 400
     vdc_max = vdc_mueV_max / omega_mueV * omega
     print(f"Vdc max. = {vdc_max}")
diff --git a/final_plots.py b/final_plots.py
index ccb8c0c540741b9c1ffee7b10f6c63a71c47d7d2..2f1fba2a0171344e08954fdb7e277af81d552850 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -29,7 +29,7 @@ from plot_pulses import load_all_pulses_full, integrate_ft
 from kondo import solveTV0_scalar
 
 # color maps with better sampling
-viridis = plt.cm.viridis.resampled(0x10000)
+viridis = mplcolors.LinearSegmentedColormap.from_list('viridis', plt.cm.viridis.colors, 0x10000)
 seismic = plt.cm.seismic.resampled(0x10000)
 
 # In this program all energies are given in units of the RTRG Kondo
@@ -235,6 +235,7 @@ def export_interp(
         o3r_scale_x = TK_VOLTAGE_O3/TK_VOLTAGE,
         o3r_scale_y = TK_VOLTAGE_O3/TK_VOLTAGE,
         o3r_scale_fixed = TK_VOLTAGE_O3/TK_VOLTAGE,
+        g_flags_bad = DataManager.SOLVER_FLAGS["improved_initial_conditions"],
         **parameters
         ):
     """
@@ -319,7 +320,7 @@ def export_interp(
             o3p = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
             )
     i_flags_good = DataManager.SOLVER_FLAGS["improved_initial_conditions"]
-    g_flags_bad = DataManager.SOLVER_FLAGS["improved_initial_conditions"]
+    #g_flags_bad = DataManager.SOLVER_FLAGS["improved_initial_conditions"]
 
     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"]]
@@ -905,11 +906,11 @@ def export_omega5_interp_deviation(
 def export_vdc0_interp(
         filename = "figdata/vdc0_interp.npz",
         omega_min = 0.1,
-        omega_max = 40,
-        vac_omega_min = 0.01,
+        omega_max = 16.5372,
+        vac_omega_min = 0,
         vac_omega_max = 10,
-        omega_res = 201,
-        vac_res = 201,
+        omega_res = 401,
+        vac_res = 301,
         korder = 2,
         ):
     """
@@ -922,10 +923,12 @@ def export_vdc0_interp(
                   x_parameter = "omega",
                   y_parameter = "vac_omega",
                   y_sort_parameter = "vac",
+                  g_flags_bad = 0,
                   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,
+                  kxorder = korder,
+                  kyorder = korder,
                   special_selection = (lambda data, *args, **kwargs: data["omega"] > 0),
                   **parameters)
 
@@ -1325,13 +1328,17 @@ def export_kogan04_pgfplots():
     """
     #omega = 3.1
     #xL = 0.05
-    omega = 3.5
-    xL = 0.166667
+    #omega = 3.5
+    #xL = 0.166667
+    omega = 3.62
+    xL = 0.2
     #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])
+    #vac_omega = np.array([0.72881, 1.13091, 1.50788, 1.6838,  3.61891])
+    # adjusted by factor 1.42:
+    vac_omega = np.array([0.73922, 1.14706, 1.52942, 1.70785, 3.6706])
     parameters = dict(
             d = 1e9,
             method = "mu",
@@ -1357,7 +1364,7 @@ def export_kogan04_pgfplots():
     #g_shift_sym = 0.015
     # for omega=3.5, xL=0.166667:
     prefactor_sym = 0.096*np.pi
-    g_shift_sym = 0.02
+    g_shift_sym = 0.021
     # for omega=4, xL=0.2:
     #prefactor_sym = 0.093*np.pi
     #g_shift_sym = 0.024
@@ -1376,7 +1383,7 @@ def export_kogan04_pgfplots():
             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, *arrays = adjust_argument_array(0, 26.5, *gdc_funcs, *gdc_funcs_asym, s1=5e-4, maxit=4, initial_resolution=54)
     vdc = np.append(-vdc[:0:-1], vdc)
 
     interp = interp_vac0(dm, **parameters)
@@ -1392,7 +1399,7 @@ def export_kogan04_pgfplots():
 
 def export_bruhat18_pgfplots(
         plot_s = 2e-3,
-        init_res = 41,
+        init_res = 101,
         ):
     dm = DataManager()
     omega1 = 4.14
@@ -1447,18 +1454,20 @@ def export_bruhat18_pgfplots(
     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, *arrays1 = adjust_argument_array(0, 21.5, *gdc_funcs1, s1=plot_s, maxit=3, initial_resolution=init_res)
+    vdc1_pat, *arrays1_pat = adjust_argument_array(0, 21.5, *g1_pat_funcs, s1=1.2*plot_s, maxit=3, 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)
+    vdc1_pat = np.append(-vdc1_pat[:0:-1], vdc1_pat) / omega1
+    g1 = np.array([np.append(a[:0:-1], a) for a in arrays1])
+    g1_pat = np.array([np.append(a[:0:-1], a) for a in arrays1_pat])
+    vdc2, *arrays2 = adjust_argument_array(0, 21.5, *gdc_funcs2, s1=plot_s, maxit=4, initial_resolution=init_res)
+    vdc2_pat, *arrays2_pat = adjust_argument_array(0, 21.5, *g2_pat_funcs, s1=1.2*plot_s, maxit=3, 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:]])
+    vdc2_pat = np.append(-vdc2_pat[:0:-1], vdc2_pat) / omega2
+    g2 = np.array([np.append(a[:0:-1], a) for a in arrays2])
+    g2_pat = np.array([np.append(a[:0:-1], a) for a in arrays2_pat])
 
-    vdc1_asym, *g1_asym = adjust_argument_array(0, 21.5, *gdc_funcs1_asym, s1=plot_s, maxit=4, initial_resolution=init_res)
+    vdc1_asym, *g1_asym = adjust_argument_array(0, 21.5, *gdc_funcs1_asym, s1=plot_s, maxit=3, 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])
 
@@ -1471,27 +1480,43 @@ def export_bruhat18_pgfplots(
     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])
+    bg2_pat = np.array([quad_vec((lambda t: splev(vdc2_pat+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])
+    bg1_pat = np.array([quad_vec((lambda t: splev((vdc1_pat+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,
+               np.array([vdc1, *g1, *bg1, *(g1 + 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),
+                   )),
+               fmt = "%.6g",
+               comments = "")
+    np.savetxt("figdata/bruhat18_19GHz_pat.dat",
+               np.array([vdc1_pat, *g1_pat, *(g1_pat + bg1_pat)]).T,
+               header = " ".join((
+                   "vdc",
                    " ".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,
+               np.array([vdc2, *g2, *bg2, *(g2 + 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),
+                   )),
+               fmt = "%.6g",
+               comments = "")
+    np.savetxt("figdata/bruhat18_12GHz_pat.dat",
+               np.array([vdc2_pat, *g2_pat, *(g2_pat + bg2_pat)]).T,
+               header = " ".join((
+                   "vdc",
                    " ".join(f"gpat{v:.3g}" for v in vac_omega2),
                    " ".join(f"gpatbg{v:.3g}" for v in vac_omega2),
                    )),
diff --git a/plot.py b/plot.py
index 07529c85c5e831f11ac3b04e63d03ec95bee9e62..af28de0ef82a738273ce4320cdb9003c750b1d2d 100644
--- a/plot.py
+++ b/plot.py
@@ -7,7 +7,7 @@ Kondo FRTRG, generate plots based on save data
 """
 
 import matplotlib.pyplot as plt
-from matplotlib.colors import LogNorm, FuncNorm
+from matplotlib.colors import LogNorm, FuncNorm, LinearSegmentedColormap
 import argparse
 import pandas as pd
 import numpy as np
@@ -25,6 +25,9 @@ from data_management import DataManager, KondoImport
 from kondo import solveTV0_scalar
 from scipy.integrate import quad, quad_vec
 
+# color maps with better sampling
+viridis = LinearSegmentedColormap.from_list('viridis', plt.cm.viridis.colors, 0x10000)
+
 # reference_values maps (omega, vdc, vac) to reference values
 REFERENCE_VALUES = {
         (10, 24, 22) : dict(
@@ -2348,10 +2351,25 @@ def overview_3d_contours(dm, show_contours=False, show_surfaces=True, save_image
     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]
-    real_cmap = lambda x: plt.cm.viridis(1 - (1.1/(x+0.1)-1)/5.12)[...,:3]
-    vdc_resolution = 400
-    omega_resolution = 400
+    real_cmap = lambda x: viridis(1 - (1.1/(x+0.1)-1)/5.12)
     omega_max = 16.5372
+    vdc = data_omega5["vdc"][0]
+    #vdc_resolution = data_omega5["vdc"].shape[1]
+    omega = data_vdc0["omega"][0]
+    #omega_resolution = data_vdc0["omega"].shape[1]
+    vac = data_omega5["vac"][:,0]
+    #assert np.allclose(vac, omega_max*data_vdc0["vac_omega"][:,0])
+
+    data = dm.list(omega=16.5372, good_flags=0x1000, bad_flags=0x2c8c, voltage_branches=4, method="mu", d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, padding=0, xL=0.5)
+    #data = data.loc[(data.vdc < 165.38) & (data.vac < 165.38)]
+    x_arr, y_arr = np.meshgrid(np.linspace(-165.372, 165.372, 601), np.linspace(0, 330.744, 601))
+    gdc_omega5_img_arr = np.pi*griddata(
+            (data.vdc-data.vac, data.vdc+data.vac),
+            data.dc_conductance,
+            (x_arr, y_arr),
+            method="cubic")
+    array = np.array(0xffff*real_cmap(gdc_omega5_img_arr[::-1]), dtype=np.uint16)
+    imwrite(f"figdata/overview_z_rotated.png", array, format="PNG-FI")
 
     fig = plt.figure()
     ax = fig.add_subplot(projection="3d")
@@ -2365,11 +2383,10 @@ def overview_3d_contours(dm, show_contours=False, show_surfaces=True, save_image
     cmap_ax.set_xlim(0, 1)
     cmap_ax.set_ylim(min_m, 1)
 
-    vdc0_selection = data_vdc0["omega"][0] <= omega_max + 2e-3
-    vdc0_y = data_vdc0["vac_omega"][:,vdc0_selection]
+    vdc0_y = data_vdc0["vac_omega"]
     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]
+    vdc0_z = data_vdc0["omega"]
+    vdc0_m = np.pi*data_vdc0["gdc_J_o3a_p"]
 
     if show_surfaces:
         ax.contourf(vdc0_m, vdc0_y, vdc0_z/TK_VOLTAGE, zdir="x", levels=fill_levels, offset=0, zorder=-10, colors=colors)
@@ -2377,9 +2394,19 @@ def overview_3d_contours(dm, show_contours=False, show_surfaces=True, save_image
         contour = ax.contour(vdc0_m, vdc0_y, vdc0_z/TK_VOLTAGE, zdir="x", levels=levels, offset=-offset, zorder=10, colors="black")
         save_contour_coordinates(contour, "figdata/overview_contour_x.dat", "vdc", "omega", "gdc")
     if save_images:
-        array = np.array(255*real_cmap(vdc0_m.T[::-1,::-1]), dtype=np.ubyte)
-        img = Image.fromarray(array)
-        img.save("figdata/overview_x.png")
+        array = np.array(0xffff*real_cmap(vdc0_m.T[::-1,::-1]), dtype=np.uint16)
+        imwrite(f"figdata/overview_x.png", array, format="PNG-FI")
+        array_new = np.empty((array.shape[1], array.shape[0]+array.shape[1], 4), dtype=np.uint16)
+        array_new[:,array.shape[0]:].fill(0xffff)
+        array_new[:,array.shape[0]:,3].fill(0)
+        array_new[:,:array.shape[0]] = array.transpose((1,0,2))
+        sheared_array = array_new.flat[:array_new.size - 4*array_new.shape[0]].reshape(array.shape[1], array.shape[0]+array.shape[1]-1, 4).transpose((1,0,2))
+        imwrite(f"figdata/overview_x_sheared.png", sheared_array, format="PNG-FI")
+        #array_new = np.empty((array.shape[0], array.shape[0]+array.shape[1], 4), dtype=np.uint16)
+        #array_new[:,array.shape[1]:].fill(0xffff)
+        #array_new[:,:array.shape[1]] = array
+        #sheared_array = array_new.flat[:array_new.size - 4*array_new.shape[0]].reshape(array.shape[0], array.shape[0]+array.shape[1]-1, 4)
+        #imwrite(f"figdata/overview_x_sheared.png", sheared_array, format="PNG-FI")
 
     omega5 = data_omega5["omega"]
     omega5_x = data_omega5["vdc"]/omega5
@@ -2391,12 +2418,10 @@ def overview_3d_contours(dm, show_contours=False, show_surfaces=True, save_image
         contour = ax.contour(omega5_x, omega5_y, omega5_m, zdir="z", levels=levels, offset=omega5/TK_VOLTAGE+offset, colors="black")
         save_contour_coordinates(contour, "figdata/overview_contour_z.dat", "vac", "vdc", "gdc")
     if save_images:
-        array = np.array(255*real_cmap(omega5_m[::-1]), dtype=np.ubyte)
-        img = Image.fromarray(array)
-        img.save("figdata/overview_z.png")
+        imwrite(f"figdata/overview_z.png", np.array(0xffff*real_cmap(omega5_m[::-1]), dtype=np.uint16), format="PNG-FI")
 
-    vac0_x = np.linspace(0, 10, vdc_resolution)
-    vac0_z = np.linspace(1e-3, omega_max, omega_resolution)
+    vac0_x = vdc/omega_max
+    vac0_z = omega
     vac0_vdc = vac0_x.reshape((-1,1)) * vac0_z.reshape((1,-1))
     vac0_m = np.pi*interp(vac0_vdc)
     if show_surfaces:
@@ -2405,9 +2430,14 @@ def overview_3d_contours(dm, show_contours=False, show_surfaces=True, save_image
         contour = ax.contour(vac0_x, vac0_m, vac0_z/TK_VOLTAGE, zdir="y", levels=levels, offset=-offset, zorder=10, colors="black")
         save_contour_coordinates(contour, "figdata/overview_contour_y.dat", "omega", "vac", "gdc")
     if save_images:
-        array = np.array(255*real_cmap(vac0_m.T[::-1]), dtype=np.ubyte)
-        img = Image.fromarray(array)
-        img.save("figdata/overview_y.png")
+        array = np.array(0xffff*real_cmap(vac0_m.T[::-1]), dtype=np.uint16)
+        imwrite(f"figdata/overview_y.png", array, format="PNG-FI")
+        array_new = np.empty((array.shape[1], array.shape[0]+array.shape[1], 4), dtype=np.uint16)
+        array_new[:,array.shape[0]:].fill(0xffff)
+        array_new[:,array.shape[0]:,3].fill(0)
+        array_new[:,:array.shape[0]] = array[::-1].transpose((1,0,2))
+        sheared_array = array_new.flat[:array_new.size - 4*array_new.shape[0]].reshape(array.shape[1], array.shape[0]+array.shape[1]-1, 4).transpose((1,0,2))[::-1]
+        imwrite(f"figdata/overview_y_sheared.png", sheared_array, format="PNG-FI")
 
     #ax.plot_surface(vdc0_x, vdc0_y, vdc0_z, rstride=1, cstride=1, facecolors=plt.cm.viridis(vdc0_m), shade=False, zorder=-10)
     #ax.set_axis_off()