diff --git a/final_plots.py b/final_plots.py
index df0c95a5a34593a643ebf4f30be5320a87405219..bf19f33d42b73726c69b49ca471b244200293214 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -18,6 +18,7 @@ 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
@@ -234,21 +235,41 @@ def export_omega5():
np.savetxt(file, np.array([vac[i]/omega, vdc[i]/omega, np.pi*gdc[i], idc[i], iac[i]]).T)
-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,
+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,
+ **parameters
):
"""
- Export interpolated data for given 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}_{order}{padding}_{method}_{suffix}
@@ -281,23 +302,20 @@ def export_omega5_interp(
"_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(
- omega = omega,
- vdc = None,
- vac = None,
- method = None,
- d = 1e9,
- xL = 0.5,
- solver_tol_rel = 1e-8,
- solver_tol_abs = 1e-10,
- voltage_branches = 4,
- )
+ 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(["vac", "vdc"])
+ 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,
@@ -321,72 +339,69 @@ def export_omega5_interp(
selected &= data["padding"] == 0
for key, value in kwargs.items():
selected &= (data[key] == value)
+ if special_selection is not None:
+ selected &= special_selection(data, order, method, padding, **kwargs)
return selected
- vac_arr = np.linspace(0, vac_max, ac_res)
- vdc_arr = np.linspace(0, vdc_max, dc_res)
- vdc_mesh, vac_mesh = np.meshgrid(vdc_arr, vac_arr)
-
def gen_data(order, method, padding=False, s=5e-5, **kwargs):
suffix = "_p" if padding else ""
reduced_data = data.loc[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 \
- or reduced_data.vdc.max() - reduced_data.vdc.min() < 10 \
- or reduced_data.vac.max() - reduced_data.vac.min() < 10:
+ if reduced_data.shape[0] < 100:
return results
+ xy_data = (x_func(reduced_data), y_func(reduced_data))
try:
results[f"gdc_{method}_{order}{suffix}"] = griddata(
- (reduced_data.vdc, reduced_data.vac),
+ xy_data,
reduced_data.dc_conductance,
- (vdc_mesh, vac_mesh),
+ mesh,
method="cubic")
- #gdc_tck = bisplrep(reduced_data.vac, reduced_data.vdc, reduced_data.dc_conductance, s=s, kx=korder, ky=korder)
- #results[f"gdc_{method}_{order}{suffix}"] = bisplev(vac_arr, vdc_arr, gdc_tck)
+ #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(
- (reduced_data.vdc, reduced_data.vac),
+ xy_data,
reduced_data.dc_current,
- (vdc_mesh, vac_mesh),
+ mesh,
method="cubic")
except:
settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method}")
try:
- idc_tck = bisplrep(reduced_data.vac, reduced_data.vdc, reduced_data.dc_current, s=s, kx=korder, ky=korder)
- results[f"idc_{method}_{order}{suffix}_spline"] = bisplev(vac_arr, vdc_arr, idc_tck)
- results[f"gdc_{method}_{order}{suffix}_ispline"] = bisplev(vac_arr, vdc_arr, idc_tck, dy=1)
+ 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(
- (reduced_data.vdc, reduced_data.vac),
+ xy_data,
2*reduced_data.ac_current_abs,
- (vdc_mesh, vac_mesh),
+ mesh,
method="cubic")
except:
settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
try:
- iac_tck = bisplrep(reduced_data.vac, reduced_data.vdc, 2*reduced_data.ac_current_abs, s=s, kx=korder, ky=korder)
- results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(vac_arr, vdc_arr, iac_tck)
- results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(vac_arr, vdc_arr, iac_tck, dx=1)
+ 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(
- (reduced_data.vdc, reduced_data.vac),
+ xy_data,
reduced_data.ac_current_phase,
- (vdc_mesh, vac_mesh),
+ 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(
- (reduced_data.vdc, reduced_data.vac),
+ xy_data,
reduced_data.gamma,
- (vdc_mesh, vac_mesh),
+ mesh,
method="cubic")
except:
settings.logger.exception(f"in griddata interpolation of gamma for {order}{suffix} {method}")
@@ -394,8 +409,7 @@ def export_omega5_interp(
np.savez(
filename,
- vdc = vdc_arr,
- vac = vac_arr,
+ **{fixed_parameter:fixed_value, x_parameter:x_mesh, y_parameter:y_mesh},
**gen_data("o2", "mu"),
**gen_data("o3", "mu"),
**gen_data("o3a", "mu"),
@@ -411,6 +425,67 @@ def export_omega5_interp(
)
+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,
+ ):
+ """
+ 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)
+ 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.05,
+ omega_max = 20,
+ vac_omega_min = 0,
+ vac_omega_max = 10,
+ omega_res = 101,
+ vac_res = 101,
+ 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
diff --git a/plot.py b/plot.py
index 8611536da8675228cb338518bbd112f202295add..9c90fd1bf934ec2fd86cd904429f3dfdf6f5ae95 100644
--- a/plot.py
+++ b/plot.py
@@ -1721,6 +1721,27 @@ def RGeq_comparison_contours(
plt.show()
+def contours_to_coordinates(contour):
+ for level, collection in zip(contour.levels, contour.collections):
+ for path in collection.get_paths():
+ yield level, path.vertices
+
+def contours_to_coordinates_unified(contour):
+ for level, vertices in contours_to_coordinates(contour):
+ array = np.empty((vertices.shape[0], 3), dtype=np.float64)
+ array[:,:2] = vertices
+ array[:,2] = level
+ yield array
+
+def save_contour_coordinates(contour, filename, x, y, z):
+ """For usage in pgf plots. Overwrites filename."""
+ with open(filename, "w") as file:
+ print(f"{x} {y} {z}", file=file)
+ for array in contours_to_coordinates_unified(contour):
+ np.savetxt(file, array, "%.6f")
+ print(file=file)
+
+
def RGeq_comparison_contours2(
dm = None,
filename = "figdata/omega5_interp.npz",
@@ -1738,13 +1759,16 @@ def RGeq_comparison_contours2(
levels = (1/np.linspace(1, 12, 12)[::-1]) / prefactor
fig, ax = plt.subplots(nrows=1, ncols=1, sharex=True, sharey=True)
for approx in approximations:
- ax.contour(vdc, vac, data[f"{observable}_mu_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+ contour = ax.contour(vdc, vac, data[f"{observable}_mu_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+ save_contour_coordinates(contour, f"figdata/contour_{observable}_mu_{approx}.dat", "vdc", "vac", "gdc")
try:
- ax.contour(vdc, vac, data[f"{observable}_J_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+ contour = ax.contour(vdc, vac, data[f"{observable}_J_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+ save_contour_coordinates(contour, f"figdata/contour_{observable}_J_{approx}.dat", "vdc", "vac", "gdc")
except KeyError:
settings.logger.warning(f"No data available for {observable}_J_{approx}")
try:
- ax.contour(vdc, vac, data[f"{observable}_J_{approx}_p"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+ contour = ax.contour(vdc, vac, data[f"{observable}_J_{approx}_p"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+ save_contour_coordinates(contour, f"figdata/contour_{observable}_J_{approx}_p.dat", "vdc", "vac", "gdc")
except KeyError:
settings.logger.warning(f"No data available for {observable}_J_{approx}_p")