2024-08-30

analysis/ICi_ICp_term_comparison.py

@@ -0,0 +1,135 @@
+import numpy as np
+from matplotlib import gridspec
+from matplotlib.lines import Line2D
+
+from charged_shells import functions as fn
+import matplotlib.pyplot as plt
+import plot_settings
+
+
+def imp_exp_term(l: int, x: float, abar: float):
+    return np.sqrt(2 * l + 1) * fn.sph_bessel_k(l, x) / fn.sph_bessel_k(l+1, x) * abar ** l
+
+
+def perm_exp_term(l: int, x: float, abar: float):
+    return np.sqrt(2 * l + 1) * x ** 2 * fn.sph_bessel_i(l, abar * x) * fn.sph_bessel_k(l, x)
+
+
+def plot_exp(kappaR: list, abar: float, save_as = None):
+    l = np.arange(11)
+    fig = plt.figure(figsize=(3, 1.8))
+    gs = gridspec.GridSpec(1, 1, figure=fig)
+    gs.update(left=0.15, right=0.95, top=0.97, bottom=0.21)
+    ax = fig.add_subplot(gs[0, 0])
+    for kr in kappaR:
+        color = next(plot_settings.COLORS)
+        ax.plot(l, imp_exp_term(l, kr, abar) / imp_exp_term(0, kr, abar), label='ICi', marker='o', ls=':', c=color)
+        ax.plot(l, perm_exp_term(l, kr, abar) / perm_exp_term(0, kr, abar), label='ICp', marker='^', ls=':', c=color)
+    # plt.legend(fontsize=10)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
+    ax.set_xlabel(r'$l$', fontsize=11)
+    ax.set_ylabel('normalized coef', fontsize=11)
+    if save_as is not None:
+        plt.savefig(save_as, dpi=300)
+    plt.show()
+
+
+def plot_exp2(kappaR: float, abar: list, save_as=None, legend=True):
+    l = np.arange(11)
+
+    line1 = Line2D([0], [0], color='black', marker='o', linestyle='None', label='$b_{i,l0}$')
+    line2 = Line2D([0], [0], color='black', marker='^', linestyle='None', label='$c_{i,l0}$')
+
+    abar_lines = []
+
+    fig = plt.figure(figsize=(3, 2))
+    gs = gridspec.GridSpec(1, 1, figure=fig)
+    gs.update(left=0.15, right=0.95, top=0.87, bottom=0.21)
+    ax = fig.add_subplot(gs[0, 0])
+    for ab in abar:
+        color = next(plot_settings.COLORS)
+        ax.plot(l, imp_exp_term(l, kappaR, ab) / imp_exp_term(0, kappaR, ab), label='rr', marker='o', ls=':', c=color)
+        ax.plot(l, perm_exp_term(l, kappaR, ab) / perm_exp_term(0, kappaR, ab), label='tt', marker='^', ls=':',
+                c=color)
+        abar_lines.append(Line2D([0], [0], color=color, linewidth=1.2, ls=':', label=rf'$\bar a={ab}$'))
+
+    # main_legend = ax.legend(frameon=False)
+    if legend:
+        extra_legend = ax.legend(handles=[line1, line2] + abar_lines, loc='upper right', fontsize=11, frameon=False)
+    # extra_legend2 = ax.legend(handles=abar_lines, loc='center right', fontsize=11, frameon=False)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
+    ax.set_xlabel(r'$l$', fontsize=11)
+    ax.set_title(f'$\kappa R={kappaR}$', fontsize=11)
+    ax.set_ylabel('normalized coef', fontsize=11)
+    if save_as is not None:
+        plt.savefig(save_as, dpi=300)
+    plt.show()
+
+
+def plot_imp_ex(abar: float):
+    kappaR = np.geomspace(0.1, 50, 1000)
+    l = np.arange(10)
+    fig, ax = plt.subplots()
+    ax.plot(kappaR, imp_exp_term(l[None, :], kappaR[:, None], abar))
+    ax.set_xscale('log')
+    plt.show()
+
+
+def plot_imp_ex_ratio(abar: float):
+    kappaR = np.geomspace(0.1, 50, 1000)
+    l = np.arange(10)
+    fig, ax = plt.subplots()
+    ax.plot(kappaR, imp_exp_term(l[None, :], kappaR[:, None], abar) / imp_exp_term(l[None, :] + 1, kappaR[:, None], abar))
+    ax.set_xscale('log')
+    plt.show()
+
+
+def plot_perm_ex(abar: float):
+    kappaR = np.geomspace(0.1, 50, 1000)
+    l = np.arange(10)
+    fig, ax = plt.subplots()
+    ax.plot(kappaR, perm_exp_term(l[None, :], kappaR[:, None], abar))
+    ax.set_xscale('log')
+    plt.show()
+
+def plot_perm_ex_ratio(abar: float):
+    kappaR = np.geomspace(0.1, 50, 1000)
+    l = np.arange(10)
+    fig, ax = plt.subplots()
+    ax.plot(kappaR, perm_exp_term(l[None, :], kappaR[:, None], abar) / perm_exp_term(l[None, :] + 1, kappaR[:, None], abar))
+    ax.set_xscale('log')
+    plt.show()
+
+
+def plot_kappaR_dep(l: int, abar: float):
+    kappaR = np.geomspace(0.1, 50, 1000)
+    fig, ax = plt.subplots()
+    ax.plot(kappaR, imp_exp_term(l, kappaR, abar))
+    ax.plot(kappaR, perm_exp_term(l, kappaR, abar))
+    ax.set_xscale('log')
+    plt.show()
+
+
+def plot_kappaR_dep_ratio(l: int, abar: float, l_diff: int = 2):
+    kappaR = np.geomspace(0.1, 50, 1000)
+    fig, ax = plt.subplots()
+    ax.plot(kappaR, imp_exp_term(l + l_diff, kappaR, abar) / imp_exp_term(l, kappaR, abar))
+    ax.plot(kappaR, perm_exp_term(l + l_diff, kappaR, abar) / perm_exp_term(l, kappaR, abar))
+    ax.set_xscale('log')
+    plt.show()
+
+
+if __name__ == '__main__':
+
+    # plot_exp([1, 10, 50], 0.4,
+    #          # save_as='/home/andraz/ChargedShells/Figures/ici_icp_expansion_a08_kr10.png'
+    #          )
+
+    plot_exp2(1, [0.5, 0.8], legend=True,
+             save_as='/home/andraz/ChargedShells/Figures/ici_icp_expansion_kr1.png'
+             )
+
+    # plot_kappaR_dep_ratio(2, 0.8, l_diff=2)
+
+    # plot_imp_ex_ratio(0.5)
+    # plot_perm_ex_ratio(0.5)

analysis/dip_path_plot.py

@@ -1,5 +1,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
+from matplotlib import gridspec
+
 from charged_shells.rotational_path import PairRotationalPath, PathEnergyPlot
 from charged_shells import expansion, interactions
 from charged_shells.parameters import ModelParams
@@ -43,10 +45,10 @@ DipolePath3.set_default_x_axis(zero_to_pi_half)
 DipolePath3.add_euler(beta2=np.pi/2, beta1=zero_to_pi_half, start_name="EP", end_name="EE")
 DipolePath3.add_euler(beta2=pi_half_to_pi, beta1=pi_half_to_pi, end_name="PP")
 DipolePath3.add_euler(beta2=pi_half_to_pi[::-1], beta1=np.pi, end_name="EP")
-DipolePath3.add_euler(beta2=pi_half_to_pi, beta1=pi_half_to_pi[::-1], end_name="EP")
-DipolePath3.add_euler(beta1=np.pi/2, beta2=pi_half_to_pi[::-1], alpha2=np.pi/2, end_name="tEE")
-DipolePath3.add_euler(beta1=np.pi/2, beta2=np.pi/2, alpha1=zero_to_pi_half[::-1], end_name="EE")
-DipolePath3.add_euler(beta1=np.pi/2, beta2=zero_to_pi_half[::-1], end_name="EP")
+DipolePath3.add_euler(beta2=pi_half_to_pi, beta1=pi_to_three_halves_pi, end_name="EP")
+DipolePath3.add_euler(beta1=3 * np.pi/2, beta2=pi_half_to_pi[::-1], alpha2=np.pi/2, end_name="tEE")
+DipolePath3.add_euler(beta1=3 * np.pi/2, beta2=np.pi/2, alpha1=zero_to_pi_half[::-1], end_name="EE")
+DipolePath3.add_euler(beta1=3 * np.pi/2, beta2=pi_half_to_pi, end_name="EP")
 
 
 def sections_plot(kappaR: float = 3, abar: float = 0.5, sigma_tilde=0.001, save_as=None):
@@ -287,10 +289,10 @@ def combined_sigma0_dependence(config_data: dict, kappaR=3., abar=0.5, sigma0=(-
 
 def combined_all(config_data: dict, save_as=None):
 
-    sigma_tilde = 0.001
+    sigma_tilde = 0.00099
     kappaR_list = [1, 3, 10]
     abar_list = [0.5, 0.4, 0.3]
-    sigma0_list = [-0.0002, 0.00, 0.0002]
+    sigma0_list = [-0.000198, 0.00, 0.000198]
     kappaR = 3
     abar = 0.5
 
@@ -309,7 +311,7 @@ def combined_all(config_data: dict, save_as=None):
 
     ex1 = expansion.MappedExpansionDipole(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
-    ex3 = ex1.clone().inverse_sign()
+    ex3 = ex1.clone().inverse_sign(exclude_00=True)
 
     path_plot = PathEnergyPlot(ex1, ex2, DipolePath3, dist=2., params=params, match_expansion_axis_to_params=0)
     energy_kappaR = path_plot.evaluate_path()
@@ -334,7 +336,7 @@ def combined_all(config_data: dict, save_as=None):
 
     ex1 = expansion.MappedExpansionDipole(np.asarray(abar_list), params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
-    ex3 = ex1.clone().inverse_sign()
+    ex3 = ex1.clone().inverse_sign(exclude_00=True)
 
     path_plot = PathEnergyPlot(ex1, ex2, DipolePath3, dist=2., params=params, match_expansion_axis_to_params=None)
     energy_abar = path_plot.evaluate_path()
@@ -355,48 +357,55 @@ def combined_all(config_data: dict, save_as=None):
     overcharged, overcharged_inv = overcharged[:int(len(overcharged) / 2)], overcharged[int(len(overcharged) / 2):]
     neutral, neutral_inv = neutral[:int(len(neutral) / 2)], neutral[int(len(neutral) / 2):]
 
-    ic_data_sigma0 = [overcharged, neutral, undercharged]
-    ic_data_sigma0_inv = [overcharged_inv, neutral_inv, undercharged_inv]
+    ic_data_sigma0 = [undercharged, neutral, overcharged]
+    ic_data_sigma0_inv = [undercharged_inv, neutral_inv, overcharged_inv]
 
     params = ModelParams(R=150, kappaR=kappaR)
     ex1 = expansion.MappedExpansionDipole(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=np.asarray(sigma0_list))
     ex2 = ex1.clone()
-    ex3 = ex1.clone().inverse_sign()
+    ex3 = ex1.clone().inverse_sign(exclude_00=True)
 
     path_plot = PathEnergyPlot(ex1, ex2, DipolePath3, dist=2., params=params, match_expansion_axis_to_params=None)
     energy_sigma0 = path_plot.evaluate_path()
     path_plot_inv = PathEnergyPlot(ex1, ex3, DipolePath3, dist=2., params=params, match_expansion_axis_to_params=None)
     energy_sigma0_inv = path_plot_inv.evaluate_path()
     x_axis_sigma0 = path_plot.rot_path.stack_x_axes()
-    labels_sigma0 = [rf'$\eta={s0/sigma_tilde}$' for s0 in sigma0_list]
-
-    fig, axs = plt.subplots(3, 1, figsize=(6, 7.8))
+    labels_sigma0 = [rf'$\eta={s0/sigma_tilde:.1f}$' for s0 in sigma0_list]
+
+    # fig, axs = plt.subplots(3, 1, figsize=(6, 7.8))
+    fig = plt.figure(figsize=(4, 3.6))
+    gs = gridspec.GridSpec(2, 1, figure=fig)
+    # gs.update(left=0.12, right=0.975, top=0.96, bottom=0.04, hspace=0.3)
+    gs.update(left=0.12, right=0.975, top=0.94, bottom=0.06, hspace=0.3)
+    # axs = [fig.add_subplot(gs[0, 0]), fig.add_subplot(gs[1, 0]), fig.add_subplot(gs[2, 0])]
+    axs = [fig.add_subplot(gs[0, 0]), fig.add_subplot(gs[1, 0])]
     for d, d_inv, en, en_inv, label, c in zip(ic_data_kappaR, ic_data_kappaR_inv, energy_kappaR.T, energy_kappaR_inv.T, labels_kappaR, COLOR_LIST):
-        axs[0].set_title('Screening', fontsize=15)
+        axs[0].set_title('Screening', fontsize=11, y=0.98)
         axs[0].plot(d[:, 0], d[:, 1], label=label, c=c)
         axs[0].plot(x_axis_kappaR, en, ls='--', c=c)
         axs[0].plot(d_inv[:, 0], d_inv[:, 1], c=c)
         axs[0].plot(x_axis_kappaR, en_inv, ls='--', c=c)
         DipolePath3.plot_style(fig, axs[0], size=None)
-        axs[0].get_legend().set_bbox_to_anchor((0.65, 1.07))
-    for d, d_inv, en, en_inv, label, c in zip(ic_data_abar, ic_data_abar_inv, energy_abar.T, energy_abar_inv.T, labels_abar, COLOR_LIST):
-        axs[1].set_title('Asymmetry', fontsize=15)
+        axs[0].get_legend().set_bbox_to_anchor((0.65, 1.03))
+    # for d, d_inv, en, en_inv, label, c in zip(ic_data_abar, ic_data_abar_inv, energy_abar.T, energy_abar_inv.T, labels_abar, COLOR_LIST):
+    #     axs[1].set_title('Eccentricity', fontsize=11, y=0.98)
+    #     axs[1].plot(d[:, 0], d[:, 1], label=label, c=c)
+    #     axs[1].plot(x_axis_abar, en, ls='--', c=c)
+    #     axs[1].plot(d_inv[:, 0], d_inv[:, 1], c=c)
+    #     axs[1].plot(x_axis_abar, en_inv, ls='--', c=c)
+    #     DipolePath3.plot_style(fig, axs[1], size=None)
+    #     axs[1].get_legend().set_bbox_to_anchor((0.65, 1.02))
+    for d, d_inv, en, en_inv, label, c in reversed(list(zip(ic_data_sigma0, ic_data_sigma0_inv, energy_sigma0.T,
+                                                            energy_sigma0_inv.T, labels_sigma0, COLOR_LIST[:3][::-1]))):
+        axs[1].set_title('Net charge', fontsize=11, y=0.98)
         axs[1].plot(d[:, 0], d[:, 1], label=label, c=c)
-        axs[1].plot(x_axis_abar, en, ls='--', c=c)
+        axs[1].plot(x_axis_sigma0, en, ls='--', c=c)
         axs[1].plot(d_inv[:, 0], d_inv[:, 1], c=c)
-        axs[1].plot(x_axis_abar, en_inv, ls='--', c=c)
+        axs[1].plot(x_axis_sigma0, en_inv, ls='--', c=c)
         DipolePath3.plot_style(fig, axs[1], size=None)
-        axs[1].get_legend().set_bbox_to_anchor((0.65, 1.05))
-    for d, d_inv, en, en_inv, label, c in zip(ic_data_sigma0, ic_data_sigma0_inv, energy_sigma0.T, energy_sigma0_inv.T, labels_sigma0, COLOR_LIST):
-        axs[2].set_title('Net charge', fontsize=15)
-        axs[2].plot(d[:, 0], d[:, 1], label=label, c=c)
-        axs[2].plot(x_axis_sigma0, en, ls='--', c=c)
-        axs[2].plot(d_inv[:, 0], d_inv[:, 1], c=c)
-        axs[2].plot(x_axis_sigma0, en_inv, ls='--', c=c)
-        DipolePath3.plot_style(fig, axs[2], size=None)
-        axs[2].get_legend().set_bbox_to_anchor((0.65, 1.04))
+        axs[1].get_legend().set_bbox_to_anchor((0.65, 1.02))
     for ax in axs:
-        ax.yaxis.set_label_coords(-0.09, 0.5)
+        ax.yaxis.set_label_coords(-0.08, 0.5)
     # axs[-1].set_xlabel('rotational path', fontsize=15)
     plt.tight_layout()
     if save_as is not None:
@@ -437,7 +446,7 @@ if __name__ == '__main__':
     #
     # combined_sigma0_dependence(config_data,
     #                      # save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('janus_charge_dep.png')
-    #                      )
+    #                     )
     #
     combined_all(config_data,
                  save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('janus_combined_dep.png')

analysis/higher_multipole_matching.py

@@ -0,0 +1,127 @@
+from analysis.peak_heigth import *
+
+Array = np.ndarray
+Expansion = expansion.Expansion
+RotConfig = Literal['ep', 'pp', 'sep']
+
+
+def get_kappaR_higher_multipole_dicts(peak: Peak, params: ModelParams, emanuele_data: Array,
+                            kappaR: Array, dist_cs: list, max_kappaR: float = 50, min_kappa: float = 0.01,
+                                      only_loaded: bool = False):
+    distances, inverse = np.unique(emanuele_data[:, 0], return_inverse=True)
+
+    energy_fn = mapping.parameter_map_two_expansions(interactions.charged_shell_energy,
+                                                     peak.kappaR_axis_in_expansion)
+
+    energy = []
+    for d in dist_cs:
+        if not only_loaded:
+            energy.append(energy_fn(peak.ex1, peak.ex2, params, dist=d))
+        else:
+            energy.append(np.full((len(kappaR),), 1.))
+
+    data_dict_ic = {}
+    data_dict_cs = {}
+    k = 0
+    for i in range(len(distances)):
+        d = np.around(distances[i], 5)
+        if d in np.around(dist_cs, 5):
+            idx, = np.nonzero(inverse == i)
+            kR = emanuele_data[idx, 1][(emanuele_data[idx, 1] <= max_kappaR) * (min_kappa <= emanuele_data[idx, 1])]
+            en = emanuele_data[idx, peak.emanuele_data_column-1][(emanuele_data[idx, 1] <= max_kappaR) * (min_kappa <= emanuele_data[idx, 1])]
+            sort = np.argsort(kR)
+            if peak.log_y:
+                data_dict_ic[d] = np.stack((kR[sort], np.abs(en)[sort])).T
+                data_dict_cs[d] = np.stack((kappaR, np.abs(energy[k]))).T
+            else:
+                data_dict_ic[d] = np.stack((kR[sort], en[sort])).T
+                data_dict_cs[d] = np.stack((kappaR, energy[k])).T
+            k += 1
+
+    return data_dict_ic, data_dict_cs
+
+
+def IC_peak_energy_plot(config_data: dict,
+                        dist: list,
+                        which: RotConfig,
+                        a_bar=0.8,
+                        min_kappaR: float = 0.01,
+                        max_kappaR: float = 50.,
+                        R: float = 150,
+                        save_as: Path = None,
+                        save_data: bool = False,
+                        quad_correction: bool = False,
+                        log_y: bool = True):
+
+    # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_11"))
+    em_data_path = (Path(config_data["emanuele_data"]).joinpath("TEST_HIGHER_MULTIPOLES"))
+    em_data = np.load(em_data_path.joinpath("higher_multiplot_energy.npz"))
+    data2 = em_data['L2']
+    data4 = em_data['L4']
+    data6 = em_data['L6']
+
+
+    num = 100 if which == 'sep' else 30
+    kappaR = np.geomspace(min_kappaR, max_kappaR, num)
+    params = ModelParams(R=R, kappaR=kappaR)
+
+    ex = expansion.MappedExpansionQuad(a_bar, params.kappaR, 0.00099, l_max=20)
+
+    if which == 'ep':
+        peak = PeakEP(ex, log_y, kappaR_axis_in_expansion=0)
+    elif which == 'pp':
+        peak = PeakPP(ex, log_y, kappaR_axis_in_expansion=0)
+    elif which == 'sep':
+        peak = PeakSEP(ex, log_y, kappaR_axis_in_expansion=0)
+    else:
+        raise ValueError
+
+    data_dict_ic2, data_dict_cs = get_kappaR_higher_multipole_dicts(peak, params, data2, kappaR, dist, max_kappaR, min_kappaR)
+    data_dict_ic4, _ = get_kappaR_higher_multipole_dicts(peak, params, data4, kappaR, dist, max_kappaR, min_kappaR)
+    data_dict_ic6, _ = get_kappaR_higher_multipole_dicts(peak, params, data6, kappaR, dist, max_kappaR, min_kappaR)
+
+    colors = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
+
+    fig = plt.figure(figsize=(2, 1.8))
+    gs = gridspec.GridSpec(1, 1, figure=fig)
+    gs.update(left=0.25, right=0.95, top=0.97, bottom=0.21)
+    ax = fig.add_subplot(gs[0, 0])
+    for (key, data2), data_cs, data4, data6 in zip(data_dict_ic2.items(), data_dict_cs.values(), data_dict_ic4.values(), data_dict_ic6.values()):
+        current_color = next(colors)
+        ax.plot(data2[:, 0], data2[:, 1], label=r'$l=2$', c=current_color)
+        ax.plot(data_cs[:, 0], data_cs[:, 1], ls='--', c='k')
+        ax.plot(data4[:, 0], data4[:, 1], ls='-', c=next(colors), label=r'$l=4$')
+        ax.plot(data6[:, 0], data6[:, 1], ls='-', c=next(colors), label=r'$l=6$')
+    ax.legend(fontsize=9, ncol=1, frameon=False, handlelength=0.7, loc='upper right',
+              bbox_to_anchor=(0.6, 0.5),
+              # bbox_to_anchor=(0.42, 0.9)
+              )
+    ax.set_title(rf'$\rho = {key:.1f}$', loc='right',
+                 # x=0.95, y=0.8,
+                 x=0.45, y=0.6,
+                 )
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
+    ax.set_xlabel(r'$\kappa R$', fontsize=11)
+
+    ax.set_ylabel(peak.y_label, fontsize=11)
+    if log_y:
+        ax.set_yscale('log')
+    ax.set_xscale('log')
+    # plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=300)
+    plt.show()
+
+
+if __name__ == '__main__':
+    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
+        config_data = json.load(config_file)
+
+    dist = [2.5]
+    IC_peak_energy_plot(config_data, dist=dist, which='ep', min_kappaR=0.1, max_kappaR=50,
+                        save_as=Path(
+                            '/home/andraz/ChargedShells/Figures/Emanuele_data/peak_ep_kappaR_higher_correction_rho25.png'),
+                        save_data=False,
+                        quad_correction=False,
+                        log_y=False
+                        )

analysis/patch_shape_comparison.py

@@ -219,33 +219,35 @@ def ic_cs_comparison2(save_data=False):
                  CSp_eta_0=np.stack((theta, potential_cs[1])).T,
                  CSp_eta_pos=np.stack((theta, potential_cs[2])).T)
 
-    fig, ax = plt.subplots(figsize=(3, 2.5))
+    fig, ax = plt.subplots(figsize=(2, 1.7))
     lines = []
     # lines.append(plt.scatter([], [], marker='x', c='k', label='CSp'))
     # lines.append(plt.scatter([], [], marker='o', facecolor='none', edgecolor='k', label='ICi'))
     for p_ic, p_cs, charge in zip(potential_ic, potential_cs, sigma0_array):
-        l, = ax.plot(theta, 1000 * p_cs.T, label=rf'$\eta = {charge/sigma_m}$', linewidth=2, zorder=0)
-        l2, = ax.plot(theta, 1000 * p_cs.T, linewidth=2, zorder=0, c='k', ls='--', dashes=(5, 5))
+        l, = ax.plot(theta, 1000 * p_cs.T, label=rf'$\eta = {charge/sigma_m}$', linewidth=1.5, zorder=0)
+        l2, = ax.plot(theta, 1000 * p_cs.T, linewidth=1.5, zorder=0, c='k', ls='--', dashes=(5, 5))
         lines.append(l)
         # ax.scatter(theta[::100], 1000 * p_ic.T[::100], marker='o', facecolors='none', edgecolors='k', s=50)
         # ax.scatter(theta[::100], 1000 * p_cs.T[::100], marker='x', c='k', s=50)
     # ax.plot(theta, potential_ic.T, label='IC', ls=':', linewidth=2, marker='o', markevery=50, mfc='none')
-    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
-    ax.set_xlabel(r'$\theta$', fontsize=15)
-    ax.set_ylabel(r'$\phi [mV]$', fontsize=15)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
+    ax.set_xlabel(r'$\theta$', fontsize=11)
+    ax.set_ylabel(r'$\phi [mV]$', fontsize=11)
     custom_ticks = [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4, np.pi]
     custom_labels = ['$0$', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$', r'$\pi$']
     plt.axhline(y=0, color='black', linestyle=':')
     plt.axvline(x=target_patch_size, color='black', linestyle=':')
-    plt.xticks(custom_ticks, custom_labels, fontsize=15)
+    plt.xticks(custom_ticks, custom_labels, fontsize=11)
 
     # lines.append(Line2D([0], [0], color='k', label='CSp'))
     # lines.append(plt.scatter([], [], marker='o', facecolor='none', edgecolor='red', label='ICi'))
-    plt.legend(handles=lines, fontsize=14, handlelength=0.8, frameon=False,
+    plt.legend(handles=lines, fontsize=9, handlelength=0.8, frameon=False,
                bbox_to_anchor=(0.57, 1.0), loc='upper center'
                )
 
-    plt.tight_layout()
+    # plt.tight_layout()
+    plt.subplots_adjust(left=0.22, right=0.97, top=0.95, bottom=0.22)
+    ax.xaxis.set_label_coords(0.5, -0.17)
     plt.savefig(Path("/home/andraz/ChargedShells/Figures/final_figures/potential_ic_cs_comparison.png"), dpi=300)
     plt.show()
 

analysis/patch_size_dependence.py

@@ -30,18 +30,20 @@ def plot_abar_dependence(*, save=False, save_data=False):
         np.savez(Path(config_data["figure_data"]).joinpath("fig_6.npz"), **data_dict)
 
     # fig, ax = plt.subplots(figsize=plt.figaspect(0.5)/1.5)
-    fig, ax = plt.subplots(figsize=(3, 2.5))
+    fig, ax = plt.subplots(figsize=(2, 1.7))
     for patch, kR, marker in zip(ps.T, kappaR, markers):
-        ax.plot(a_bar, patch, label=rf'$\kappa R$ = {kR}', linewidth=2.5,
+        ax.plot(a_bar, patch, label=rf'$\kappa R$ = {kR}', linewidth=1.5,
                 # marker=marker, markerfacecolor='none', markersize=8, markevery=10
                 )
-    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
-    ax.set_xlabel(r'$\bar a$', fontsize=15)
-    ax.set_ylabel(r'$\theta_0$', fontsize=15)
-    plt.legend(handlelength=0.6, fontsize=13, frameon=False, loc='upper right', bbox_to_anchor=(0.56, 0.58))
-    plt.tight_layout()
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
+    ax.set_xlabel(r'$\bar a$', fontsize=11)
+    ax.set_ylabel(r'$\theta_0$', fontsize=11)
+    plt.legend(handlelength=0.6, fontsize=9, frameon=False, loc='upper right', bbox_to_anchor=(0.55, 0.63))
+    # plt.tight_layout()
+    plt.subplots_adjust(left=0.22, right=0.97, top=0.95, bottom=0.22)
     plt.axhline(y=mark_ps, color='black', linestyle=':')
     plt.axvline(x=0.5, color='black', linestyle=':')
+    ax.xaxis.set_label_coords(0.5, -0.17)
     if save:
         plt.savefig(Path("/home/andraz/ChargedShells/Figures/final_figures/patch_size_potential.png"), dpi=300)
     plt.show()

analysis/peak_heigth.py

@@ -1,3 +1,4 @@
+from matplotlib import gridspec
 from charged_shells import expansion, interactions, mapping, functions
 from charged_shells.parameters import ModelParams
 import numpy as np
@@ -28,7 +29,7 @@ class PeakEP(Peak):
 
     def __init__(self, ex: Expansion, log_y: bool = False, kappaR_axis_in_expansion: int = None):
         self.emanuele_data_column = 4
-        self.y_label = r'$|U_{EP}|$' if log_y else r'$U_{EP}$'
+        self.y_label = r'$|V_\mathrm{EP}|$' if log_y else r'$V_\mathrm{EP}$'
         self.ex1 = ex.clone()
         self.ex2 = ex.clone()
         self.ex1.rotate_euler(alpha=0, beta=np.pi / 2, gamma=0)
@@ -40,7 +41,7 @@ class PeakPP(Peak):
 
     def __init__(self, ex: Expansion, log_y: bool = False, kappaR_axis_in_expansion: int = None):
         self.emanuele_data_column = 3
-        self.y_label = r'$U_{PP}$'
+        self.y_label = r'$V_\mathrm{PP}$'
         self.ex1 = ex.clone()
         self.ex2 = ex.clone()
         self.log_y = log_y
@@ -51,11 +52,11 @@ class PeakSEP(Peak):
 
     def __init__(self, ex: Expansion, log_y: bool = False, kappaR_axis_in_expansion: int = None):
         self.emanuele_data_column = 5
-        self.y_label = r'$|U_{sEP}|$' if log_y else r'$U_{sEP}$'
+        self.y_label = r'$|V_\mathrm{sEP}|$' if log_y else r'$V_\mathrm{sEP}$'
         self.ex1 = ex.clone()
         self.ex2 = ex.clone()
         self.ex1.rotate_euler(alpha=0, beta=np.pi / 4, gamma=0)
-        self.ex2.rotate_euler(alpha=0, beta=np.pi / 4, gamma=0)
+        self.ex2.rotate_euler(alpha=0, beta=3 * np.pi / 4, gamma=0)
         self.log_y = log_y
         self.kappaR_axis_in_expansion = kappaR_axis_in_expansion
 
@@ -88,12 +89,16 @@ def get_charge_energy_dicts(peak: Peak, params: ModelParams, emanuele_data: Arra
     return data_dict_ic, data_dict_cs
 
 
-def get_kappaR_energy_dicts(peak: Peak, params: ModelParams, emanuele_data: Array, kappaR: Array, abar_cs: list, max_kappaR: float = 50):
+def get_kappaR_energy_dicts(peak: Peak, params: ModelParams, emanuele_data: Array,
+                            kappaR: Array, abar_cs: list, max_kappaR: float = 50, only_loaded: bool = False):
     abar_ic, inverse = np.unique(emanuele_data[:, 1], return_inverse=True)
 
     energy_fn = mapping.parameter_map_two_expansions(partial(interactions.charged_shell_energy, dist=2),
                                                      peak.kappaR_axis_in_expansion)
-    energy = energy_fn(peak.ex1, peak.ex2, params)
+    if not only_loaded:
+        energy = energy_fn(peak.ex1, peak.ex2, params)
+    else:
+        energy = np.full((len(abar_cs), len(kappaR)), 1.)
 
     data_dict_ic = {}
     data_dict_cs = {}
@@ -177,7 +182,7 @@ def IC_peak_energy_plot(config_data: dict,
     kappaR = np.geomspace(min_kappaR, max_kappaR, num)
     params = ModelParams(R=R, kappaR=kappaR)
 
-    ex = expansion.MappedExpansionQuad(np.array(a_bar)[:, None], params.kappaR[None, :], 0.001, l_max=20)
+    ex = expansion.MappedExpansionQuad(np.array(a_bar)[:, None], params.kappaR[None, :], 0.00099, l_max=20)
 
     if which == 'ep':
         peak = PeakEP(ex, log_y, kappaR_axis_in_expansion=1)
@@ -190,27 +195,43 @@ def IC_peak_energy_plot(config_data: dict,
 
     data_dict_ic, data_dict_cs = get_kappaR_energy_dicts(peak, params, data, kappaR, a_bar, max_kappaR)
 
+    corrected_data = {}
+    for key, data in data_dict_ic.items():
+        corrected_data[key] = data[:, 1] * higher_multipole_ic_correction(data[:, 0], float(key), higher_l=4)
+
+    corrected_data2 = {}
+    for key, data in data_dict_ic.items():
+        corrected_data2[key] = data[:, 1] * higher_multipole_ic_correction(data[:, 0], float(key), higher_l=6)
+
     colors = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
 
-    fig, ax = plt.subplots(figsize=(4.25, 3))
-    for (key, data1), data2 in zip(data_dict_ic.items(), data_dict_cs.values()):
-        current_color = next(colors)
-        ax.plot(data1[:, 0], data1[:, 1], label=rf'$\bar a = {key:.2f}$', c=current_color)
-        ax.plot(data2[:, 0], data2[:, 1], ls='--', c=current_color)
-    ax.legend(fontsize=14, ncol=1, frameon=False, handlelength=0.7, loc='upper right',
-              bbox_to_anchor=(0.42, 0.42),
+    fig = plt.figure(figsize=(2, 1.8))
+    gs = gridspec.GridSpec(1, 1, figure=fig)
+    gs.update(left=0.25, right=0.95, top=0.97, bottom=0.21)
+    ax = fig.add_subplot(gs[0, 0])
+    i = 0
+    for (key, data1), data2, cd, cd2 in zip(data_dict_ic.items(), data_dict_cs.values(), corrected_data.values(), corrected_data2.values()):
+        i += 1
+        if i == 1:
+            current_color = next(colors)
+            ax.plot(data1[:, 0], data1[:, 1], label=rf'$\bar a = {key:.2f}$', c=current_color)
+            ax.plot(data2[:, 0], data2[:, 1], ls='--', c=current_color)
+            ax.plot(data1[:, 0], cd, ls='-', c='r', label=r'$l=4$')
+            ax.plot(data1[:, 0], cd2, ls='-', c='g', label=r'$l=6$')
+    ax.legend(fontsize=9, ncol=1, frameon=False, handlelength=0.7, loc='upper right',
+              bbox_to_anchor=(0.7, 0.52),
               # bbox_to_anchor=(0.42, 0.9)
               )
-    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
-    ax.set_xlabel(r'$\kappa R$', fontsize=15)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
+    ax.set_xlabel(r'$\kappa R$', fontsize=11)
 
-    ax.set_ylabel(peak.y_label, fontsize=15)
+    ax.set_ylabel(peak.y_label, fontsize=11)
     if log_y:
         ax.set_yscale('log')
     ax.set_xscale('log')
-    plt.tight_layout()
+    # plt.tight_layout()
     if save_as is not None:
-        plt.savefig(save_as, dpi=600)
+        plt.savefig(save_as, dpi=300)
     plt.show()
 
 
@@ -347,7 +368,7 @@ def IC_peak_energy_kappaR_combined_plot(config_data: dict,
     params = ModelParams(R=R, kappaR=kappaR)
 
     ex = expansion.MappedExpansionQuad(np.sort(np.array(a_bar))[:, None],  # sorting necessary as it is also in energy_dicts()
-                                       params.kappaR[None, :], 0.001, l_max=20)
+                                       params.kappaR[None, :], 0.00099, l_max=20)
 
     peak_ep = PeakEP(ex, log_y=True, kappaR_axis_in_expansion=1)
     peak_pp = PeakPP(ex, log_y=True, kappaR_axis_in_expansion=1)
@@ -357,35 +378,41 @@ def IC_peak_energy_kappaR_combined_plot(config_data: dict,
     data_ic = []
     data_cs = []
     for peak in peaks:
-        dict_ic, dict_cs = get_kappaR_energy_dicts(peak, params, data, kappaR, a_bar, max_kappaR)
+        dict_ic, dict_cs = get_kappaR_energy_dicts(peak, params, data, kappaR, a_bar, max_kappaR, only_loaded=False)
         data_ic.append(dict_ic)
         data_cs.append(dict_cs)
 
-    legend_coords = [(0.58, 0.45), (0.58, 0.43), (0.58, 0.9)]
+    legend_coords = [(0.55, 0.45), (0.55, 0.43), (0.55, 0.63)]
 
-    fig, axs = plt.subplots(3, 1, figsize=(3, 7.8))
+    # fig, axs = plt.subplots(3, 1, figsize=(2, 5.4))
+    fig = plt.figure(figsize=(2, 5.4))
+    gs = gridspec.GridSpec(3, 1, figure=fig)
+    gs.update(left=0.25, right=0.95, top=0.99, bottom=0.07, hspace=0.3)
+    axs = [fig.add_subplot(gs[0, 0]), fig.add_subplot(gs[1, 0]), fig.add_subplot(gs[2, 0])]
     for ax, data_dict_ic, data_dict_cs, peak, lc in zip(axs, data_ic, data_cs, peaks, legend_coords):
         colors = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
         for ab in a_bar:
             key = np.around(ab, 5)
             current_color = next(colors)
-            ax.plot(data_dict_ic[key][:, 0], data_dict_ic[key][:, 1], label=rf'$\bar a = {key:.2f}$', c=current_color)
-            ax.plot(data_dict_cs[key][:, 0], data_dict_cs[key][:, 1], ls='--', c=current_color)
-        ax.axvline(x=10, c='k', ls=':')
-        ax.axvline(x=3, c='k', ls=':')
-        ax.axvline(x=1, c='k', ls=':')
-        ax.legend(fontsize=13, ncol=1, frameon=False, handlelength=0.7, loc='upper right', bbox_to_anchor=lc)
-        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+            ax.plot(data_dict_ic[key][:, 0], data_dict_ic[key][:, 1], label=rf'$\bar a = {key:.1f}$', c=current_color, linewidth=1.5)
+            ax.plot(data_dict_cs[key][:, 0], data_dict_cs[key][:, 1], ls='--', c=current_color, linewidth=1.5)
+        ax.axvline(x=10, c='k', ls=':', linewidth=1)
+        ax.axvline(x=3, c='k', ls=':', linewidth=1)
+        ax.axvline(x=1, c='k', ls=':', linewidth=1)
+        ax.legend(fontsize=9, ncol=1, frameon=False, handlelength=0.7, loc='upper right', bbox_to_anchor=lc)
+        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
         # if ax == axs[-1]:
-        ax.set_xlabel(r'$\kappa R$', fontsize=15)
+        ax.set_xlabel(r'$\kappa R$', fontsize=11)
 
-        ax.set_ylabel(peak.y_label, fontsize=15)
+        ax.set_ylabel(peak.y_label, fontsize=11)
         if peak.log_y:
             ax.set_yscale('log')
         ax.set_xscale('log')
-        ax.yaxis.set_label_coords(-0.2, 0.5)
-    plt.subplots_adjust(left=0.3)
-    plt.tight_layout()
+        ax.yaxis.set_label_coords(-0.25, 0.5)
+        ax.xaxis.set_label_coords(0.5, -0.12)
+    # plt.subplots_adjust(left=0.1)
+    # plt.tight_layout()
+    # plt.subplots_adjust(left=0.25, right=0.97, top=0.85, bottom=0.15)
     if save_as is not None:
         plt.savefig(save_as, dpi=300)
     plt.show()
@@ -403,8 +430,8 @@ def IC_peak_energy_charge_combined_plot(config_data: dict,
     data = em_data['changeZc']
 
     params = ModelParams(R=R, kappaR=3)
-    sigma0 = np.linspace(-0.0003, 0.0003, 300)
-    sigma_tilde = 0.001
+    sigma0 = np.linspace(-0.000297, 0.000297, 300)
+    sigma_tilde = 0.00099
 
     ex = expansion.MappedExpansionQuad(np.sort(np.array(a_bar)),  # sorting necessary as it is also in energy_dicts()
                                        kappaR=3, sigma0=sigma0, l_max=20, sigma_tilde=sigma_tilde)
@@ -421,64 +448,74 @@ def IC_peak_energy_charge_combined_plot(config_data: dict,
         data_ic.append(dict_ic)
         data_cs.append(dict_cs)
 
-    fig, axs = plt.subplots(3, 1, figsize=(3, 7.8))
+    # fig, axs = plt.subplots(3, 1, figsize=(3, 7.8))
+    fig = plt.figure(figsize=(2, 5.4))
+    gs = gridspec.GridSpec(3, 1, figure=fig)
+    gs.update(left=0.25, right=0.95, top=0.99, bottom=0.07, hspace=0.3)
+    axs = [fig.add_subplot(gs[0, 0]), fig.add_subplot(gs[1, 0]), fig.add_subplot(gs[2, 0])]
     for ax, data_dict_ic, data_dict_cs, peak in zip(axs, data_ic, data_cs, peaks):
         colors = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
         for ab in a_bar:
             key = np.around(ab, 5)
             current_color = next(colors)
-            ax.plot(data_dict_ic[key][:, 0], data_dict_ic[key][:, 1], label=rf'$\bar a = {key:.2f}$', c=current_color)
-            ax.plot(data_dict_cs[key][:, 0], data_dict_cs[key][:, 1], ls='--', c=current_color)
-        ax.legend(fontsize=13, ncol=1, frameon=False, handlelength=0.7, loc='upper right',
+            ax.plot(data_dict_ic[key][:, 0], data_dict_ic[key][:, 1], label=rf'$\bar a = {key:.1f}$', c=current_color, linewidth=1.5)
+            ax.plot(data_dict_cs[key][:, 0], data_dict_cs[key][:, 1], ls='--', c=current_color, linewidth=1.5)
+        ax.legend(fontsize=9, ncol=1, frameon=False, handlelength=0.7, loc='upper right',
                   bbox_to_anchor=(0.77, 1.03),
                   )
-        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
         # if ax == axs[-1]:
-        ax.set_xlabel(r'$\eta$', fontsize=15)
+        ax.set_xlabel(r'$\eta$', fontsize=11)
 
-        ax.set_ylabel(peak.y_label, fontsize=15)
+        ax.set_ylabel(peak.y_label, fontsize=11)
         if peak.log_y:
             ax.set_yscale('log')
         # ax.set_xscale('log')
-        ax.yaxis.set_label_coords(-0.2, 0.5)
-    plt.subplots_adjust(left=0.3)
-    plt.tight_layout()
+        ax.yaxis.set_label_coords(-0.18, 0.5)
+        ax.xaxis.set_label_coords(0.5, -0.12)
+    # plt.subplots_adjust(left=0.3)
+    # plt.tight_layout()
     if save_as is not None:
         plt.savefig(save_as, dpi=300)
     plt.show()
 
 
+def higher_multipole_ic_correction(kappaR, abar, higher_l: int = 4):
+    return (functions.sph_bessel_k(3, kappaR) * functions.sph_bessel_i(2, kappaR * abar) * abar ** (higher_l - 2) /
+            (functions.sph_bessel_k(higher_l + 1, kappaR) * functions.sph_bessel_i(higher_l, kappaR * abar)))
+
+
 def main():
     with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
         config_data = json.load(config_file)
 
-    # a_bar = [0.2, 0.5, 0.8]
-    # IC_peak_energy_plot(config_data, a_bar=a_bar, which='ep', max_kappaR=50,
-    #                     # save_as=Path('/home/andraz/ChargedShells/Figures/Emanuele_data/peak_pp_kappaR.png'),
+    a_bar = [0.2, 0.5, 0.8]
+    # IC_peak_energy_plot(config_data, a_bar=a_bar, which='pp', min_kappaR=0.01, max_kappaR=50,
+    #                     save_as=Path('/home/andraz/ChargedShells/Figures/Emanuele_data/peak_pp_kappaR_higher_correction_abar02.png'),
     #                     save_data=False,
     #                     quad_correction=False,
-    #                     log_y=True
+    #                     log_y=False
     #                     )
     # IC_peak_energy_kappaR_combined_plot(config_data, a_bar,
     #                                     save_as=Path(
     #                                         '/home/andraz/ChargedShells/Figures/final_figures/peak_combined_kappaR.png')
     #                                     )
 
-    # a_bar = [0.1, 0.2, 0.3]
+    a_bar = [0.1, 0.2, 0.3]
     # IC_peak_energy_charge_plot(config_data, a_bar=a_bar, which='sep',
     #                            # save_as=Path('/home/andraz/ChargedShells/Figures/Emanuele_data/peak_sep_charge.png'),
     #                            )
-    # IC_peak_energy_charge_combined_plot(config_data, a_bar,
-    #                                     save_as=Path('/home/andraz/ChargedShells/Figures/final_figures/peak_combined_charge.png')
-    #                                     )
+    IC_peak_energy_charge_combined_plot(config_data, a_bar,
+                                        save_as=Path('/home/andraz/ChargedShells/Figures/final_figures/peak_combined_charge.png')
+                                        )
 
     kappaR = [0.01, 3.02407, 30]
-    IC_peak_energy_abar_plot(config_data, kappaR=kappaR, which='sep', min_abar=0.2, max_abar=0.8,
-                        save_as=Path('/home/andraz/ChargedShells/Figures/Emanuele_data/peak_sep_abar.png'),
-                        save_data=False,
-                        quad_correction=False,
-                        log_y=False
-                        )
+    # IC_peak_energy_abar_plot(config_data, kappaR=kappaR, which='sep', min_abar=0.2, max_abar=0.8,
+    #                     save_as=Path('/home/andraz/ChargedShells/Figures/Emanuele_data/peak_sep_abar.png'),
+    #                     save_data=False,
+    #                     quad_correction=False,
+    #                     log_y=False
+    #                     )
 
 
 if __name__ == '__main__':

analysis/peak_heigth_janus.py

@@ -0,0 +1,183 @@
+from matplotlib import gridspec
+from charged_shells import expansion, interactions, mapping, functions
+from charged_shells.parameters import ModelParams
+import numpy as np
+from typing import Literal
+from pathlib import Path
+from functools import partial
+import json
+from plot_settings import *
+from dataclasses import dataclass
+import peak_heigth
+
+Array = np.ndarray
+Expansion = expansion.Expansion
+
+
+class JanusPeakPP(peak_heigth.Peak):
+
+    def __init__(self, ex: Expansion, log_y: bool = False, kappaR_axis_in_expansion: int = None):
+        self.emanuele_data_column = 3
+        self.y_label = r'$V_\mathrm{PP, 1}$'
+        self.ex1 = ex.clone()
+        self.ex2 = ex.clone()
+        self.ex1.rotate_euler(beta=np.pi, alpha=0, gamma=0)
+        self.ex2.rotate_euler(beta=np.pi, alpha=0, gamma=0)
+        self.log_y = log_y
+        self.kappaR_axis_in_expansion = kappaR_axis_in_expansion
+
+
+class JanusPeakPPinv(peak_heigth.Peak):
+
+    def __init__(self, ex: Expansion, log_y: bool = False, kappaR_axis_in_expansion: int = None):
+        self.emanuele_data_column = 4
+        self.y_label = r'$V_\mathrm{PP, 2}$'
+        self.ex1 = ex.clone()
+        self.ex2 = ex.clone().inverse_sign(exclude_00=True)
+        self.ex1.rotate_euler(beta=np.pi, alpha=0, gamma=0)
+        self.ex2.rotate_euler(beta=np.pi, alpha=0, gamma=0)
+        self.log_y = log_y
+        self.kappaR_axis_in_expansion = kappaR_axis_in_expansion
+
+
+class JanusPeakEP(peak_heigth.Peak):
+
+    def __init__(self, ex: Expansion, log_y: bool = False, kappaR_axis_in_expansion: int = None):
+        self.emanuele_data_column = 4
+        self.y_label = r'$V_{EP, 1}$'
+        self.ex1 = ex.clone()
+        self.ex2 = ex.clone()
+        self.ex1.rotate_euler(beta=np.pi/2, alpha=0, gamma=0)
+        self.ex2.rotate_euler(beta=np.pi/2, alpha=0, gamma=0)
+        self.log_y = log_y
+        self.kappaR_axis_in_expansion = kappaR_axis_in_expansion
+
+
+class JanusPeakEPinv(peak_heigth.Peak):
+
+    def __init__(self, ex: Expansion, log_y: bool = False, kappaR_axis_in_expansion: int = None):
+        self.emanuele_data_column = 4
+        self.y_label = r'$V_{EP, 2}$'
+        self.ex1 = ex.clone()
+        self.ex2 = ex.clone().inverse_sign(exclude_00=True)
+        self.ex1.rotate_euler(beta=np.pi/2, alpha=0, gamma=0)
+        self.ex2.rotate_euler(beta=np.pi/2, alpha=0, gamma=0)
+        self.log_y = log_y
+        self.kappaR_axis_in_expansion = kappaR_axis_in_expansion
+
+
+def get_charge_energy_dicts(peak: peak_heigth.Peak, params: ModelParams, emanuele_data: Array, sigma0: Array, abar_cs: list, sigma_tilde=0.001):
+    abar_ic, inverse, counts = np.unique(emanuele_data[:, 1], return_counts=True, return_inverse=True)
+    # abar_ic = [0.1, 0.2, 0.3]
+
+    energy_fn = mapping.parameter_map_two_expansions(partial(interactions.charged_shell_energy, dist=2),
+                                                     peak.kappaR_axis_in_expansion)
+    energy = energy_fn(peak.ex1, peak.ex2, params)
+
+    data_dict_ic = {}
+    data_dict_cs = {}
+    k = 0
+    for i in range(len(abar_ic)):
+        ab = np.around(abar_ic[i], 5)
+        if ab in np.around(abar_cs, 5):
+            idx, = np.nonzero(inverse == i)
+            charge = emanuele_data[idx, 0]
+            en = emanuele_data[idx, peak.emanuele_data_column]
+            sort = np.argsort(charge)
+            if peak.log_y:
+                data_dict_ic[ab] = np.stack((charge[sort], np.abs(en)[sort])).T
+                data_dict_cs[ab] = np.stack((sigma0 / sigma_tilde, np.abs(energy[k]))).T
+            else:
+                data_dict_ic[ab] = np.stack((charge[sort], en[sort])).T
+                data_dict_cs[ab] = np.stack((sigma0 / sigma_tilde, energy[k])).T
+            k += 1
+
+    return data_dict_ic, data_dict_cs
+
+
+def IC_peak_energy_charge_combined_plot(config_data: dict,
+                        a_bar: list,
+                        R: float = 150,
+                        save_as: Path = None,
+                        log_y: bool = False):
+
+    # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_11"))
+    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS")).joinpath("FIG_5_JANUS_NEW_CHARGE")
+    em_data = np.load(em_data_path.joinpath("pair_energy_PP_janus.npz"))
+    data = em_data['changezc']
+
+    params = ModelParams(R=R, kappaR=3)
+    sigma0 = np.linspace(-0.00099, 0.00099, 300)
+    sigma_tilde = 0.00099
+
+    ex = expansion.MappedExpansionDipole(np.sort(np.array(a_bar)),  # sorting necessary as it is also in energy_dicts()
+                                       kappaR=3, sigma0=sigma0, l_max=20, sigma_tilde=sigma_tilde)
+
+    peak_pp = JanusPeakPP(ex, log_y)
+    peak_pp_inv = JanusPeakPPinv(ex, log_y)
+    peaks = [peak_pp, peak_pp_inv]
+
+    # peak_ep = JanusPeakEP(ex, log_y)
+    # peak_ep_inv = JanusPeakEPinv(ex, log_y)
+    # peaks = [peak_ep, peak_ep_inv]
+
+    data_ic = []
+    data_cs = []
+    for peak in peaks:
+        dict_ic, dict_cs = get_charge_energy_dicts(peak, params, data, sigma0, a_bar, sigma_tilde)
+        data_ic.append(dict_ic)
+        data_cs.append(dict_cs)
+
+    # fig, axs = plt.subplots(3, 1, figsize=(3, 7.8))
+    fig = plt.figure(figsize=(4, 1.7))
+    gs = gridspec.GridSpec(1, 2, figure=fig)
+    gs.update(left=0.125, right=0.975, top=0.99, bottom=0.21, wspace=0.35)
+    axs = [fig.add_subplot(gs[0, 0]), fig.add_subplot(gs[0, 1])]
+    for ax, data_dict_ic, data_dict_cs, peak in zip(axs, data_ic, data_cs, peaks):
+        colors = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
+        for ab in a_bar:
+            key = np.around(ab, 5)
+            current_color = next(colors)
+            ax.plot(data_dict_ic[key][:, 0], data_dict_ic[key][:, 1], label=rf'$\bar a = {key:.1f}$', c=current_color, linewidth=1.5)
+            ax.plot(data_dict_cs[key][:, 0], data_dict_cs[key][:, 1], ls='--', c=current_color, linewidth=1.5,
+                    # label=rf'$\bar a = {key:.1f}$'
+                    )
+        ax.legend(fontsize=9, ncol=1, frameon=False, handlelength=0.7, loc='upper right',
+                  bbox_to_anchor=(0.75, 1.03),
+                  )
+        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
+        # if ax == axs[-1]:
+        ax.set_xlabel(r'$\eta$', fontsize=11)
+
+        ax.set_ylabel(peak.y_label, fontsize=11)
+        if peak.log_y:
+            ax.set_yscale('log')
+        # ax.set_xscale('log')
+        ax.yaxis.set_label_coords(-0.18, 0.5)
+        ax.xaxis.set_label_coords(0.5, -0.12)
+    axs[0].set_ylim(-23, 23)
+    # plt.subplots_adjust(left=0.3)
+    # plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=300)
+    plt.show()
+
+
+def main():
+    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
+        config_data = json.load(config_file)
+
+    a_bar = [0.2, 0.5, 0.8]
+    # IC_peak_energy_kappaR_combined_plot(config_data, a_bar,
+    #                                     save_as=Path(
+    #                                         '/home/andraz/ChargedShells/Figures/final_figures/peak_combined_kappaR.png')
+    #                                     )
+
+    # a_bar = [0.1, 0.2, 0.3]
+    # IC_peak_energy_charge_combined_plot(config_data, a_bar,
+    #                                     save_as=Path('/home/andraz/ChargedShells/Figures/final_figures/janus_peak_combined_charge.png')
+    #                                     )
+
+if __name__ == '__main__':
+    main()
+

analysis/quad_path_plot.py

@@ -1,4 +1,5 @@
 import numpy as np
+from matplotlib import gridspec
 from matplotlib.lines import Line2D
 
 from charged_shells.rotational_path import PairRotationalPath, PathEnergyPlot
@@ -19,7 +20,7 @@ QuadPath.set_default_x_axis(zero_to_pi_half)
 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=pi_half_to_pi, 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")
 
@@ -211,7 +212,7 @@ def IC_sigma0_dependence(config_data: dict, save_as=None):
 def combined_distance_dependence(config_data: dict, dist: Array = 2 * np.array([1., 1.15, 1.3, 1.45]),
                                  kappaR: float = 3,
                                  abar: float = 0.5,
-                                 sigma_tilde=0.001,
+                                 sigma_tilde=0.00099,
                                  save_as=None):
 
     # em_data_path = Path(config_data["emanuele_data"]).joinpath("FIG_12")
@@ -237,8 +238,11 @@ def combined_distance_dependence(config_data: dict, dist: Array = 2 * np.array([
     labels = [rf'$\rho/R ={d}$' for d in dist]
 
     # additional legend
-    line1 = Line2D([0], [0], color='black', linewidth=1, label='ICi')
-    line2 = Line2D([0], [0], color='black', linestyle='--', linewidth=1, label='CSp')
+    line1 = Line2D([0], [0], color='black', linewidth=1.2, label='ICi')
+    line2 = Line2D([0], [0], color='black', linestyle='--', linewidth=1.2, label='CSp')
+
+    # np.savetxt("/home/andraz/Downloads/calibrate_CSp.dat", np.stack((x_axis, plots[0])).T)
+    # np.savetxt("/home/andraz/Downloads/calibrate_ICi.dat", ic_data[0])
 
     fig, ax = plt.subplots()
     for i, (d, en, label, c) in enumerate(zip(ic_data, plots, labels, COLORS)):
@@ -247,7 +251,7 @@ def combined_distance_dependence(config_data: dict, dist: Array = 2 * np.array([
             ax.plot(x_axis, en, ls='--', c=c)
     QuadPath.plot_style(fig, ax)
     main_legend = ax.get_legend()
-    extra_legend = ax.legend(handles=[line1, line2], loc='upper left', fontsize=15, frameon=False)
+    extra_legend = ax.legend(handles=[line1, line2], loc='upper left', fontsize=11, frameon=False)
     ax.add_artist(main_legend)
     ax.add_artist(extra_legend)
     if save_as is not None:
@@ -390,10 +394,10 @@ def combined_sigma0_dependence(config_data: dict, kappaR=3., abar=0.5, sigma0=(0
 
 def combined_all(config_data: dict, save_as=None):
 
-    sigma_tilde = 0.001
+    sigma_tilde = 0.00099
     kappaR_list = [1, 3, 10]
     abar_list = [0.5, 0.4, 0.3]
-    sigma0_list = [0.0002, 0.00, -0.0002]
+    sigma0_list = [0.000198, 0.00, -0.000198]
     kappaR = 3
     abar = 0.5
 
@@ -448,29 +452,38 @@ def combined_all(config_data: dict, save_as=None):
     energy_sigma0 = path_plot.evaluate_path()
     x_axis_sigma0 = path_plot.rot_path.stack_x_axes()
 
-    labels_sigma0 = [rf'$\eta={s0/sigma_tilde}$' for s0 in sigma0_list]
+    labels_sigma0 = [rf'$\eta={s0/sigma_tilde:.1f}$' for s0 in sigma0_list]
+    
+    line1 = Line2D([0], [0], color='black', linewidth=1.2, label='ICi')
+    line2 = Line2D([0], [0], color='black', linestyle='--', linewidth=1.2, label='CSp')
 
-    fig, axs = plt.subplots(3, 1, figsize=(6, 7.8))
+    # fig, axs = plt.subplots(3, 1, figsize=(6, 7.8))
+    fig = plt.figure(figsize=(4, 5.4))
+    gs = gridspec.GridSpec(3, 1, figure=fig)
+    gs.update(left=0.12, right=0.975, top=0.96, bottom=0.04, hspace=0.3)
+    axs = [fig.add_subplot(gs[0, 0]), fig.add_subplot(gs[1, 0]), fig.add_subplot(gs[2, 0])]
     for d, en, label, c in zip(ic_data_kappaR, energy_kappaR.T, labels_kappaR, COLOR_LIST):
-        axs[0].set_title('Screening', fontsize=15)
+        axs[0].set_title('Screening', fontsize=11, y=0.98)
         axs[0].plot(d[:, 0], d[:, 1], label=label, c=c)
         axs[0].plot(x_axis_kappaR, en, ls='--', c=c)
+        extra_legend = axs[0].legend(handles=[line1, line2], loc='upper left', fontsize=11, frameon=False)
+        axs[0].add_artist(extra_legend)
     QuadPath.plot_style(fig, axs[0], size=None)
     for d, en, label, c in zip(ic_data_abar, energy_abar.T, labels_abar, COLOR_LIST):
-        axs[1].set_title('Asymmetry', fontsize=15)
+        axs[1].set_title('Eccentricity', fontsize=11, y=0.98)
         axs[1].plot(d[:, 0], d[:, 1], label=label, c=c)
         axs[1].plot(x_axis_abar, en, ls='--', c=c)
     QuadPath.plot_style(fig, axs[1], size=None)
     for d, en, label, c in zip(ic_data_sigma0, energy_sigma0.T, labels_sigma0, COLOR_LIST):
-        axs[2].set_title('Net charge', fontsize=15)
+        axs[2].set_title('Net charge', fontsize=11, y=0.98)
         axs[2].plot(d[:, 0], d[:, 1], label=label, c=c)
         axs[2].plot(x_axis_sigma0, en, ls='--', c=c)
     for ax in axs:
-        ax.yaxis.set_label_coords(-0.09, 0.5)
+        ax.yaxis.set_label_coords(-0.08, 0.5)
     # axs[-1].set_xlabel('rotational path', fontsize=15)
     QuadPath.plot_style(fig, axs[2], size=None)
     for ax in axs:
-        ax.get_legend().set_bbox_to_anchor((0.65, 1))
+        ax.get_legend().set_bbox_to_anchor((0.6, 1))
     if save_as is not None:
         plt.savefig(save_as, dpi=300)
     plt.show()
@@ -525,13 +538,13 @@ def main():
 
     # combined_rescaled_distance_dependence(config_data)
 
-    combined_distance_dependence(config_data,
-                                 save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('quad_dist_dep.png')
-                                 )
+    # combined_distance_dependence(config_data,
+    #                              save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('quad_dist_dep.png')
+    #                              )
 
-    # combined_all(config_data,
-    #              save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('quad_combined_dep.png')
-    #              )
+    combined_all(config_data,
+                 save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('quad_combined_dep.png')
+                 )
 
 
 if __name__ == '__main__':

charged_shells/mapping.py

@@ -1,4 +1,4 @@
-from typing import Callable, TypeVar
+from typing import Callable, TypeVar, Protocol
 import numpy as np
 from charged_shells.expansion import Expansion
 from charged_shells.parameters import ModelParams
@@ -50,7 +50,11 @@ def unravel_expansion_over_axis(ex: Expansion, axis: int | None, param_list_len:
 
 
 SingleExpansionFn = Callable[[Expansion, ModelParams], T]
-TwoExpansionsFn = Callable[[Expansion, Expansion, ModelParams], T]
+# TwoExpansionsFn = Callable[[Expansion, Expansion, ModelParams], T]
+
+class TwoExpansionsFn(Protocol):
+    def __call__(self, ex1: Expansion, ex2: Expansion, params: ModelParams, **kwargs) -> T:
+        ...
 
 
 def parameter_map_single_expansion(f: SingleExpansionFn,
@@ -76,14 +80,14 @@ def parameter_map_two_expansions(f: TwoExpansionsFn,
                                  match_expansion_axis_to_params: int = None,
                                  ) -> TwoExpansionsFn:
 
-    def mapped_f(ex1: Expansion, ex2: Expansion, params: ModelParams):
+    def mapped_f(ex1: Expansion, ex2: Expansion, params: ModelParams, **kwargs):
         params_list = unravel_params(params)
         expansion_list1 = unravel_expansion_over_axis(ex1, match_expansion_axis_to_params, len(params_list))
         expansion_list2 = unravel_expansion_over_axis(ex2, match_expansion_axis_to_params, len(params_list))
 
         results = []
         for exp1, exp2, prms in zip(expansion_list1, expansion_list2, params_list):
-            results.append(f(exp1, exp2, prms))
+            results.append(f(exp1, exp2, prms, **kwargs))
 
         if match_expansion_axis_to_params is not None:
             # return the params-matched axis to where it belongs

charged_shells/rotational_path.py

@@ -68,18 +68,23 @@ class PairRotationalPath:
     def plot_style(self, fig, ax,
                    energy_units: interactions.EnergyUnit = 'kT',
                    legend: bool = True,
-                   size: tuple | None = (6, 2.5),
+                   size: tuple | None = (4, 1.7),
                    legend_loc: str = 'upper left'):
         ax.axhline(y=0, c='k', linestyle=':')
         if legend:
-            ax.legend(fontsize=15, frameon=False, loc=legend_loc, bbox_to_anchor=(0.6, 1))
-        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+            ax.legend(fontsize=10, frameon=False, loc=legend_loc, bbox_to_anchor=(0.57, 1))
+        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=11)
         # ax.set_xlabel('angle', fontsize=15)
-        ax.set_ylabel(f'$U [{energy_units}]$', fontsize=15)
-        ax.set_xticks(list(self.ticks.keys()), list(self.ticks.values()), fontsize=15)
+        if energy_units == 'kT':
+            energy_units = 'k_B T'
+        ax.set_ylabel(f'$V [{energy_units}]$', fontsize=11)
+        ax.set_xticks(list(self.ticks.keys()), list(self.ticks.values()), fontsize=11)
         if size is not None:
             fig.set_size_inches(size)
-        plt.tight_layout()
+        for line in ax.get_lines():
+            line.set_linewidth(1.2)
+        # plt.tight_layout()
+        plt.subplots_adjust(left=0.12, right=0.97, top=0.95, bottom=0.15)
 
 
 @dataclass

config.json

@@ -1,5 +1,5 @@
 {
-  "emanuele_data":  "/home/andraz/ChargedShells/data_Emanuele/2024-06-24",
+  "emanuele_data":  "/home/andraz/ChargedShells/data_Emanuele/2024-08-26",
   "figures": "/home/andraz/ChargedShells/Figures",
   "figure_data": "/home/andraz/ChargedShells/Figures/Data",
   "expansion_data": "/home/andraz/ChargedShells/ExpansionData/"