|
@@ -1,10 +1,11 @@
|
|
|
import numpy as np
|
|
|
from charged_shells.rotational_path import PairRotationalPath, PathEnergyPlot
|
|
|
-from charged_shells import expansion, patch_size, interactions
|
|
|
+from charged_shells import expansion
|
|
|
from charged_shells.parameters import ModelParams
|
|
|
import matplotlib.pyplot as plt
|
|
|
from pathlib import Path
|
|
|
import json
|
|
|
+import quadrupole_model_mappings
|
|
|
|
|
|
Array = np.ndarray
|
|
|
|
|
@@ -14,12 +15,12 @@ pi_half_to_pi = np.linspace(np.pi/2, np.pi, 100, endpoint=True)
|
|
|
|
|
|
QuadPath = PairRotationalPath()
|
|
|
QuadPath.set_default_x_axis(zero_to_pi_half)
|
|
|
-QuadPath.add_euler(beta1=np.pi/2, beta2=zero_to_pi_half)
|
|
|
-QuadPath.add_euler(beta1=zero_to_pi_half[::-1], beta2=zero_to_pi_half[::-1])
|
|
|
-QuadPath.add_euler(beta1=zero_to_pi_half)
|
|
|
-QuadPath.add_euler(beta1=zero_to_pi_half[::-1], beta2=zero_to_pi_half)
|
|
|
-QuadPath.add_euler(beta1=np.pi/2, beta2=zero_to_pi_half, alpha2=np.pi/2)
|
|
|
-QuadPath.add_euler(beta1=np.pi/2, beta2=np.pi/2, alpha1=zero_to_pi_half[::-1])
|
|
|
+QuadPath.add_euler(beta1=np.pi/2, beta2=zero_to_pi_half, start_name="EP", end_name="EE")
|
|
|
+QuadPath.add_euler(beta1=zero_to_pi_half[::-1], beta2=zero_to_pi_half[::-1], end_name="PP")
|
|
|
+QuadPath.add_euler(beta1=zero_to_pi_half, end_name="EP")
|
|
|
+QuadPath.add_euler(beta1=zero_to_pi_half[::-1], beta2=zero_to_pi_half, end_name="EP")
|
|
|
+QuadPath.add_euler(beta1=np.pi/2, beta2=zero_to_pi_half, alpha2=np.pi/2, end_name="tEE")
|
|
|
+QuadPath.add_euler(beta1=np.pi/2, beta2=np.pi/2, alpha1=zero_to_pi_half[::-1], end_name="EE")
|
|
|
|
|
|
DipolePath = PairRotationalPath()
|
|
|
DipolePath.set_default_x_axis(zero_to_pi_half)
|
|
@@ -46,34 +47,58 @@ DipolePath2.add_euler(beta2=zero_to_pi_half, beta1=pi_half_to_pi[::-1], alpha1=n
|
|
|
DipolePath2.add_euler(beta2=pi_half_to_pi, beta1=zero_to_pi_half[::-1], alpha1=np.pi)
|
|
|
|
|
|
|
|
|
-def model_comparison(save_as=None):
|
|
|
- params = ModelParams(R=150, kappaR=3)
|
|
|
- abar = 0.5
|
|
|
+def model_comparison(config_data: dict, save_as=None, save_data=False):
|
|
|
+ kappaR = 3
|
|
|
+ params = ModelParams(R=150, kappaR=kappaR)
|
|
|
+ a_bar = 0.5
|
|
|
sigma_tilde = 0.001
|
|
|
|
|
|
- ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
|
|
|
+ ex1 = expansion.MappedExpansionQuad(a_bar, params.kappaR, sigma_tilde, l_max=30)
|
|
|
ex2 = ex1.clone()
|
|
|
|
|
|
+ # matching other models to the mapped CSp model based on equal patch size in potential
|
|
|
+ ex_gauss = quadrupole_model_mappings.ic_to_gauss(sigma_tilde, a_bar, params, l_max=30, sigma0=0)
|
|
|
+ ex_gauss2 = ex_gauss.clone()
|
|
|
+ ex_cap = quadrupole_model_mappings.ic_to_cap(sigma_tilde, a_bar, params, l_max=30, sigma0=0)
|
|
|
+ ex_cap2 = ex_cap.clone()
|
|
|
+
|
|
|
+ # path plots for all models
|
|
|
path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params)
|
|
|
energy = path_plot.evaluate_path()
|
|
|
x_axis = path_plot.rot_path.stack_x_axes()
|
|
|
|
|
|
+ path_plot_gauss = PathEnergyPlot(ex_gauss, ex_gauss2, QuadPath, dist=2., params=params)
|
|
|
+ energy_gauss = path_plot_gauss.evaluate_path()
|
|
|
+
|
|
|
+ path_plot_cap = PathEnergyPlot(ex_cap, ex_cap2, QuadPath, dist=2., params=params)
|
|
|
+ energy_cap = path_plot_cap.evaluate_path()
|
|
|
+
|
|
|
# peak_energy_sanity_check
|
|
|
- ex1new = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
|
|
|
- ex2new = ex1new.clone()
|
|
|
- pp_energy = interactions.charged_shell_energy(ex1new, ex2new, params)
|
|
|
- print(f'PP energy: {pp_energy}')
|
|
|
+ # ex1new = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
|
|
|
+ # ex2new = ex1new.clone()
|
|
|
+ # pp_energy = interactions.charged_shell_energy(ex1new, ex2new, params)
|
|
|
+ # print(f'PP energy: {pp_energy}')
|
|
|
|
|
|
# Emanuele data
|
|
|
- with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
|
|
|
- config_data = json.load(config_file)
|
|
|
em_data = np.load(Path(config_data["emanuele_data"]).joinpath("FIG_7").joinpath("pathway.npz"))['arr_0']
|
|
|
+ # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_8").joinpath("FIXEDCHARGE")
|
|
|
+ # .joinpath("FIX_A").joinpath("ECC_0.25"))
|
|
|
+ # em_data = np.load(em_data_path.joinpath(f"EMME_{kappaR}.").joinpath("pathway.npz"))['arr_0']
|
|
|
+
|
|
|
+ if save_data:
|
|
|
+ np.savez(Path(config_data["figure_data"]).joinpath(f"fig_7_kR{kappaR}.npz"),
|
|
|
+ ICi=em_data,
|
|
|
+ CSp=np.stack((x_axis, np.squeeze(energy))).T,
|
|
|
+ CSp_gauss=np.stack((x_axis, np.squeeze(energy_gauss))).T,
|
|
|
+ CSp_cap=np.stack((x_axis, np.squeeze(energy_cap))).T)
|
|
|
|
|
|
fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
|
|
|
ax.plot(em_data[:, 0], em_data[:, 1], label='ICi')
|
|
|
ax.plot(x_axis, np.squeeze(energy), label='CSp')
|
|
|
+ # ax.plot(x_axis, np.squeeze(energy_gauss), label='CSp - Gauss')
|
|
|
+ # ax.plot(x_axis, np.squeeze(energy_cap), label='CSp - cap')
|
|
|
# ax.plot(x_axis, em_data[:, 1] / np.squeeze(energy), label='CSp')
|
|
|
- PathEnergyPlot.plot_style(fig, ax)
|
|
|
+ path_plot.plot_style(fig, ax)
|
|
|
if save_as is not None:
|
|
|
plt.savefig(save_as, dpi=600)
|
|
|
plt.show()
|
|
@@ -113,7 +138,7 @@ def sigma0_dependence(sigma0: Array, kappaR: float, abar: float, sigma_tilde=0.0
|
|
|
|
|
|
path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=None)
|
|
|
|
|
|
- path_plot.plot(labels=[rf'$\tilde \sigma_0={s0}$' for s0 in sigma0],
|
|
|
+ path_plot.plot(labels=[rf'$\eta={s0 / sigma_tilde}$' for s0 in sigma0],
|
|
|
# norm_euler_angles={'beta2': np.pi},
|
|
|
save_as=save_as)
|
|
|
|
|
@@ -136,7 +161,7 @@ def distance_dependence(dist: Array, kappaR: float, abar: float, sigma_tilde=0.0
|
|
|
fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
|
|
|
for pl, lbl in zip(plots, labels):
|
|
|
ax.plot(x_axis, pl, label=lbl)
|
|
|
- PathEnergyPlot.plot_style(fig, ax)
|
|
|
+ QuadPath.plot_style(fig, ax)
|
|
|
ax.set_ylabel(r'$U \kappa\rho e^{\kappa\rho}$', fontsize=15)
|
|
|
if save_as is not None:
|
|
|
plt.savefig(save_as, dpi=600)
|
|
@@ -157,7 +182,7 @@ def IC_kappaR_dependence(config_data: dict, save_as=None):
|
|
|
ax.plot(kR1[:, 0], kR1[:, 1], label=labels[0])
|
|
|
ax.plot(kR3[:, 0], kR3[:, 1], label=labels[1])
|
|
|
ax.plot(kR10[:, 0], kR10[:, 1], label=labels[2])
|
|
|
- PathEnergyPlot.plot_style(fig, ax)
|
|
|
+ QuadPath.plot_style(fig, ax)
|
|
|
if save_as is not None:
|
|
|
plt.savefig(save_as, dpi=600)
|
|
|
plt.show()
|
|
@@ -176,7 +201,7 @@ def IC_abar_dependence(config_data: dict, save_as=None):
|
|
|
ax.plot(a03[:, 0], a03[:, 1], label=labels[0])
|
|
|
ax.plot(a04[:, 0], a04[:, 1], label=labels[1])
|
|
|
ax.plot(a05[:, 0], a05[:, 1], label=labels[2])
|
|
|
- PathEnergyPlot.plot_style(fig, ax)
|
|
|
+ QuadPath.plot_style(fig, ax)
|
|
|
if save_as is not None:
|
|
|
plt.savefig(save_as, dpi=600)
|
|
|
plt.show()
|
|
@@ -188,13 +213,13 @@ def IC_sigma0_dependence(config_data: dict, save_as=None):
|
|
|
neutral = np.load(em_data_path.joinpath("ZC_-560").joinpath("pathway.npz"))['arr_0']
|
|
|
overchargerd = np.load(em_data_path.joinpath("ZC_-842.74").joinpath("pathway.npz"))['arr_0']
|
|
|
|
|
|
- labels = [rf'$\tilde \sigma_0={s0}$' for s0 in [-0.001, 0, 0.001]]
|
|
|
+ labels = [rf'$\eta={eta}$' for eta in [-0.1, 0, 0.1]]
|
|
|
|
|
|
fig, ax = plt.subplots()
|
|
|
ax.plot(overchargerd[:, 0], overchargerd[:, 1], label=labels[0])
|
|
|
ax.plot(neutral[:, 0], neutral[:, 1], label=labels[1])
|
|
|
ax.plot(undercharged[:, 0], undercharged[:, 1], label=labels[2])
|
|
|
- PathEnergyPlot.plot_style(fig, ax)
|
|
|
+ QuadPath.plot_style(fig, ax)
|
|
|
if save_as is not None:
|
|
|
plt.savefig(save_as, dpi=600)
|
|
|
plt.show()
|
|
@@ -205,8 +230,8 @@ def main():
|
|
|
with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
|
|
|
config_data = json.load(config_file)
|
|
|
|
|
|
- # model_comparison(
|
|
|
- # # save_as=Path(config_data["figures"]).joinpath("Emanuele_data").joinpath('IC_CS_quadrupole_path.pdf')
|
|
|
+ # model_comparison(config_data, save_data=False,
|
|
|
+ # save_as=Path(config_data["figures"]).joinpath("Emanuele_data").joinpath('IC_CS_quadrupole_path.pdf')
|
|
|
# )
|
|
|
|
|
|
# kappaR_dependence(np.array([1, 3, 10]), 0.5,
|
|
@@ -214,25 +239,26 @@ def main():
|
|
|
# )
|
|
|
#
|
|
|
# abar_dependence(np.array([0.3, 0.4, 0.5]), 3,
|
|
|
- # save_as=Path("/home/andraz/ChargedShells/Figures/quadrupole_abar_dep.png"))
|
|
|
+ # save_as=Path("/home/andraz/ChargedShells/Figures/quadrupole_abar_dep.png")
|
|
|
+ # )
|
|
|
|
|
|
- sigma0_dependence(np.array([0.001, 0.002, 0.003]), 3, 0.5,
|
|
|
- # save_as=Path("/home/andraz/ChargedShells/Figures/quadrupole_charge_dep_abar05_kappaR3.png")
|
|
|
- )
|
|
|
+ # sigma0_dependence(np.array([-0.001, 0.00, 0.001]), 3, 0.5,
|
|
|
+ # save_as=Path("/home/andraz/ChargedShells/Figures/quadrupole_charge_dep_abar05_kappaR3.png")
|
|
|
+ # )
|
|
|
|
|
|
# distance_dependence(dist=np.array([2, 3, 4, 6, 10, 20]), kappaR=3, abar=0.5,
|
|
|
# # save_as=Path(config_data["figures"]).joinpath('quadrupole_distance_dep.png')
|
|
|
# )
|
|
|
|
|
|
# IC_kappaR_dependence(config_data,
|
|
|
- # # save_as=Path(config_data["figures"]).joinpath("Emanuele_data"). joinpath('quadrupole_kappaR_dep.png')
|
|
|
+ # save_as=Path(config_data["figures"]).joinpath("Emanuele_data").joinpath('quadrupole_kappaR_dep.png')
|
|
|
# )
|
|
|
#
|
|
|
# IC_abar_dependence(config_data, save_as=Path(config_data["figures"]).joinpath("Emanuele_data").
|
|
|
# joinpath('quadrupole_abar_dep.png'))
|
|
|
#
|
|
|
- # IC_sigma0_dependence(config_data, save_as=Path(config_data["figures"]).joinpath("Emanuele_data").
|
|
|
- # joinpath('quadrupole_charge_dep_abar05_kappaR3.png'))
|
|
|
+ IC_sigma0_dependence(config_data, save_as=Path(config_data["figures"]).joinpath("Emanuele_data").
|
|
|
+ joinpath('quadrupole_charge_dep_abar05_kappaR3.png'))
|
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|