From d8ab37ef50ac6d875a1367ff0af90342a4d3ac10 Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Thu, 29 Dec 2022 00:58:55 +0100
Subject: [PATCH] more plots, small changes in some plots
---
final_plots.py | 3 +-
gen_data.py | 2 +-
plot.py | 216 +++++++++++++++++++++++++++++++++++++++++--------
plot_pulses.py | 40 +++++++--
4 files changed, 216 insertions(+), 45 deletions(-)
diff --git a/final_plots.py b/final_plots.py
index b48fab8..0150b1d 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -455,11 +455,12 @@ def export_omega5_interp(
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=4, xL=0.5)
+ 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)
diff --git a/gen_data.py b/gen_data.py
index 09478ed..1ef10e0 100644
--- a/gen_data.py
+++ b/gen_data.py
@@ -196,7 +196,7 @@ python -m frtrg_kondo.gen_data \\
parallel_group.add_argument("--steps", metavar="int", type=int, nargs='+',
help = "Number of steps, provided for each independent parameter swipe dimension")
parallel_group.add_argument("--scale", type=str, nargs='+', default="linear",
- choices = ("linear", "log"),
+ choices = ("linear", "lin", "log"),
help = "Scale used for swipes (must get same number of options as --steps)")
parallel_group.add_argument("--threads", type=int, metavar="int", default=4,
help = "Number parallel processes (set to 0 to use all CPUs)")
diff --git a/plot.py b/plot.py
index f6e3ffd..614054d 100644
--- a/plot.py
+++ b/plot.py
@@ -7,7 +7,7 @@ Kondo FRTRG, generate plots based on save data
"""
import matplotlib.pyplot as plt
-import matplotlib.colors as mplcolors
+from matplotlib.colors import LogNorm, FuncNorm
import argparse
import pandas as pd
import numpy as np
@@ -529,10 +529,30 @@ def plot_comparison(dm, omega, d=1e9, **parameters):
Only data points which exist for both methods with equal Vdc, Vac, D and
frequency are considered.
"""
- for name in ("omega", "method"):
+ #for name in ("omega", "method"):
+ # parameters.pop(name, None)
+ #results_J = dm.list(omega=omega, method='J', d=d, **parameters).sort_values(["vdc","vac"])
+ #results_mu = dm.list(omega=omega, method='mu', d=d, **parameters).sort_values(["vdc","vac"])
+ ## For comparison with resonant dc voltage:
+ for name in ("omega", "method", "voltage_branches"):
parameters.pop(name, None)
- results_J = dm.list(omega=omega, method='J', d=d, **parameters).sort_values(["vdc","vac"])
- results_mu = dm.list(omega=omega, method='mu', d=d, **parameters).sort_values(["vdc","vac"])
+ parameters["include_Ga"] = False
+ parameters["integral_method"] = -15
+ parameters["truncation_order"] = 3
+ parameters["solver_tol_rel"] = 1e-8
+ parameters["solver_tol_abs"] = 1e-10
+ results_J = dm.list(omega=omega, method='mu', d=d, voltage_branches=4, **parameters).sort_values(["vdc","vac"])
+ results_mu = dm.list(omega=omega, method='J', d=d, voltage_branches=0, **parameters).sort_values(["vdc","vac"])
+ ## For comparison of different numbers of voltage branches:
+ #for name in ("omega", "method", "voltage_branches"):
+ # parameters.pop(name, None)
+ #parameters["include_Ga"] = True
+ #parameters["integral_method"] = -15
+ #parameters["truncation_order"] = 3
+ #parameters["solver_tol_rel"] = 1e-8
+ #parameters["solver_tol_abs"] = 1e-10
+ #results_J = dm.list(omega=omega, method='mu', d=d, voltage_branches=4, **parameters).sort_values(["vdc","vac"])
+ #results_mu = dm.list(omega=omega, method='mu', d=d, voltage_branches=7, **parameters).sort_values(["vdc","vac"])
results_J.vac = results_J.vac.round(6)
results_J.vdc = results_J.vdc.round(6)
results_mu.vac = results_mu.vac.round(6)
@@ -586,10 +606,21 @@ def plot_comparison_relative(dm, omega, d=1e9, **parameters):
Only data points which exist for both methods with equal Vdc, Vac, D and
frequency are considered.
"""
- for name in ("omega", "method"):
+ #for name in ("omega", "method"):
+ # parameters.pop(name, None)
+ #results_J = dm.list(omega=omega, method='J', d=d, **parameters).sort_values(["vdc","vac"])
+ #results_mu = dm.list(omega=omega, method='mu', d=d, **parameters).sort_values(["vdc","vac"])
+ ## For comparison with resonant dc voltage:
+ for name in ("omega", "method", "voltage_branches"):
parameters.pop(name, None)
- results_J = dm.list(omega=omega, method='J', d=d, **parameters).sort_values(["vdc","vac"])
- results_mu = dm.list(omega=omega, method='mu', d=d, **parameters).sort_values(["vdc","vac"])
+ parameters["include_Ga"] = False
+ parameters["integral_method"] = -15
+ parameters["truncation_order"] = 3
+ parameters["solver_tol_rel"] = 1e-8
+ parameters["solver_tol_abs"] = 1e-10
+ results_J = dm.list(omega=omega, method='mu', d=d, voltage_branches=4, **parameters).sort_values(["vdc","vac"])
+ results_mu = dm.list(omega=omega, method='J', d=d, voltage_branches=0, **parameters).sort_values(["vdc","vac"])
+
results_J.vac = results_J.vac.round(6)
results_J.vdc = results_J.vdc.round(6)
results_mu.vac = results_mu.vac.round(6)
@@ -644,55 +675,55 @@ def plot_floquet_matrices(kondo : KondoImport, norm_min=1e-6):
idx = kondo.voltage_branches if gamma.ndim == 3 else ...
try:
axes[0].set_title("Γ")
- img = axes[0].imshow(np.abs(gamma[idx]), norm=mplcolors.LogNorm(norm_min))
+ img = axes[0].imshow(np.abs(gamma[idx]), norm=LogNorm(norm_min))
fig.colorbar(img, ax=axes[0])
except AttributeError:
pass
try:
axes[1].set_title("Z")
- img = axes[1].imshow(np.abs(kondo.z[idx]), norm=mplcolors.LogNorm(norm_min))
+ img = axes[1].imshow(np.abs(kondo.z[idx]), norm=LogNorm(norm_min))
fig.colorbar(img, ax=axes[1])
except AttributeError:
pass
try:
axes[2].set_title("ΓL")
- img = axes[2].imshow(np.abs(kondo.gammaL), norm=mplcolors.LogNorm(norm_min))
+ img = axes[2].imshow(np.abs(kondo.gammaL), norm=LogNorm(norm_min))
fig.colorbar(img, ax=axes[2])
except AttributeError:
pass
try:
axes[3].set_title("δΓL")
- img = axes[3].imshow(np.abs(kondo.deltaGammaL), norm=mplcolors.LogNorm(norm_min))
+ img = axes[3].imshow(np.abs(kondo.deltaGammaL), norm=LogNorm(norm_min))
fig.colorbar(img, ax=axes[3])
except AttributeError:
pass
try:
axes[4].set_title("δΓ")
- img = axes[4].imshow(np.abs(kondo.deltaGamma[1 if gamma.ndim == 3 else ...]), norm=mplcolors.LogNorm(norm_min))
+ img = axes[4].imshow(np.abs(kondo.deltaGamma[1 if gamma.ndim == 3 else ...]), norm=LogNorm(norm_min))
fig.colorbar(img, ax=axes[4])
except AttributeError:
pass
try:
axes[5].set_title("yL")
- img = axes[5].imshow(np.abs(kondo.yL), norm=mplcolors.LogNorm(norm_min))
+ img = axes[5].imshow(np.abs(kondo.yL), norm=LogNorm(norm_min))
fig.colorbar(img, ax=axes[5])
except AttributeError:
pass
try:
axes[6].set_title("G2")
- img = axes[6].imshow(np.abs(kondo.g2[:,:,idx].transpose(0,2,1,3).reshape((4*kondo.nmax+2, 4*kondo.nmax+2))), norm=mplcolors.LogNorm(norm_min))
+ img = axes[6].imshow(np.abs(kondo.g2[:,:,idx].transpose(0,2,1,3).reshape((4*kondo.nmax+2, 4*kondo.nmax+2))), norm=LogNorm(norm_min))
fig.colorbar(img, ax=axes[6])
except AttributeError:
pass
try:
axes[7].set_title("G3")
- img = axes[7].imshow(np.abs(kondo.g3[:,:,idx].transpose(0,2,1,3).reshape((4*kondo.nmax+2, 4*kondo.nmax+2))), norm=mplcolors.LogNorm(norm_min))
+ img = axes[7].imshow(np.abs(kondo.g3[:,:,idx].transpose(0,2,1,3).reshape((4*kondo.nmax+2, 4*kondo.nmax+2))), norm=LogNorm(norm_min))
fig.colorbar(img, ax=axes[7])
except AttributeError:
pass
try:
axes[8].set_title("I")
- img = axes[8].imshow(np.abs(kondo.current.transpose(0,2,1,3).reshape((4*kondo.nmax+2, 4*kondo.nmax+2))), norm=mplcolors.LogNorm(norm_min))
+ img = axes[8].imshow(np.abs(kondo.current.transpose(0,2,1,3).reshape((4*kondo.nmax+2, 4*kondo.nmax+2))), norm=LogNorm(norm_min))
fig.colorbar(img, ax=axes[8])
except AttributeError:
pass
@@ -809,14 +840,14 @@ def fixed_omega_dconvergence(dm, omega=16.5372, solver_tol_rel=1e-10, solver_tol
vac_grid,
idc_enhanced_grid,
c=d_grid,
- norm=mplcolors.LogNorm(),
+ norm=LogNorm(),
cmap=plt.cm.viridis_r,
marker="x")
plot = ax_iac.scatter(
vdc_grid,
2*iac_enhanced_grid,
c=d_grid,
- norm=mplcolors.LogNorm(),
+ norm=LogNorm(),
cmap=plt.cm.viridis_r,
marker="x")
cb = fig.colorbar(plot, ax=[ax_idc, ax_iac])
@@ -872,7 +903,7 @@ def check_convergence(dm, **parameters):
x_mu = x[mu]
x_mod = x[mod]
drtol = table.d*table.solver_tol_rel
- norm = mplcolors.LogNorm(max(drtol.min(), 0.05), min(drtol.max(), 500))
+ norm = LogNorm(max(drtol.min(), 0.05), min(drtol.max(), 500))
c_j = drtol[j]
c_mu = drtol[mu]
c_mod = drtol[mod]
@@ -948,7 +979,7 @@ def check_convergence_nmax(dm, **parameters):
mod = (table.solver_flags & DataManager.SOLVER_FLAGS["simplified_initial_conditions"]) != 0
j = (~mod) & (table.method == "J")
mu = (~mod) & (table.method == "mu")
- norm = mplcolors.LogNorm(table.voltage_branches.min(), table.voltage_branches.max())
+ norm = LogNorm(table.voltage_branches.min(), table.voltage_branches.max())
nmax_j = table.nmax[j]
nmax_mu = table.nmax[mu]
nmax_mod = table.nmax[mod]
@@ -1039,7 +1070,7 @@ def check_convergence_fit(dm, **parameters):
x_mu = x[mu]
x_mod = x[mod]
drtol = table.d*table.solver_tol_rel
- norm = mplcolors.LogNorm(max(drtol.min(), 0.05), min(drtol.max(), 500))
+ norm = LogNorm(max(drtol.min(), 0.05), min(drtol.max(), 500))
c_j = drtol[j]
c_mu = drtol[mu]
c_mod = drtol[mod]
@@ -1559,6 +1590,7 @@ def compare_RGeq_interp(
ref_method = "mu",
observable = "gdc",
suffix = "",
+ ref_filename = None,
**trashparams):
data = np.load(filename)
prefactor = 1
@@ -1590,7 +1622,11 @@ def compare_RGeq_interp(
ax1.set_xlim(vdc[0], vdc[-1])
ax1.set_ylim(vac[0], vac[-1])
- ref_data = data[f"{observable}_{ref_method}_o3p{suffix}"]
+ if ref_filename is not None:
+ refdata = np.load(ref_filename)
+ ref_data = refdata[f"{observable}_{ref_method}_o3p{suffix}"]
+ else:
+ ref_data = data[f"{observable}_{ref_method}_o3p{suffix}"]
for ax in axes.flat:
ax.contour(vdc, vac, ref_data, levels, **ref_contour_kwargs)
@@ -1646,13 +1682,33 @@ def compare_RGeq_interp(
# ax.contour(vdc, vac, ref_data, missing_levels, **contour_kwargs)
+class GNorm(plt.Normalize):
+ def __init__(self, vmin=0.1, vmax=1, clip=False):
+ assert vmin > 0
+ assert vmax > vmin
+ self._vmin = vmin
+ self._vmax = vmax
+ self._clip = clip
+
+ def __call__(self, value, clip=None):
+ return (1 - self._vmin/value) / (1 - self._vmin/self._vmax)
+
+def GNorm(vmin, vmax):
+ assert vmax > vmin > 0
+ scale = 1 - vmin/vmax
+ return FuncNorm(((lambda x: (1 - vmin/x)/scale), (lambda x: vmin/(1 - scale*x))), vmin, vmax)
+
+
def RGeq_comparison_order2(
dm = None,
filename = "figdata/omega5_interp.npz",
- method = "J",
- compare = "o3a_p",
- reference = "o3a",
+ method = "mu",
+ compare = "o3a",
+ #compare = "o3a",
+ reference = "o3p",
reference_method = "mu",
+ #ref_filename = "figdata/omega5_interp_vb7.npz",
+ ref_filename = "figdata/omega5_interp.npz",
observable = "gdc",
prefactor = np.pi,
**trashparams):
@@ -1665,14 +1721,18 @@ def RGeq_comparison_order2(
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, sharex=True, sharey=True)
data = np.load(filename)
cmp = prefactor*data[f"{observable}_{method}_{compare}"]
- ref = prefactor*data[f"{observable}_{reference_method}_{reference}"]
+ if ref_filename is None:
+ ref = prefactor*data[f"{observable}_{reference_method}_{reference}"]
+ else:
+ ref_data = np.load(ref_filename)
+ ref = prefactor*ref_data[f"{observable}_{reference_method}_{reference}"]
#gdc_o3 = np.pi*data["gdc_mu_o3"]
- #norm = mplcolors.LogNorm(0.08, 0.4)
- norm = plt.Normalize(0.08, 0.4)
+ #norm = LogNorm(0.08, 0.4)
+ norm = GNorm(0.08, 1)
cmap = plt.cm.viridis
omega = 16.5372
- vdc = data["vdc"] / omega
- vac = data["vac"] / omega
+ vdc = data["vdc"][0] / omega
+ vac = data["vac"][:,0] / omega
extent1 = (1.5*vdc[0] - 0.5*vdc[1], 1.5*vdc[-1] - 0.5*vdc[-2], 1.5*vac[0] - 0.5*vac[1], 1.5*vac[-1] - 0.5*vac[-2])
extent2 = (1.5*vdc[0] - 0.5*vdc[1], -1.5*vdc[-1] + 0.5*vdc[-2], 1.5*vac[0] - 0.5*vac[1], 1.5*vac[-1] - 0.5*vac[-2])
img1 = ax1.imshow(ref, cmap=cmap, norm=norm, extent=extent1, origin="lower")
@@ -1703,16 +1763,34 @@ def RGeq_comparison_contours(
reference = "o3p",
observable = "gdc",
prefactor = np.pi,
+ comp1_filename = None,
+ comp2_filename = None,
**trashparams):
+ """
+ Example usage:
+ >>> RGeq_comparison_contours(comp1_filename="figdata/omega5_interp_vb7.npz", compare1="o3a", compare2="o3a", reference="o3a")
+ """
contour_kwargs = dict(colors="#000000", linewidths=0.8, zorder=2)
ref_contour_kwargs = dict(colors="#ff0000", linewidths=0.8, zorder=1)
data = np.load(filename)
- cmp1 = prefactor*data[f"{observable}_{method}_{compare1}"]
- cmp2 = prefactor*data[f"{observable}_{method}_{compare2}"]
+ if comp1_filename is not None:
+ comp1_data = np.load(comp1_filename)
+ assert np.allclose(data["vdc"], comp1_data["vdc"])
+ assert np.allclose(data["vac"], comp1_data["vac"])
+ cmp1 = prefactor*comp1_data[f"{observable}_{method}_{compare1}"]
+ else:
+ cmp1 = prefactor*data[f"{observable}_{method}_{compare1}"]
+ if comp2_filename is not None:
+ comp2_data = np.load(comp2_filename)
+ assert np.allclose(data["vdc"], comp2_data["vdc"])
+ assert np.allclose(data["vac"], comp2_data["vac"])
+ cmp2 = prefactor*comp2_data[f"{observable}_{method}_{compare2}"]
+ else:
+ cmp2 = prefactor*data[f"{observable}_{method}_{compare2}"]
ref = prefactor*data[f"{observable}_{method}_{reference}"]
omega = 16.5372
- vdc = data["vdc"] / omega
- vac = data["vac"] / omega
+ vdc = data["vdc"][0] / omega
+ vac = data["vac"][:,0] / omega
levels = 1/np.linspace(1, 12, 12)[::-1]
#levels = 0.05 + np.linspace(0.03**0.1, 0.9**.1, 20)**10
#levels = np.logspace(-2,0,50)
@@ -1994,7 +2072,7 @@ def matrix_truncation(
kondos[idx] = kondo
fig, axes = plt.subplots(nrows=2, ncols=4, sharex=True, sharey=True)
- norm = mplcolors.LogNorm(1e-9, 1)
+ norm = LogNorm(1e-9, 1)
cmap = plt.cm.viridis
#z_ref = (z_matrices[("mu", True, False)] + z_matrices[("J", True, True)])/2
z_ref_j = z_matrices[("J", True, True)]
@@ -2263,5 +2341,73 @@ def overview_3d_contours(dm, show_contours=False, show_surfaces=True, save_image
plt.show()
+def RGeq_comparison_new(
+ dm = None,
+ filename = "figdata/omega5_interp.npz",
+ method = "mu",
+ compare = "o2",
+ reference = "o3a",
+ reference_method = "mu",
+ #ref_filename = "figdata/omega5_interp_vb7.npz",
+ ref_filename = "figdata/omega5_interp.npz",
+ observable = "gdc",
+ prefactor = np.pi,
+ **trashparams):
+ """
+ Compare two different results for the differential conductance
+ (or another observable) side-by-side and by plotting their absolute
+ and relative difference.
+ Used to generate some draft figures for the thesis.
+ """
+ fig, ((ax_dc, ax_trash), (ax_img, ax_ac)) = plt.subplots(nrows=2, ncols=2, sharex="col", sharey="row")
+ data = np.load(filename)
+ cmp = prefactor*data[f"{observable}_{method}_{compare}"]
+ if ref_filename is None:
+ ref = prefactor*data[f"{observable}_{reference_method}_{reference}"]
+ else:
+ ref_data = np.load(ref_filename)
+ ref = prefactor*ref_data[f"{observable}_{reference_method}_{reference}"]
+ #gdc_o3 = np.pi*data["gdc_mu_o3"]
+ #norm = LogNorm(0.08, 0.4)
+
+ omega = 16.5372
+ vdc = data["vdc"][0] / omega
+ vac = data["vac"][:,0] / omega
+
+ skip = 40
+ plot_number_dc = data["vdc"].shape[0] // skip + 1
+ plot_number_ac = data["vac"].shape[1] // skip + 1
+ vac_max = vac[-1]
+ vdc_max = vdc[-1]
+ for i in range(plot_number_dc):
+ color = plt.cm.viridis(data["vac"][i*skip,0]/(omega*vac_max))
+ ax_dc.plot(data["vdc"][i*skip]/omega, cmp[i*skip], color=color)
+ ax_dc.plot(data["vdc"][i*skip]/omega, ref[i*skip], color=color, linestyle=":")
+ for i in range(plot_number_ac):
+ color = plt.cm.viridis(data["vdc"][0,i*skip]/(omega*vdc_max))
+ ax_ac.plot(cmp[:,i*skip], data["vac"][:,i*skip]/omega, color=color)
+ ax_ac.plot(ref[:,i*skip], data["vac"][:,i*skip]/omega, color=color, linestyle=":")
+
+ #ax_dc.plot(data["vdc"][::skip].T/omega, cmp[::skip].T, color="blue")
+ #ax_dc.plot(data["vdc"][::skip].T/omega, ref[::skip].T, ":", color="red")
+ #ax_ac.plot(cmp[:,::skip], data["vac"][:,::skip]/omega, color="blue")
+ #ax_ac.plot(ref[:,::skip], data["vac"][:,::skip]/omega, color="red")
+
+ extent = (1.5*vdc[0] - 0.5*vdc[1], 1.5*vdc[-1] - 0.5*vdc[-2], 1.5*vac[0] - 0.5*vac[1], 1.5*vac[-1] - 0.5*vac[-2])
+ ccmap = plt.cm.seismic
+ absdiff = cmp - ref
+ absdiff_max = np.abs(absdiff).max()
+ reldiff = absdiff/ref
+ reldiff_max = np.abs(reldiff).max()
+ img = ax_img.imshow(reldiff, cmap=ccmap, norm=plt.Normalize(-reldiff_max, reldiff_max), extent=extent, origin="lower")
+ #img4 = ax2.imshow(absdiff, cmap=ccmap, norm=plt.Normalize(-absdiff_max, absdiff_max), extent=extent2, origin="lower")
+ fig.colorbar(img, ax=ax_trash, location="right")
+ ax_img.set_xlabel(r"$V_\mathrm{avg}~(\Omega)$")
+ ax_img.set_ylabel(r"$V_\mathrm{osc}~(\Omega)$")
+ ax_ac.set_xlabel(r"$G~(2e^2/h)$")
+ ax_dc.set_ylabel(r"$G~(2e^2/h)$")
+ plt.show()
+
+
if __name__ == '__main__':
main()
diff --git a/plot_pulses.py b/plot_pulses.py
index f0869c3..43f3c4c 100644
--- a/plot_pulses.py
+++ b/plot_pulses.py
@@ -108,18 +108,28 @@ def integrate_ft(t0, t1, fourier_coef):
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)
+def plot_transported_charge(dm, omega=None, **parameters):
+ if omega<=0:
+ omega = None
+ table = load_all_pulses(dm, omega=omega)
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")
+ fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1)
+ ax1.set_xlabel("pulse duration t (1/Tk)")
+ ax2.set_xlabel("pulse phase φ/2π")
+ ax1.set_ylabel("charge per pulse Q (e)")
+ ax2.set_ylabel("charge per pulse Q (e)")
+
+ phase_norm = plt.Normalize(0, 1.25)
+ duration_norm = LogNorm(0.01, 0.2)
+ ax1.scatter(table.pulse_duration, table.dc_current/table.omega*2*np.pi, c=table.pulse_phase, s=100/table.omega**2, marker="o", norm=phase_norm)
+ ax2.scatter(table.pulse_phase, table.dc_current/table.omega*2*np.pi, c=table.pulse_duration, s=100/table.omega**2, marker="o", norm=duration_norm)
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"):
+ for (row, kondo) in load_all_pulses_full(dm, pulse_type="gauss_symmetric", omega=omega):
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)
@@ -130,7 +140,21 @@ def plot_transported_charge(dm, **parameters):
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)
+ sym_pulse_duration = np.array(sym_pulse_duration)
+ sorting = np.argsort(sym_pulse_duration)
+ sym_pulse_duration = sym_pulse_duration[sorting]
+ sym_pulse_omega = np.array(sym_pulse_omega)[sorting]
+ sym_pulse_phase = np.array(sym_pulse_phase)[sorting]
+ sym_pulse_charge = np.array(sym_pulse_charge)[sorting]
+
+ ax1.scatter(sym_pulse_duration, sym_pulse_charge, c=sym_pulse_phase, marker="^", s=80/np.array(sym_pulse_omega)**2, norm=phase_norm)
+ for phase in (0.25, 0.5, 0.75, 1):
+ selection = np.abs(sym_pulse_phase - phase) < 1e-3
+ ax1.plot(sym_pulse_duration[selection], sym_pulse_charge[selection], color="black")
+ ax2.scatter(sym_pulse_phase, sym_pulse_charge, c=sym_pulse_duration, marker="^", s=80/np.array(sym_pulse_omega)**2, norm=duration_norm)
+ selection = np.abs(sym_pulse_duration - 0.04) < 1e-4
+ sorting = np.argsort(sym_pulse_phase[selection])
+ ax2.plot(sym_pulse_phase[selection][sorting], sym_pulse_charge[selection][sorting], color="black")
plt.show()
@@ -213,8 +237,8 @@ def parameter_overview(dm, **trash):
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")
+ ax_gauss.plot(gauss_duration[gauss_good], gauss_phase[gauss_good], "x")
+ ax_gauss.plot(gauss_duration[~gauss_good], gauss_phase[~gauss_good], "+")
plt.show()
--
GitLab