analysis/ICi_ICp_term_comparison.py

@@ -0,0 +1,134 @@
+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=FIGURES_PATH.joinpath('ici_icp_expansion_a08_kr10.png')
+    #          )
+
+    plot_exp2(1, [0.5, 0.8], legend=True,
+             # save_as=FIGURES_PATH.joinpath('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/config.py

@@ -0,0 +1,10 @@
+import json
+from pathlib import Path
+
+with open(Path(__file__).resolve().parent.joinpath('config.json')) as f:
+    config = json.load(f)
+
+ICI_DATA_PATH = Path(config["ICi_data_path"])
+FIGURES_PATH = Path(config["figures_path"])
+FIGURES_DATA_PATH = Path(config["figures_data_path"])
+EXPANSION_DATA_PATH = Path(config["expansion_data_path"])

analysis/dip_path_plot.py

@@ -1,10 +1,9 @@
-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 import interactions, charge_distributions
 from charged_shells.parameters import ModelParams
-from pathlib import Path
-import json
+from config import *
 from plot_settings import *
 
 Array = np.ndarray
@@ -43,16 +42,16 @@ 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):
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionDipole(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, DipolePath3, dist=2., params=params, match_expansion_axis_to_params=None)
@@ -63,7 +62,7 @@ def sections_plot(kappaR: float = 3, abar: float = 0.5, sigma_tilde=0.001, save_
 def kappaR_dependence(kappaR: Array, abar: float, sigma_tilde=0.001, save_as=None):
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionDipole(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, DipolePath3, dist=2., params=params, match_expansion_axis_to_params=0)
@@ -76,7 +75,7 @@ def kappaR_dependence(kappaR: Array, abar: float, sigma_tilde=0.001, save_as=Non
 def abar_dependence(abar: Array, kappaR: float, sigma_tilde=0.001, save_as=None):
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionDipole(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, DipolePath3, dist=2., params=params, match_expansion_axis_to_params=None)
@@ -89,7 +88,7 @@ def abar_dependence(abar: Array, kappaR: float, sigma_tilde=0.001, save_as=None)
 def sigma0_dependence(sigma0: Array, kappaR: float, abar: float, sigma_tilde=0.001, save_as=None):
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionDipole(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=sigma0)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=sigma0)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, DipolePath3, dist=2., params=params, match_expansion_axis_to_params=None)
@@ -99,13 +98,13 @@ def sigma0_dependence(sigma0: Array, kappaR: float, abar: float, sigma_tilde=0.0
                    save_as=save_as)
 
 
-def model_comparison(config_data: dict, save_as=None, save_data=False):
+def model_comparison(save_as=None, save_data=False):
     kappaR = 3
     params = ModelParams(R=150, kappaR=kappaR)
     a_bar = 0.5
     sigma_tilde = 0.001
 
-    ex1 = expansion.MappedExpansionDipole(a_bar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(a_bar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
     ex3 = ex1.clone().inverse_sign()
 
@@ -135,10 +134,10 @@ def model_comparison(config_data: dict, save_as=None, save_data=False):
     # pp_energy = interactions.charged_shell_energy(ex1new, ex2new, params)
     # print(f'PP energy: {pp_energy}')
 
-    # Emanuele data
-    em_data = np.load(Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_A").joinpath("pathway.npz"))['arr_0']
-    # 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")
+    # ICi data
+    em_data = np.load(ICI_DATA_PATH.joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_A").joinpath("pathway.npz"))['arr_0']
+    # em_data = np.load(ICI_DATA_PATH.joinpath("FIG_7").joinpath("pathway.npz"))['arr_0']
+    # em_data_path = (ICI_DATA_PATH.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']
     em_data, em_data_inv = em_data[:int(len(em_data) / 2)], em_data[int(len(em_data) / 2):]
@@ -164,9 +163,9 @@ def model_comparison(config_data: dict, save_as=None, save_data=False):
     plt.show()
 
 
-def combined_kappaR_dependence(config_data: dict, kappaR: list[int], abar: float, sigma_tilde=0.001, save_as=None):
+def combined_kappaR_dependence(kappaR: list[int], abar: float, sigma_tilde=0.001, save_as=None):
 
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_C")
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_C")
                     .joinpath("FIX_A").joinpath(f"ECC_{np.round(abar/2, 4)}"))
 
     ic_data = []
@@ -179,7 +178,7 @@ def combined_kappaR_dependence(config_data: dict, kappaR: list[int], abar: float
 
     params = ModelParams(R=150, kappaR=np.asarray(kappaR))
 
-    ex1 = expansion.MappedExpansionDipole(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
     ex3 = ex1.clone().inverse_sign()
 
@@ -204,9 +203,9 @@ def combined_kappaR_dependence(config_data: dict, kappaR: list[int], abar: float
     plt.show()
     
     
-def combined_abar_dependence(config_data: dict, kappaR: int, abar: list[float], sigma_tilde=0.001, save_as=None):
+def combined_abar_dependence(kappaR: int, abar: list[float], sigma_tilde=0.001, save_as=None):
     
-    em_data_path = Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_B").joinpath("FIX_M")
+    em_data_path = ICI_DATA_PATH.joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_B").joinpath("FIX_M")
 
     ic_data = []
     ic_data_inv = []
@@ -218,7 +217,7 @@ def combined_abar_dependence(config_data: dict, kappaR: int, abar: list[float],
 
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionDipole(np.asarray(abar), params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(np.asarray(abar), params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
     ex3 = ex1.clone().inverse_sign()
 
@@ -243,13 +242,13 @@ def combined_abar_dependence(config_data: dict, kappaR: int, abar: list[float],
     plt.show()
 
 
-def combined_sigma0_dependence(config_data: dict, kappaR=3., abar=0.5, sigma0=(-0.0002, 0.00, 0.0002), sigma_tilde=0.001, save_as=None):
+def combined_sigma0_dependence(kappaR=3., abar=0.5, sigma0=(-0.0002, 0.00, 0.0002), sigma_tilde=0.001, save_as=None):
 
-    em_data_path = Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_D").joinpath("CHANGE_ZC")
+    em_data_path = ICI_DATA_PATH.joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_D").joinpath("CHANGE_ZC")
     undercharged = np.load(em_data_path.joinpath("ZC_-56").joinpath("pathway.npz"))['arr_0']
     overcharged = np.load(em_data_path.joinpath("ZC_56").joinpath("pathway.npz"))['arr_0']
 
-    neutral_path = Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_B").joinpath("FIX_M")
+    neutral_path = ICI_DATA_PATH.joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_B").joinpath("FIX_M")
     neutral = np.load(neutral_path.joinpath(f"ECC_{np.round(abar/2, 4)}").joinpath(f"EMME_{int(kappaR)}.").joinpath("pathway.npz"))['arr_0']
 
     undercharged, undercharged_inv = undercharged[:int(len(undercharged) / 2)], undercharged[int(len(undercharged) / 2):]
@@ -260,7 +259,7 @@ def combined_sigma0_dependence(config_data: dict, kappaR=3., abar=0.5, sigma0=(-
     ic_data_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))
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=np.asarray(sigma0))
     ex2 = ex1.clone()
     ex3 = ex1.clone().inverse_sign(exclude_00=True)
 
@@ -285,16 +284,16 @@ def combined_sigma0_dependence(config_data: dict, kappaR=3., abar=0.5, sigma0=(-
     plt.show()
 
 
-def combined_all(config_data: dict, save_as=None):
+def combined_all(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
 
-    em_data_kappaR = (Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_C")
+    em_data_kappaR = (ICI_DATA_PATH.joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_C")
                       .joinpath("FIX_A").joinpath(f"ECC_{np.round(abar / 2, 4)}"))
     
     ic_data_kappaR = []
@@ -307,9 +306,9 @@ def combined_all(config_data: dict, save_as=None):
 
     params = ModelParams(R=150, kappaR=np.asarray(kappaR_list))
 
-    ex1 = expansion.MappedExpansionDipole(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(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()
@@ -318,7 +317,7 @@ def combined_all(config_data: dict, save_as=None):
     x_axis_kappaR = path_plot.rot_path.stack_x_axes()
     labels_kappaR = [rf'$\kappa R={kR}$' for kR in [1, 3, 10]]
 
-    em_data_abar = Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_B").joinpath("FIX_M")
+    em_data_abar = ICI_DATA_PATH.joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_B").joinpath("FIX_M")
 
     ic_data_abar = []
     ic_data_abar_inv = []
@@ -332,9 +331,9 @@ def combined_all(config_data: dict, save_as=None):
 
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionDipole(np.asarray(abar_list), params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(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()
@@ -343,10 +342,10 @@ def combined_all(config_data: dict, save_as=None):
     x_axis_abar = path_plot.rot_path.stack_x_axes()
     labels_abar = [rf'$\bar a={a}$' for a in abar_list]
 
-    em_data_charge = Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_D").joinpath("CHANGE_ZC")
+    em_data_charge = ICI_DATA_PATH.joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_D").joinpath("CHANGE_ZC")
     undercharged = np.load(em_data_charge.joinpath("ZC_-56").joinpath("pathway.npz"))['arr_0']
     overcharged = np.load(em_data_charge.joinpath("ZC_56").joinpath("pathway.npz"))['arr_0']
-    neutral_path = Path(config_data["emanuele_data"]).joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_B").joinpath("FIX_M")
+    neutral_path = ICI_DATA_PATH.joinpath("FIG_5_JANUS").joinpath("FIG_5_JANUS_B").joinpath("FIX_M")
     neutral = np.load(
         neutral_path.joinpath(f"ECC_{np.round(abar / 2, 4)}").joinpath(f"EMME_{int(kappaR)}.").joinpath("pathway.npz"))[
         'arr_0']
@@ -355,50 +354,57 @@ 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))
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(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()
+    # plt.tight_layout()
     if save_as is not None:
         plt.savefig(save_as, dpi=300)
     plt.show()
@@ -406,9 +412,6 @@ def combined_all(config_data: dict, save_as=None):
 
 if __name__ == '__main__':
 
-    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-        config_data = json.load(config_file)
-
     # sections_plot(save_as=Path("/home/andraz/ChargedShells/Figures/dipole_test_path.png"))
 
     # kappaR_dependence(np.array([3, 5]), 0.5,
@@ -423,22 +426,22 @@ if __name__ == '__main__':
     #                   save_as=Path("/home/andraz/ChargedShells/Figures/quadrupole_charge_dep_abar05_kappaR3.png")
     #                   )
 
-    # model_comparison(config_data,
-    #                  # save_as=Path(config_data["figures"]).joinpath("Emanuele_data").joinpath('IC_CS_janus_path.pdf')
+    # model_comparison(
+    #                  # save_as=FIGURES_PATH.joinpath("ICi_data").joinpath('IC_CS_janus_path.pdf')
     #                  )
 
-    # combined_kappaR_dependence(config_data, kappaR=[1, 3, 10], abar=0.5,
-    #                      # save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('janus_kappaR_dep.png')
+    # combined_kappaR_dependence(kappaR=[1, 3, 10], abar=0.5,
+    #                      # save_as=FIGURES_PATH.joinpath("final_figures").joinpath('janus_kappaR_dep.png')
     #                      )
     #
-    # combined_abar_dependence(config_data, kappaR=3, abar=[0.3, 0.4, 0.5],
-    #                          # save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('janus_abar_dep.png')
+    # combined_abar_dependence(kappaR=3, abar=[0.3, 0.4, 0.5],
+    #                          # save_as=FIGURES_PATH.joinpath("final_figures").joinpath('janus_abar_dep.png')
     #                          )
     #
-    # combined_sigma0_dependence(config_data,
-    #                      # save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('janus_charge_dep.png')
-    #                      )
+    # combined_sigma0_dependence(
+    #                      # save_as=FIGURES_PATH.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')
+    combined_all(
+                 # save_as=FIGURES_PATH.joinpath("final_figures").joinpath('janus_combined_dep.png')
                  )

analysis/energy_gap.py

@@ -1,14 +1,13 @@
 import matplotlib.pyplot as plt
 import numpy as np
-from charged_shells import expansion, interactions, mapping
+from charged_shells import expansion, interactions, mapping, charge_distributions
 from charged_shells.parameters import ModelParams
 from functools import partial
-from pathlib import Path
 from matplotlib import cm
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 from matplotlib.colors import TwoSlopeNorm
 from matplotlib.ticker import FuncFormatter
-import json
+from config import *
 
 Expansion = expansion.Expansion
 
@@ -32,7 +31,7 @@ def abar_kappaR_dependence(save_as=None):
     kappaR = np.linspace(0.01, 25, 25)
     a_bar = np.array([0.2, 0.4, 0.6, 0.8])
     params = ModelParams(R=150, kappaR=kappaR)
-    ex = expansion.MappedExpansionQuad(a_bar[:, None], kappaR[None, :], 0.001)
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar[:, None], kappaR[None, :], 0.001)
     # ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
     gap = energy_gap(ex, params, match_expansion_axis_to_params=1)
 
@@ -54,7 +53,7 @@ def abar_kappaR_dependence2(save_as=None):
     kappaR = np.array([1, 3, 10, 30])
     a_bar = np.linspace(0.2, 0.8, 30)
     params = ModelParams(R=150, kappaR=kappaR)
-    ex = expansion.MappedExpansionQuad(a_bar[:, None], kappaR[None, :], 0.001)
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar[:, None], kappaR[None, :], 0.001)
     # ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
     gap = energy_gap(ex, params, match_expansion_axis_to_params=1)
 
@@ -77,7 +76,7 @@ def charge_kappaR_dependence(a_bar, min_charge, max_charge, save_as=None, cmap=c
     sigma_tilde = 0.001
     params = ModelParams(R=150, kappaR=kappaR)
     charge = np.linspace(min_charge, max_charge, 100)
-    ex = expansion.MappedExpansionQuad(a_bar, kappaR, sigma_tilde, sigma0=charge)
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar, kappaR, sigma_tilde, sigma0=charge)
     # ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
     gap = energy_gap(ex, params, match_expansion_axis_to_params=0)
 
@@ -112,7 +111,7 @@ def charge_kappaR_dependence_heatmap(a_bar, min_charge, max_charge, save_as=None
     kappaR = np.linspace(0.01, 10, 50)
     params = ModelParams(R=150, kappaR=kappaR)
     charge = np.linspace(min_charge, max_charge, 100)
-    ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001, sigma0=charge)
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar, kappaR, 0.001, sigma0=charge)
     # ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
     gap = energy_gap(ex, params, match_expansion_axis_to_params=0)
 
@@ -152,8 +151,8 @@ def charge_kappaR_dependence_heatmap(a_bar, min_charge, max_charge, save_as=None
     plt.show()
 
 
-def IC_gap_plot(config_data: dict, save_as=None):
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
+def IC_gap_plot(save_as=None):
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_10"))
     em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
     for k in list(em_data.keys()):
         data = em_data[k]
@@ -163,8 +162,8 @@ def IC_gap_plot(config_data: dict, save_as=None):
         print('\n')
 
 
-def IC_gap_kappaR(config_data: dict, save_as=None):
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
+def IC_gap_kappaR(save_as=None):
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_10"))
     em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
     data = em_data['fixA']
 
@@ -180,8 +179,8 @@ def IC_gap_kappaR(config_data: dict, save_as=None):
     plt.show()
 
 
-def IC_gap_abar(config_data: dict, save_as=None):
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
+def IC_gap_abar(save_as=None):
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_10"))
     em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
     data = em_data['fixM']
 
@@ -197,9 +196,9 @@ def IC_gap_abar(config_data: dict, save_as=None):
     plt.show()
 
 
-def IC_gap_charge_at_abar(a_bar, config_data: dict, save_as=None, cmap=cm.coolwarm, which_change='changezp',
+def IC_gap_charge_at_abar(a_bar, save_as=None, cmap=cm.coolwarm, which_change='changezp',
                           eta_min: float = None, eta_max: float = None):
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_10"))
     em_data = np.load(em_data_path.joinpath("relative_gap_ZC.npz"))
     data = em_data[which_change]
 
@@ -264,14 +263,12 @@ def IC_gap_charge_at_abar(a_bar, config_data: dict, save_as=None, cmap=cm.coolwa
 
 def test_gap(a_bar, kappaR, charge):
     params = ModelParams(R=150, kappaR=kappaR)
-    ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001, sigma0=charge)
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar, kappaR, 0.001, sigma0=charge)
     gap = energy_gap(ex, params, match_expansion_axis_to_params=None)
     print(gap)
 
 
 def main():
-    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-        config_data = json.load(config_file)
 
     # test_gap(0.3, 10, charge=-0.003)
 
@@ -286,13 +283,13 @@ def main():
     #                          save_as=Path("/home/andraz/ChargedShells/Figures/full_amplitude_heatmap_abar05.png"),
     #                          cmap=cm.bwr)
 
-    # IC_gap_plot(config_data)
+    # IC_gap_plot()
 
-    # IC_gap_kappaR(config_data)
-    # IC_gap_abar(config_data)
+    # IC_gap_kappaR()
+    # IC_gap_abar()
 
-    IC_gap_charge_at_abar(0.3, config_data, which_change='changezc', eta_min=-2, eta_max=2,
-                          save_as=Path("/home/andraz/ChargedShells/Figures/Emanuele_data/IC_full_amplitude_charge_abar03.png")
+    IC_gap_charge_at_abar(0.3, which_change='changezc', eta_min=-2, eta_max=2,
+                          save_as=FIGURES_PATH.joinpath('ICi_data').joinpath('IC_full_amplitude_charge_abar03.png')
                           )
 
 

analysis/expanison_coef_export.py

@@ -1,22 +1,19 @@
-from charged_shells import expansion, parameters
-from pathlib import Path
-import json
+from charged_shells import parameters, charge_distributions
+from config import *
 import numpy as np
 
 
 def main():
-    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-        config_data = json.load(config_file)
 
     abar = np.array([0.2, 0.3, 0.4, 0.5, 0.5, 0.7, 0.8])
     params = parameters.ModelParams(kappaR=3, R=150)
-    ex = expansion.MappedExpansionDipole(abar, params.kappaR, 0.001, 30)
+    ex = charge_distributions.create_mapped_dipolar_expansion(abar, params.kappaR, 0.001, 30)
     print(ex.coefs.shape)
     l_arr, m_arr = ex.lm_arrays
     lm = np.stack((l_arr, m_arr)).T
     print(lm.shape)
 
-    np.savez(Path(config_data["expansion_data"]).joinpath("janus.npz"), abar=abar, lm=lm, coefs=ex.coefs)
+    np.savez(ICI_DATA_PATH.joinpath("janus.npz"), abar=abar, lm=lm, coefs=ex.coefs)
 
 
 

analysis/expansion_plot.py

@@ -1,9 +1,8 @@
 from charged_shells import expansion, parameters
 import numpy as np
 import matplotlib.pyplot as plt
-from pathlib import Path
 import plotly.graph_objects as go
-import json
+from config import *
 import quadrupole_model_mappings
 
 Expansion = expansion.Expansion
@@ -28,7 +27,7 @@ def plot_theta_profile_multiple(ex_list: list[Expansion], label_list, phi: float
     ax.set_ylabel(r'$\sigma$', fontsize=13)
     plt.legend(fontsize=12)
     plt.tight_layout()
-    plt.savefig(Path("/home/andraz/ChargedShells/Figures/charge_shape_comparison.png"), dpi=600)
+    plt.savefig(FIGURES_PATH.joinpath("charge_shape_comparison.png"), dpi=600)
     plt.show()
 
 
@@ -65,8 +64,6 @@ def plot_charge_3d(ex: Expansion, num_theta=100, num_phi=100, save_as: Path = No
 
 
 def main():
-    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-        config_data = json.load(config_file)
 
     params = parameters.ModelParams(kappaR=3, R=150)
     # ex = expansion.MappedExpansionQuad(0.328, params.kappaR, 0.001, 30)
@@ -77,7 +74,7 @@ def main():
     ex = quadrupole_model_mappings.ic_to_cap(0.001, 0.328, params, l_max=50)
     # print(np.real(ex.coefs))
     # plot_theta_profile(ex, num=1000, theta_end=2 * np.pi, phi=0)
-    plot_charge_3d(ex, save_as=Path(config_data["figures"]).joinpath("model_3D_cap.png"))
+    plot_charge_3d(ex, save_as=FIGURES_PATH.joinpath("model_3D_cap.png"))
 
     # new_coeffs = expanison.expansion_rotation(Quaternion(np.arange(20).reshape(5, 4)).normalized, ex.coeffs, ex.l_array)
     # print(new_coeffs.shape)

analysis/higher_multipole_matching.py

@@ -0,0 +1,129 @@
+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, abar: float,
+                            kappaR: Array, dist_cs: list, max_kappaR: float = 50, min_kappa: float = 0.01,
+                                      only_loaded: bool = False):
+
+    correct_abar_indices, = np.nonzero(np.isclose(emanuele_data[:, 1], abar))
+    emanuele_data = emanuele_data[correct_abar_indices]
+
+    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)
+            relevant_kR_indices = (emanuele_data[idx, 2] <= max_kappaR) * (min_kappa <= emanuele_data[idx, 2])
+            kR = emanuele_data[idx, 2][relevant_kR_indices]
+            en = emanuele_data[idx, peak.ici_data_column][relevant_kR_indices]
+            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(dist: list,
+                        which: RotConfig,
+                        abar=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 = (ICI_DATA_PATH.joinpath("FIG_11"))
+    em_data_path = (ICI_DATA_PATH.joinpath("TEST_HIGHER_MULTIPOLES"))
+    em_data = np.load(em_data_path.joinpath("higher_multiplot_energy_2.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 = charge_distributions.create_mapped_quad_expansion(abar, 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,abar, kappaR, dist, max_kappaR, min_kappaR)
+    data_dict_ic4, _ = get_kappaR_higher_multipole_dicts(peak, params, data4, abar, kappaR, dist, max_kappaR, min_kappaR)
+    data_dict_ic6, _ = get_kappaR_higher_multipole_dicts(peak, params, data6, abar, kappaR, dist, max_kappaR, min_kappaR)
+
+    colors = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
+
+    fig = plt.figure(figsize=(2, 2))
+    gs = gridspec.GridSpec(1, 1, figure=fig)
+    gs.update(left=0.25, right=0.95, top=0.9, 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.5, 0.5),
+              # bbox_to_anchor=(0.42, 0.9)
+              )
+    ax.set_title(f'$\\bar a = {abar:.1f}$, $\\rho = {key:.1f}$', loc='center',
+                 # x=0.95, y=0.8,
+                 y=0.96,
+                 # fontsize=9
+                 )
+    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__':
+
+    dist = [3.]
+    IC_peak_energy_plot(abar=0.5, dist=dist, which='pp', min_kappaR=0.1, max_kappaR=50,
+                        save_as=FIGURES_PATH.joinpath("ICi_data").joinpath("peak_pp_kappaR_higher_correction_abar05_rho30.png"),
+                        save_data=False,
+                        quad_correction=False,
+                        log_y=False
+                        )

analysis/patch_shape_comparison.py

@@ -1,9 +1,8 @@
-from charged_shells import expansion, potentials, patch_size
+from charged_shells import potentials, patch_size, charge_distributions
 import numpy as np
 from charged_shells.parameters import ModelParams
-from pathlib import Path
 from matplotlib.lines import Line2D
-import json
+from config import *
 import quadrupole_model_mappings
 from plot_settings import *
 
@@ -19,8 +18,8 @@ def ic_cs_comparison():
     phi = 0.
     dist = 1
 
-    ex_cp = expansion.MappedExpansionQuad(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30,
-                                          sigma0=sigma0)
+    ex_cp = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30,
+                                                                    sigma0=sigma0)
     potential_cp = potentials.charged_shell_potential(theta, phi, dist, ex_cp, params)
 
     potential_ic = potentials.inverse_patchy_particle_potential(theta, dist, a_bar, -2 * sigma_m + sigma0,
@@ -52,10 +51,10 @@ def models_comparison(save_data=False):
     sigma0 = 0  # taking this total charge parameter nonzero causes cap model to fail in finding the appropriate theta0
 
     def fn(x):
-        return expansion.MappedExpansionQuad(a_bar=x, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30, sigma0=sigma0)
+        return charge_distributions.create_mapped_quad_expansion(a_bar=x, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30, sigma0=sigma0)
 
     a_bar = patch_size.inverse_potential_patch_size(target_patch_size, fn, 0.5, params)
-    ex_point = expansion.MappedExpansionQuad(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30, sigma0=sigma0)
+    ex_point = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30, sigma0=sigma0)
 
     ex_gauss = quadrupole_model_mappings.ic_to_gauss(sigma_m, a_bar, params, l_max=30, sigma0=sigma0)
     ex_cap = quadrupole_model_mappings.ic_to_cap(sigma_m, a_bar, params, l_max=50, sigma0=sigma0)
@@ -73,7 +72,7 @@ def models_comparison(save_data=False):
     # print(potential)
 
     if save_data:
-        with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
+        with open(Path("/analysis/config.json")) as config_file:
             config_data = json.load(config_file)
         np.savez(Path(config_data["figure_data"]).joinpath("fig_6_shape_models.npz"),
                  ICi=np.stack((theta, potential_ic)).T,
@@ -114,12 +113,12 @@ def abar_comparison():
     dist = 1
 
     def fn1(x):
-        return expansion.MappedExpansionQuad(a_bar=x, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30)
+        return charge_distributions.create_mapped_quad_expansion(a_bar=x, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30)
 
     pots = []
     for ps in target_patch_sizes:
         a_bar = patch_size.inverse_potential_patch_size(ps, fn1, 0.5, params)
-        ex_point = expansion.MappedExpansionQuad(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30)
+        ex_point = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30)
         pots.append(potentials.charged_shell_potential(theta, phi, dist, ex_point, params))
 
     fig, ax = plt.subplots(figsize=(4.125, 3))
@@ -149,8 +148,8 @@ def charge_comparsion():
     phi = 0.
     dist = 1
 
-    ex_point = expansion.MappedExpansionQuad(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m,
-                                             l_max=30, sigma0=charges)
+    ex_point = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m,
+                                                                       l_max=30, sigma0=charges)
     pot = potentials.charged_shell_potential(theta, phi, dist, ex_point, params)
 
     fig, ax = plt.subplots(figsize=(4.125, 3))
@@ -165,7 +164,7 @@ def charge_comparsion():
     plt.xticks(custom_ticks, custom_labels, fontsize=15)
     plt.legend(fontsize=12)
     plt.tight_layout()
-    # plt.savefig(Path("/home/andraz/ChargedShells/Figures/potential_charge_comparison.pdf"), dpi=300)
+    # plt.savefig(FIGURES_PATH.joinpath('potential_charge_comparison.pdf'), dpi=300)
     plt.show()
 
 
@@ -184,32 +183,33 @@ def ic_cs_comparison2(save_data=False):
     dist = 1
 
     def fn1(x):
-        return expansion.MappedExpansionQuad(a_bar=x, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30, sigma0=0)
+        return charge_distributions.create_mapped_quad_expansion(a_bar=x, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30, sigma0=0)
     # a_bar is determined only for a neutral particle (patch size only well-defined in this case)
     # a_bar = patch_size.inverse_potential_patch_size(target_patch_size, fn1, 0.5, params)
     # print(a_bar)
 
     a_bar = 0.5
-    target_patch_size = patch_size.potential_patch_size(expansion.MappedExpansionQuad(a_bar=a_bar,
-                                                                                      kappaR=params.kappaR,
-                                                                                      sigma_tilde=sigma_m,
-                                                                                      l_max=30,
-                                                                                      sigma0=0),
-                                                        params)
+    target_patch_size = patch_size.potential_patch_size(
+        charge_distributions.create_mapped_quad_expansion(a_bar=a_bar,
+                                                                kappaR=params.kappaR,
+                                                                sigma_tilde=sigma_m,
+                                                                l_max=30,
+                                                                sigma0=0),
+        params)
 
     potential_ic = []
     potential_cs = []
     for sigma0 in sigma0_array:
 
-        ex_point = expansion.MappedExpansionQuad(a_bar=a_bar, kappaR=params.kappaR,
-                                                 sigma_tilde=sigma_m, l_max=30, sigma0=sigma0)
+        ex_point = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar, kappaR=params.kappaR,
+                                                                           sigma_tilde=sigma_m, l_max=30, sigma0=sigma0)
 
         potential_ic.append(potentials.inverse_patchy_particle_potential(theta, dist, a_bar, -2 * sigma_m + sigma0,
                                                                     (sigma_m, sigma_m), params, 30))
         potential_cs.append(potentials.charged_shell_potential(theta, phi, dist, ex_point, params))
 
     if save_data:
-        with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
+        with open(Path("/analysis/config.json")) as config_file:
             config_data = json.load(config_file)
         np.savez(Path(config_data["figure_data"]).joinpath("fig_6a.npz"),
                  ICi_eta_neg=np.stack((theta, potential_ic[0])).T,
@@ -219,34 +219,36 @@ 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'$\psi [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.savefig(Path("/home/andraz/ChargedShells/Figures/final_figures/potential_ic_cs_comparison.png"), dpi=300)
+    # 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(FIGURES_PATH.joinpath('final_figures').joinpath('potential_ic_cs_comparison.png'), dpi=300)
     plt.show()
 
 

analysis/patch_size_dependence.py

@@ -1,9 +1,8 @@
-from pathlib import Path
 import numpy as np
-from charged_shells import expansion, patch_size
+from charged_shells import patch_size, charge_distributions
 from charged_shells.parameters import ModelParams
 import time
-import json
+from config import *
 from plot_settings import *
 
 
@@ -11,11 +10,11 @@ def plot_abar_dependence(*, save=False, save_data=False):
     a_bar = np.linspace(0.2, 0.8, 101)
     kappaR = np.array([1, 3, 10])
     params = ModelParams(R=150, kappaR=kappaR)
-    ex = expansion.MappedExpansionQuad(a_bar=a_bar[:, None], sigma_tilde=0.001, l_max=30, kappaR=kappaR[None, :])
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar[:, None], sigma_tilde=0.001, l_max=30, kappaR=kappaR[None, :])
 
     ps = patch_size.potential_patch_size(ex, params, match_expansion_axis_to_params=1)
 
-    marked_point_ex = expansion.MappedExpansionQuad(a_bar=0.5, sigma_tilde=0.001, l_max=30, kappaR=3)
+    marked_point_ex = charge_distributions.create_mapped_quad_expansion(a_bar=0.5, sigma_tilde=0.001, l_max=30, kappaR=3)
     mark_ps = patch_size.potential_patch_size(marked_point_ex, ModelParams(R=150, kappaR=3),
                                               match_expansion_axis_to_params=None)
 
@@ -25,32 +24,32 @@ def plot_abar_dependence(*, save=False, save_data=False):
         data_dict = {}
         for key, patch in zip(["kR=01", "kR=1", "kR=3", "kR=10", "kR=30"], ps.T):
             data_dict[key] = np.stack((a_bar, patch)).T
-        with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-            config_data = json.load(config_file)
-        np.savez(Path(config_data["figure_data"]).joinpath("fig_6.npz"), **data_dict)
+        np.savez(FIGURES_DATA_PATH.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.savefig(FIGURES_PATH.joinpath('final_figures').joinpath('patch_size_potential.png'), dpi=300)
     plt.show()
 
 
 def plot_abar_charge_dependence(*, save=False):
     a_bar = np.linspace(0.2, 0.8, 101)
     kappaR = np.array([0.1, 1, 3, 10, 30])
-    ex = expansion.MappedExpansionQuad(a_bar=a_bar[:, None], sigma_tilde=0.001, l_max=30, kappaR=kappaR[None, :])
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar[:, None], sigma_tilde=0.001, l_max=30, kappaR=kappaR[None, :])
 
     ps = patch_size.charge_patch_size(ex)
 
@@ -65,14 +64,14 @@ def plot_abar_charge_dependence(*, save=False):
     plt.legend(fontsize=14)
     plt.tight_layout()
     if save:
-        plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_charge.pdf"), dpi=300)
+        plt.savefig(FIGURES_PATH.joinpath('patch_size_charge.pdf'), dpi=300)
     plt.show()
 
 
 def plot_kappaR_charge_dependence(*, normalized=False, save=False):
     a_bar = np.array([0.2, 0.4, 0.6, 0.8])
     kappaR = np.linspace(0.01, 30, 100)
-    ex = expansion.MappedExpansionQuad(a_bar=a_bar[None, :], sigma_tilde=0.001, l_max=30, kappaR=kappaR[:, None])
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar[None, :], sigma_tilde=0.001, l_max=30, kappaR=kappaR[:, None])
 
     ps = patch_size.charge_patch_size(ex)
     if normalized:
@@ -89,7 +88,7 @@ def plot_kappaR_charge_dependence(*, normalized=False, save=False):
     plt.legend(fontsize=14)
     plt.tight_layout()
     if save:
-        plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_charge_nonmonotonicity.pdf"), dpi=300)
+        plt.savefig(FIGURES_PATH.joinpath('patch_size_charge_nonmonotonicity.pdf'), dpi=300)
     plt.show()
 
 
@@ -100,7 +99,7 @@ def plot_sigma0_dependence(*, save=False):
     # sigma0 = np.linspace(-0.001, 0.0015, 6)
     sigma0 = np.array([-0.0002, 0, 0.0002])
     params = ModelParams(R=150, kappaR=kappaR)
-    ex = expansion.MappedExpansionQuad(a_bar=0.3, sigma_tilde=0.001, l_max=30, kappaR=kappaR, sigma0=sigma0)
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar=0.3, sigma_tilde=0.001, l_max=30, kappaR=kappaR, sigma0=sigma0)
 
     t0 = time.perf_counter()
     ps = patch_size.potential_patch_size(ex, params, match_expansion_axis_to_params=0, skip_nozero_cases=True)
@@ -119,7 +118,7 @@ def plot_sigma0_dependence(*, save=False):
     plt.legend(fontsize=14, loc='upper right')
     plt.tight_layout()
     if save:
-        plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential_charge_abar03.pdf"), dpi=300)
+        plt.savefig(FIGURES_PATH.joinpath('patch_size_potential_charge_abar03.pdf'), dpi=300)
     plt.show()
 
 
@@ -130,8 +129,8 @@ def plot_sigma0_dependence_relative(*, save=False):
     # sigma0 = np.linspace(-0.001, 0.0015, 6)
     sigma0 = np.array([-0.0002, 0, 0.0002])
     params = ModelParams(R=150, kappaR=kappaR)
-    ex = expansion.MappedExpansionQuad(a_bar=0.8, sigma_tilde=0.001, l_max=30, kappaR=kappaR, sigma0=sigma0)
-    ex_neutral = expansion.MappedExpansionQuad(a_bar=0.3, sigma_tilde=0.001, l_max=30, kappaR=kappaR, sigma0=0)
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar=0.8, sigma_tilde=0.001, l_max=30, kappaR=kappaR, sigma0=sigma0)
+    ex_neutral = charge_distributions.create_mapped_quad_expansion(a_bar=0.3, sigma_tilde=0.001, l_max=30, kappaR=kappaR, sigma0=0)
 
     ps = patch_size.potential_patch_size(ex, params, match_expansion_axis_to_params=0, skip_nozero_cases=True)
     ps_neutral = patch_size.potential_patch_size(ex_neutral, params, match_expansion_axis_to_params=0)
@@ -149,7 +148,7 @@ def plot_sigma0_dependence_relative(*, save=False):
     plt.legend(fontsize=14)
     plt.tight_layout()
     if save:
-        plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential_charge_abar03_relative.pdf"), dpi=300)
+        plt.savefig(FIGURES_PATH.joinpath('patch_size_potential_charge_abar03_relative.pdf'), dpi=300)
     plt.show()
 
 

analysis/peak_heigth.py

@@ -1,10 +1,10 @@
-from charged_shells import expansion, interactions, mapping, functions
+from matplotlib import gridspec
+from charged_shells import expansion, interactions, mapping, functions, charge_distributions
 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 config import *
 from plot_settings import *
 from dataclasses import dataclass
 
@@ -19,7 +19,7 @@ class Peak:
     y_label: str
     ex1: Expansion
     ex2: Expansion
-    emanuele_data_column: int
+    ici_data_column: int
     log_y: bool
     kappaR_axis_in_expansion: int = None
 
@@ -27,8 +27,8 @@ class Peak:
 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.ici_data_column = 4
+        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)
@@ -39,8 +39,8 @@ class PeakEP(Peak):
 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.ici_data_column = 3
+        self.y_label = r'$V_\mathrm{PP}$'
         self.ex1 = ex.clone()
         self.ex2 = ex.clone()
         self.log_y = log_y
@@ -50,12 +50,12 @@ class PeakPP(Peak):
 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.ici_data_column = 5
+        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
 
@@ -75,7 +75,7 @@ def get_charge_energy_dicts(peak: Peak, params: ModelParams, emanuele_data: Arra
         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]
+            en = emanuele_data[idx, peak.ici_data_column]
             sort = np.argsort(charge)
             if peak.log_y:
                 data_dict_ic[ab] = np.stack((charge[sort] / 280 + 2, np.abs(en)[sort])).T
@@ -88,12 +88,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):
-    abar_ic, inverse = np.unique(emanuele_data[:, 1], return_inverse=True)
+def get_kappaR_energy_dicts(peak: Peak, params: ModelParams, ici_data: Array,
+                            kappaR: Array, abar_cs: list, max_kappaR: float = 50, only_loaded: bool = False):
+    abar_ic, inverse = np.unique(ici_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 = {}
@@ -102,8 +106,8 @@ def get_kappaR_energy_dicts(peak: Peak, params: ModelParams, emanuele_data: Arra
         ab = np.around(abar_ic[i], 5)
         if ab in np.around(abar_cs, 5):
             idx, = np.nonzero(inverse == i)
-            kR = emanuele_data[idx, 2][emanuele_data[idx, 2] <= max_kappaR]
-            en = emanuele_data[idx, peak.emanuele_data_column][emanuele_data[idx, 1] <= max_kappaR]
+            kR = ici_data[idx, 2][ici_data[idx, 2] <= max_kappaR]
+            en = ici_data[idx, peak.ici_data_column][ici_data[idx, 1] <= max_kappaR]
             sort = np.argsort(kR)
             if peak.log_y:
                 data_dict_ic[ab] = np.stack((kR[sort], np.abs(en)[sort])).T
@@ -116,9 +120,9 @@ def get_kappaR_energy_dicts(peak: Peak, params: ModelParams, emanuele_data: Arra
     return data_dict_ic, data_dict_cs
 
 
-def get_abar_energy_dicts(peak: Peak, params: ModelParams, emanuele_data: Array, a_bar: Array, kR_cs: list, max_abar: float = 0.8):
+def get_abar_energy_dicts(peak: Peak, params: ModelParams, ici_data: Array, a_bar: Array, kR_cs: list, max_abar: float = 0.8):
 
-    kR_ic, inverse = np.unique(emanuele_data[:, 2], return_inverse=True)
+    kR_ic, inverse = np.unique(ici_data[:, 2], return_inverse=True)
 
     energy_fn = mapping.parameter_map_two_expansions(partial(interactions.charged_shell_energy, dist=2),
                                                      peak.kappaR_axis_in_expansion)
@@ -131,8 +135,8 @@ def get_abar_energy_dicts(peak: Peak, params: ModelParams, emanuele_data: Array,
         kr = np.around(kR_ic[i], 5)
         if kr in np.around(kR_cs, 5):
             idx, = np.nonzero(inverse == i)
-            abar = emanuele_data[idx, 1][emanuele_data[idx, 1] <= max_abar]
-            en = emanuele_data[idx, peak.emanuele_data_column][emanuele_data[idx, 1] <= max_abar]
+            abar = ici_data[idx, 1][ici_data[idx, 1] <= max_abar]
+            en = ici_data[idx, peak.ici_data_column][ici_data[idx, 1] <= max_abar]
             sort = np.argsort(abar)
             if peak.log_y:
                 data_dict_ic[kr] = np.stack((abar[sort], np.abs(en)[sort])).T
@@ -151,14 +155,11 @@ def save_data(data_dict_ic, data_dict_cs, a_bar: list):
         cs_save = {str(key): val for key, val in data_dict_cs.items()}
         ic_save['abar'] = a_bar
         cs_save['abar'] = a_bar
-        with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-            config_data = json.load(config_file)
-        np.savez(Path(config_data["figure_data"]).joinpath("fig_11_IC.npz"), **ic_save)
-        np.savez(Path(config_data["figure_data"]).joinpath("fig_11_CS.npz"), **cs_save)
+        np.savez(FIGURES_DATA_PATH.joinpath("fig_11_IC.npz"), **ic_save)
+        np.savez(FIGURES_DATA_PATH.joinpath("fig_11_CS.npz"), **cs_save)
 
 
-def IC_peak_energy_plot(config_data: dict,
-                        a_bar: list,
+def IC_peak_energy_plot(a_bar: list,
                         which: RotConfig,
                         min_kappaR: float = 0.01,
                         max_kappaR: float = 50.,
@@ -168,8 +169,8 @@ def IC_peak_energy_plot(config_data: dict,
                         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("FIG_4_PANELS_BDF"))
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_11"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_PANELS_BDF"))
     em_data = np.load(em_data_path.joinpath("newpair_energy.npz"))
     data = em_data['fixA']
 
@@ -177,7 +178,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 = charge_distributions.create_mapped_quad_expansion(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,31 +191,47 @@ 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()
 
 
-def IC_peak_energy_abar_plot(config_data: dict,
+def IC_peak_energy_abar_plot(
                         kappaR: list,
                         which: RotConfig,
                         min_abar: float = 0.01,
@@ -225,8 +242,8 @@ def IC_peak_energy_abar_plot(config_data: dict,
                         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("FIG_4_PANELS_BDF"))
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_11"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_PANELS_BDF"))
     em_data = np.load(em_data_path.joinpath("newpair_energy.npz"))
     data = em_data['fixA']
 
@@ -234,7 +251,7 @@ def IC_peak_energy_abar_plot(config_data: dict,
     a_bar = np.linspace(min_abar, max_abar, num)
     params = ModelParams(R=R, kappaR=np.array(kappaR))
 
-    ex = expansion.MappedExpansionQuad(np.array(a_bar)[None, :], params.kappaR[:, None], 0.001, l_max=20)
+    ex = charge_distributions.create_mapped_quad_expansion(np.array(a_bar)[None, :], params.kappaR[:, None], 0.001, l_max=20)
 
     if which == 'ep':
         peak = PeakEP(ex, log_y, kappaR_axis_in_expansion=0)
@@ -272,8 +289,7 @@ def IC_peak_energy_abar_plot(config_data: dict,
     plt.show()
 
 
-def IC_peak_energy_charge_plot(config_data: dict,
-                        a_bar: list,
+def IC_peak_energy_charge_plot(a_bar: list,
                         which: RotConfig,
                         max_kappaR: float = 30.,
                         R: float = 150,
@@ -282,8 +298,8 @@ def IC_peak_energy_charge_plot(config_data: dict,
                         quad_correction: bool = False,
                         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_4_PANELS_BDF"))
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_11"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_PANELS_BDF"))
     em_data = np.load(em_data_path.joinpath("newpair_energy.npz"))
     data = em_data['changeZc']
 
@@ -291,7 +307,7 @@ def IC_peak_energy_charge_plot(config_data: dict,
     sigma0 = np.linspace(-0.0003, 0.0003, 300)
     sigma_tilde = 0.001
 
-    ex = expansion.MappedExpansionQuad(np.array(a_bar), kappaR=3, sigma0=sigma0, l_max=20, sigma_tilde=sigma_tilde)
+    ex = charge_distributions.create_mapped_quad_expansion(np.array(a_bar), kappaR=3, sigma0=sigma0, l_max=20, sigma_tilde=sigma_tilde)
 
     if which == 'ep':
         peak = PeakEP(ex, log_y)
@@ -329,16 +345,15 @@ def IC_peak_energy_charge_plot(config_data: dict,
     plt.show()
 
 
-def IC_peak_energy_kappaR_combined_plot(config_data: dict,
-                        a_bar: list,
+def IC_peak_energy_kappaR_combined_plot(a_bar: list,
                         R: float = 150,
                         save_as: Path = None,
                         min_kappaR: float = 0.01,
                         max_kappaR: float = 50.,
                         ):
 
-    # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_11"))
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_4_PANELS_BDF"))
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_11"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_PANELS_BDF"))
     em_data = np.load(em_data_path.joinpath("newpair_energy.npz"))
     data = em_data['fixA']
 
@@ -346,8 +361,8 @@ def IC_peak_energy_kappaR_combined_plot(config_data: dict,
     kappaR = np.geomspace(min_kappaR, max_kappaR, num)
     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)
+    ex = charge_distributions.create_mapped_quad_expansion(np.sort(np.array(a_bar))[:, None],  # sorting necessary as it is also in energy_dicts()
+                                       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,57 +372,62 @@ 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()
 
 
-def IC_peak_energy_charge_combined_plot(config_data: dict,
-                        a_bar: list,
+def IC_peak_energy_charge_combined_plot(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_4_PANELS_BDF"))
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_11"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_PANELS_BDF"))
     em_data = np.load(em_data_path.joinpath("newpair_energy.npz"))
     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)
+    ex = charge_distributions.create_mapped_quad_expansion(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_ep = PeakEP(ex, log_y)
     peak_pp = PeakPP(ex, log_y)
@@ -421,64 +441,71 @@ 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(a_bar=a_bar, which='pp', min_kappaR=0.01, max_kappaR=50,
+    #                     save_as=FIGURES_PATH.joinpath('ICi_data').joinpath('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')
+    # IC_peak_energy_kappaR_combined_plot(a_bar,
+    #                                     save_as=FIGURES_PATH.joinpath('final_figures').joinpath('peak_combined_kappaR.png')
     #                                     )
 
-    # 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'),
+    a_bar = [0.1, 0.2, 0.3]
+    # IC_peak_energy_charge_plot(a_bar=a_bar, which='sep',
+    #                            # save_as=FIGURES_PATH.joinpath('ICi_data').joinpath('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(a_bar,
+                                        save_as=FIGURES_PATH.joinpath('final_figures').joinpath('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(kappaR=kappaR, which='sep', min_abar=0.2, max_abar=0.8,
+    #                     save_as=FIGURES_PATH.joinpath('ICi_data').joinpath('peak_sep_abar.png'),
+    #                     save_data=False,
+    #                     quad_correction=False,
+    #                     log_y=False
+    #                     )
 
 
 if __name__ == '__main__':

analysis/peak_heigth_janus.py

@@ -0,0 +1,173 @@
+from matplotlib import gridspec
+from charged_shells import expansion, interactions, mapping, functions, charge_distributions
+from charged_shells.parameters import ModelParams
+import numpy as np
+from typing import Literal
+from functools import partial
+from config import *
+from plot_settings import *
+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.ici_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.ici_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.ici_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.ici_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, ici_data: Array, sigma0: Array, abar_cs: list, sigma_tilde=0.001):
+    abar_ic, inverse, counts = np.unique(ici_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 = ici_data[idx, 0]
+            en = ici_data[idx, peak.ici_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(a_bar: list,
+                        R: float = 150,
+                        save_as: Path = None,
+                        log_y: bool = False):
+
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_11"))
+    em_data_path = (ICI_DATA_PATH.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 = charge_distributions.create_mapped_dipolar_expansion(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():
+
+    a_bar = [0.2, 0.5, 0.8]
+    # a_bar = [0.1, 0.2, 0.3]
+    IC_peak_energy_charge_combined_plot(a_bar,
+                                        save_as=FIGURES_PATH.joinpath('final_figures').joinpath('janus_peak_combined_charge.png')
+                                        )
+
+if __name__ == '__main__':
+    main()
+

analysis/plot_settings.py

@@ -1,8 +1,6 @@
 from matplotlib import cm
 import matplotlib.pyplot as plt
 from itertools import cycle
-import matplotlib.legend as mlegend
-from matplotlib.legend_handler import HandlerBase
 
 # all the common imports and settings for plotting
 

analysis/quad_path_plot.py

@@ -1,12 +1,10 @@
 import numpy as np
+from matplotlib import gridspec
 from matplotlib.lines import Line2D
-
 from charged_shells.rotational_path import PairRotationalPath, PathEnergyPlot
-from charged_shells import expansion
+from charged_shells import charge_distributions
 from charged_shells.parameters import ModelParams
-from pathlib import Path
-import json
-import quadrupole_model_mappings
+from config import *
 from plot_settings import *
 
 Array = np.ndarray
@@ -19,7 +17,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")
 
@@ -30,7 +28,7 @@ def model_comparison(config_data: dict, save_as=None, save_data=False):
     a_bar = 0.5
     sigma_tilde = 0.001
 
-    ex1 = expansion.MappedExpansionQuad(a_bar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(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
@@ -57,9 +55,9 @@ def model_comparison(config_data: dict, save_as=None, save_data=False):
     # print(f'PP energy: {pp_energy}')
 
     # Emanuele data
-    em_data = np.load(Path(config_data["emanuele_data"]).joinpath("FIG_3C").joinpath("pathway.npz"))['arr_0']
-    # 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")
+    em_data = np.load(ICI_DATA_PATH.joinpath("FIG_3C").joinpath("pathway.npz"))['arr_0']
+    # em_data = np.load(ICI_DATA_PATH.joinpath("FIG_7").joinpath("pathway.npz"))['arr_0']
+    # em_data_path = (ICI_DATA_PATH.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']
 
@@ -86,7 +84,7 @@ def model_comparison(config_data: dict, save_as=None, save_data=False):
 def kappaR_dependence(kappaR: Array, abar: float, sigma_tilde=0.001, save_as=None):
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=0)
@@ -99,7 +97,7 @@ def kappaR_dependence(kappaR: Array, abar: float, sigma_tilde=0.001, save_as=Non
 def abar_dependence(abar: Array, kappaR: float, sigma_tilde=0.001, save_as=None):
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=None)
@@ -112,7 +110,7 @@ def abar_dependence(abar: Array, kappaR: float, sigma_tilde=0.001, save_as=None)
 def sigma0_dependence(sigma0: Array, kappaR: float, abar: float, sigma_tilde=0.001, save_as=None):
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=sigma0)
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=sigma0)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=None)
@@ -125,7 +123,7 @@ def sigma0_dependence(sigma0: Array, kappaR: float, abar: float, sigma_tilde=0.0
 def distance_dependence(dist: Array, kappaR: float, abar: float, sigma_tilde=0.001, save_as=None):
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     plots = []
@@ -147,11 +145,11 @@ def distance_dependence(dist: Array, kappaR: float, abar: float, sigma_tilde=0.0
     plt.show()
 
 
-def IC_kappaR_dependence(config_data: dict, save_as=None):
+def IC_kappaR_dependence(save_as=None):
 
-    # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_8").joinpath("FIXEDCHARGE")
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_8").joinpath("FIXEDCHARGE")
     #                 .joinpath("FIX_A").joinpath("ECC_0.25"))
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE")
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE")
                     .joinpath("FIX_A").joinpath("ECC_0.25"))
     kR1 = np.load(em_data_path.joinpath("EMME_1.").joinpath("pathway.npz"))['arr_0']
     kR3 = np.load(em_data_path.joinpath("EMME_3.").joinpath("pathway.npz"))['arr_0']
@@ -169,10 +167,10 @@ def IC_kappaR_dependence(config_data: dict, save_as=None):
     plt.show()
 
 
-def IC_abar_dependence(config_data: dict, save_as=None):
+def IC_abar_dependence(save_as=None):
 
-    # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_8").joinpath("FIXEDCHARGE").joinpath("FIX_M"))
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE").joinpath("FIX_M"))
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_8").joinpath("FIXEDCHARGE").joinpath("FIX_M"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE").joinpath("FIX_M"))
     a03 = np.load(em_data_path.joinpath("ECC_0.15").joinpath("EMME_3.").joinpath("pathway.npz"))['arr_0']
     a04 = np.load(em_data_path.joinpath("ECC_0.2").joinpath("EMME_3.").joinpath("pathway.npz"))['arr_0']
     a05 = np.load(em_data_path.joinpath("ECC_0.25").joinpath("EMME_3.").joinpath("pathway.npz"))['arr_0']
@@ -189,9 +187,9 @@ def IC_abar_dependence(config_data: dict, save_as=None):
     plt.show()
 
 
-def IC_sigma0_dependence(config_data: dict, save_as=None):
-    # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_8").joinpath("CHARGE_ZC"))
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_4_Panels_ACE").joinpath("CHARGE_ZC"))
+def IC_sigma0_dependence(save_as=None):
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_8").joinpath("CHARGE_ZC"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_Panels_ACE").joinpath("CHARGE_ZC"))
     undercharged = np.load(em_data_path.joinpath("ZC_-277.27").joinpath("pathway.npz"))['arr_0']
     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']
@@ -208,24 +206,24 @@ def IC_sigma0_dependence(config_data: dict, save_as=None):
     plt.show()
 
 
-def combined_distance_dependence(config_data: dict, dist: Array = 2 * np.array([1., 1.15, 1.3, 1.45]),
+def combined_distance_dependence(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")
-    em_data_path = Path(config_data["emanuele_data"]).joinpath("FIG_3D_LONG_DIST")
+    # em_data_path = ICI_DATA_PATH.joinpath("FIG_12")
+    em_data_path = ICI_DATA_PATH.joinpath("FIG_3D_LONG_DIST")
     em_data = np.load(em_data_path.joinpath("pathway_fig12A.npz"))
 
-    em_data_d2 = np.load(Path(config_data["emanuele_data"]).joinpath("FIG_3C").joinpath("pathway.npz"))['arr_0']
+    em_data_d2 = np.load(ICI_DATA_PATH.joinpath("FIG_3C").joinpath("pathway.npz"))['arr_0']
 
     ic_data = [em_data_d2]
     for key, d in em_data.items():
         ic_data.append(d)
 
     params = ModelParams(R=150, kappaR=kappaR)
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     plots = []
@@ -237,8 +235,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 +248,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:
@@ -255,15 +256,14 @@ def combined_distance_dependence(config_data: dict, dist: Array = 2 * np.array([
     plt.show()
 
 
-def combined_rescaled_distance_dependence(config_data: dict,
-                                          dist: Array = 2 * np.array([1, 1.5, 2, 3, 5, 10]),
+def combined_rescaled_distance_dependence(dist: Array = 2 * np.array([1, 1.5, 2, 3, 5, 10]),
                                           kappaR: float = 3,
                                           abar: float = 0.5,
                                           sigma_tilde=0.001,
                                           save_as=None):
 
-    # em_data_path = Path(config_data["emanuele_data"]).joinpath("FIG_12")
-    em_data_path = Path(config_data["emanuele_data"]).joinpath("FIG_3D_LONG_DIST")
+    # em_data_path = ICI_DATA_PATH.joinpath("FIG_12")
+    em_data_path = ICI_DATA_PATH.joinpath("FIG_3D_LONG_DIST")
     em_data = np.load(em_data_path.joinpath("pathway_fig12B.npz"))
     ic_data = []
     for key, d in em_data.items():
@@ -273,7 +273,7 @@ def combined_rescaled_distance_dependence(config_data: dict,
         ic_data.append(d)
 
     params = ModelParams(R=150, kappaR=kappaR)
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     plots = []
@@ -297,11 +297,11 @@ def combined_rescaled_distance_dependence(config_data: dict,
     plt.show()
 
 
-def combined_kappaR_dependence(config_data: dict, kappaR: list[int], abar: float, sigma_tilde=0.001, save_as=None):
+def combined_kappaR_dependence(kappaR: list[int], abar: float, sigma_tilde=0.001, save_as=None):
 
-    # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_8").joinpath("FIXEDCHARGE")
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_8").joinpath("FIXEDCHARGE")
     #                 .joinpath("FIX_A").joinpath(f"ECC_{np.round(abar/2, 4)}"))
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE")
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE")
                     .joinpath("FIX_A").joinpath(f"ECC_{np.round(abar/2, 4)}"))
 
     ic_data = []
@@ -310,7 +310,7 @@ def combined_kappaR_dependence(config_data: dict, kappaR: list[int], abar: float
 
     params = ModelParams(R=150, kappaR=np.asarray(kappaR))
 
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=0)
@@ -329,10 +329,10 @@ def combined_kappaR_dependence(config_data: dict, kappaR: list[int], abar: float
     plt.show()
 
 
-def combined_abar_dependence(config_data: dict, kappaR: int, abar: list[float], sigma_tilde=0.001, save_as=None):
+def combined_abar_dependence(kappaR: int, abar: list[float], sigma_tilde=0.001, save_as=None):
 
-    # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_8").joinpath("FIXEDCHARGE").joinpath("FIX_M"))
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE").joinpath("FIX_M"))
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_8").joinpath("FIXEDCHARGE").joinpath("FIX_M"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE").joinpath("FIX_M"))
 
     ic_data = []
     for ab in abar:
@@ -341,7 +341,7 @@ def combined_abar_dependence(config_data: dict, kappaR: int, abar: list[float],
 
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionQuad(np.asarray(abar), params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(np.asarray(abar), params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=None)
@@ -360,16 +360,16 @@ def combined_abar_dependence(config_data: dict, kappaR: int, abar: list[float],
     plt.show()
 
 
-def combined_sigma0_dependence(config_data: dict, kappaR=3., abar=0.5, sigma0=(0.0002, 0.00, -0.0002), sigma_tilde=0.001, save_as=None):
-    # em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_8").joinpath("CHARGE_ZC"))
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_4_Panels_ACE").joinpath("CHARGE_ZC"))
+def combined_sigma0_dependence(kappaR=3., abar=0.5, sigma0=(0.0002, 0.00, -0.0002), sigma_tilde=0.001, save_as=None):
+    # em_data_path = (ICI_DATA_PATH.joinpath("FIG_8").joinpath("CHARGE_ZC"))
+    em_data_path = (ICI_DATA_PATH.joinpath("FIG_4_Panels_ACE").joinpath("CHARGE_ZC"))
     undercharged = np.load(em_data_path.joinpath("ZC_-503").joinpath("pathway.npz"))['arr_0']
     neutral = np.load(em_data_path.joinpath("ZC_-560").joinpath("pathway.npz"))['arr_0']
     overchargerd = np.load(em_data_path.joinpath("ZC_-617").joinpath("pathway.npz"))['arr_0']
     ic_data = [undercharged, neutral, overchargerd]
 
     params = ModelParams(R=150, kappaR=kappaR)
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=np.asarray(sigma0))
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=np.asarray(sigma0))
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=None)
@@ -388,17 +388,16 @@ def combined_sigma0_dependence(config_data: dict, kappaR=3., abar=0.5, sigma0=(0
     plt.show()
 
 
-def combined_all(config_data: dict, save_as=None):
+def combined_all(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
 
-
-    em_data_kappaR = (Path(config_data["emanuele_data"]).joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE")
+    em_data_kappaR = (ICI_DATA_PATH.joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE")
                     .joinpath("FIX_A").joinpath(f"ECC_{np.round(abar/2, 4)}"))
 
     ic_data_kappaR = []
@@ -407,7 +406,7 @@ def combined_all(config_data: dict, save_as=None):
 
     params = ModelParams(R=150, kappaR=np.asarray(kappaR_list))
 
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=0)
@@ -415,8 +414,7 @@ def combined_all(config_data: dict, save_as=None):
     x_axis_kappaR = path_plot.rot_path.stack_x_axes()
     labels_kappaR = [rf'$\kappa R={kR}$' for kR in [1, 3, 10]]
 
-
-    em_data_abar = (Path(config_data["emanuele_data"]).joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE").joinpath("FIX_M"))
+    em_data_abar = ICI_DATA_PATH.joinpath("FIG_4_Panels_ACE").joinpath("FIXEDCHARGE").joinpath("FIX_M")
 
     ic_data_abar = []
     for ab in abar_list:
@@ -425,7 +423,7 @@ def combined_all(config_data: dict, save_as=None):
 
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex1 = expansion.MappedExpansionQuad(np.asarray(abar_list), params.kappaR, sigma_tilde, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(np.asarray(abar_list), params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=None)
@@ -434,43 +432,52 @@ def combined_all(config_data: dict, save_as=None):
     labels_abar = [rf'$\bar a={a}$' for a in abar_list]
 
 
-    em_data_charge = (Path(config_data["emanuele_data"]).joinpath("FIG_4_Panels_ACE").joinpath("CHARGE_ZC"))
+    em_data_charge = (ICI_DATA_PATH.joinpath("FIG_4_Panels_ACE").joinpath("CHARGE_ZC"))
     undercharged = np.load(em_data_charge.joinpath("ZC_-503").joinpath("pathway.npz"))['arr_0']
     neutral = np.load(em_data_charge.joinpath("ZC_-560").joinpath("pathway.npz"))['arr_0']
     overchargerd = np.load(em_data_charge.joinpath("ZC_-617").joinpath("pathway.npz"))['arr_0']
     ic_data_sigma0 = [undercharged, neutral, overchargerd]
 
     params = ModelParams(R=150, kappaR=kappaR)
-    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=np.asarray(sigma0_list))
+    ex1 = charge_distributions.create_mapped_quad_expansion(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=np.asarray(sigma0_list))
     ex2 = ex1.clone()
 
     path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=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()
@@ -478,11 +485,8 @@ def combined_all(config_data: dict, save_as=None):
 
 def main():
 
-    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-        config_data = json.load(config_file)
-
     # model_comparison(config_data, save_data=False,
-    #     save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('quad_og_comparison.png')
+    #     save_as=FIGURES_PATH.joinpath("final_figures").joinpath('quad_og_comparison.png')
     # )
 
     # kappaR_dependence(np.array([1, 3, 10]), 0.5,
@@ -498,40 +502,40 @@ def main():
     #                   )
 
     # 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')
+    #                     # save_as=FIGURES_PATH.joinpath('quadrupole_distance_dep.png')
     #                     )
 
-    # IC_kappaR_dependence(config_data,
-    #                      save_as=Path(config_data["figures"]).joinpath("Emanuele_data").joinpath('IC_quadrupole_kappaR_dep.png')
+    # IC_kappaR_dependence(
+    #                      save_as=FIGURES_PATH.joinpath("ICi_data").joinpath('IC_quadrupole_kappaR_dep.png')
     #                      )
     #
-    # IC_abar_dependence(config_data, save_as=Path(config_data["figures"]).joinpath("Emanuele_data").
+    # IC_abar_dependence(save_as=FIGURES_PATH.joinpath("ICi_data").
     #                    joinpath('IC_quadrupole_abar_dep.png'))
     #
-    # IC_sigma0_dependence(config_data, save_as=Path(config_data["figures"]).joinpath("Emanuele_data").
+    # IC_sigma0_dependence(save_as=FIGURES_PATH.joinpath("ICi_data").
     #                      joinpath('IC_quadrupole_charge_dep_abar05_kappaR3.png'))
 
-    # combined_kappaR_dependence(config_data, kappaR=[1, 3, 10], abar=0.5,
-    #                      save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('quad_kappaR_dep.png')
+    # combined_kappaR_dependence(kappaR=[1, 3, 10], abar=0.5,
+    #                      save_as=FIGURES_PATH.joinpath("final_figures").joinpath('quad_kappaR_dep.png')
     #                      )
 
-    # combined_sigma0_dependence(config_data,
-    #                      save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('quad_charge_dep.png')
+    # combined_sigma0_dependence(
+    #                      save_as=FIGURES_PATH.joinpath("final_figures").joinpath('quad_charge_dep.png')
     #                      )
 
-    # combined_abar_dependence(config_data, kappaR=3, abar=[0.3, 0.4, 0.5],
-    #                          save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('quad_abar_dep.png')
+    # combined_abar_dependence(kappaR=3, abar=[0.3, 0.4, 0.5],
+    #                          save_as=FIGURES_PATH.joinpath("final_figures").joinpath('quad_abar_dep.png')
     #                          )
 
-    # combined_rescaled_distance_dependence(config_data)
+    # combined_rescaled_distance_dependence()
 
-    combined_distance_dependence(config_data,
-                                 save_as=Path(config_data["figures"]).joinpath("final_figures").joinpath('quad_dist_dep.png')
-                                 )
+    # combined_distance_dependence(
+    #                              save_as=FIGURES_PATH.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(
+                 # save_as=FIGURES_PATH.joinpath("final_figures").joinpath('quad_combined_dep.png')
+                 )
 
 
 if __name__ == '__main__':

analysis/quadrupole_model_mappings.py

@@ -1,4 +1,4 @@
-from charged_shells import patch_size, expansion, parameters
+from charged_shells import patch_size, parameters, charge_distributions
 from charged_shells import functions as fn
 from scipy.special import eval_legendre
 import numpy as np
@@ -19,35 +19,35 @@ def point_to_cap_magnitude(sigma_tilde: Array, a_bar: Array, theta0: Array, kapp
 
 
 def ic_to_gauss(sigma_tilde, a_bar, params: ModelParams, l_max: int = 30,
-                sigma0: float = 0) -> expansion.GaussianCharges:
+                sigma0: float = 0) -> charge_distributions.create_gaussian_charge_expansion:
 
-    ex_mapped = expansion.MappedExpansionQuad(a_bar=a_bar, kappaR=params.kappaR,
-                                              sigma_tilde=sigma_tilde, l_max=30, sigma0=sigma0)
+    ex_mapped = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar, kappaR=params.kappaR,
+                                                                        sigma_tilde=sigma_tilde, l_max=30, sigma0=sigma0)
     target_patch_size = patch_size.potential_patch_size(ex_mapped, params)
-    sigma0_mapped = expansion.net_charge_map(sigma0, params.kappaR)
+    sigma0_mapped = charge_distributions.net_charge_map(sigma0, params.kappaR)
 
     def fn_gauss(x):
-        return expansion.GaussianCharges(lambda_k=x, omega_k=np.array([[0, 0], [np.pi, 0]]),
-                                         sigma1=sigma_tilde, l_max=l_max, sigma0=sigma0_mapped)
+        return charge_distributions.create_gaussian_charge_expansion(lambda_k=x, omega_k=np.array([[0, 0], [np.pi, 0]]),
+                                                                   sigma1=sigma_tilde, l_max=l_max, sigma0=sigma0_mapped)
 
     lbd = patch_size.inverse_potential_patch_size(target_patch_size, fn_gauss, 5, params)
     gauss_sigma = point_to_gauss_magnitude(sigma_tilde, a_bar, lbd, params.kappaR)
-    return expansion.GaussianCharges(lambda_k=lbd, omega_k=np.array([[0, 0], [np.pi, 0]]),
-                                     sigma1=gauss_sigma, l_max=l_max, sigma0=sigma0_mapped)
+    return charge_distributions.create_gaussian_charge_expansion(lambda_k=lbd, omega_k=np.array([[0, 0], [np.pi, 0]]),
+                                                               sigma1=gauss_sigma, l_max=l_max, sigma0=sigma0_mapped)
 
 
-def ic_to_cap(sigma_tilde, a_bar, params: ModelParams, l_max: int = 30, sigma0: float = 0) -> expansion.SphericalCap:
+def ic_to_cap(sigma_tilde, a_bar, params: ModelParams, l_max: int = 30, sigma0: float = 0) -> charge_distributions.create_spherical_cap_expansion:
 
-    ex_mapped = expansion.MappedExpansionQuad(a_bar=a_bar, kappaR=params.kappaR,
-                                              sigma_tilde=sigma_tilde, l_max=30, sigma0=sigma0)
+    ex_mapped = charge_distributions.create_mapped_quad_expansion(a_bar=a_bar, kappaR=params.kappaR,
+                                                                        sigma_tilde=sigma_tilde, l_max=30, sigma0=sigma0)
     target_patch_size = patch_size.potential_patch_size(ex_mapped, params)
-    sigma0_mapped = expansion.net_charge_map(sigma0, params.kappaR)
+    sigma0_mapped = charge_distributions.net_charge_map(sigma0, params.kappaR)
 
     def fn_cap(x):
-        return expansion.SphericalCap(theta0_k=x, sigma1=sigma_tilde, l_max=l_max,
-                                      omega_k=np.array([[0, 0], [np.pi, 0]]), sigma0=sigma0_mapped)
+        return charge_distributions.create_spherical_cap_expansion(theta0_k=x, sigma1=sigma_tilde, l_max=l_max,
+                                                                omega_k=np.array([[0, 0], [np.pi, 0]]), sigma0=sigma0_mapped)
 
     theta0 = patch_size.inverse_potential_patch_size(target_patch_size, fn_cap, 0.5, params)
     cap_sigma = point_to_cap_magnitude(sigma_tilde, a_bar, theta0, params.kappaR)
-    return expansion.SphericalCap(theta0_k=theta0, sigma1=cap_sigma,
-                                  omega_k=np.array([[0, 0], [np.pi, 0]]), l_max=l_max, sigma0=0)
+    return charge_distributions.create_spherical_cap_expansion(theta0_k=theta0, sigma1=cap_sigma,
+                                                            omega_k=np.array([[0, 0], [np.pi, 0]]), l_max=l_max, sigma0=0)

analysis/sEE_minimum.py

@@ -1,10 +1,9 @@
 import matplotlib.pyplot as plt
 import numpy as np
-from charged_shells import expansion, interactions, mapping
+from charged_shells import expansion, interactions, mapping, charge_distributions
 from charged_shells.parameters import ModelParams
 from functools import partial
-from pathlib import Path
-import json
+from config import *
 from matplotlib import cm
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 from labellines import labelLine, labelLines
@@ -41,7 +40,7 @@ def contours():
     kappaR = np.linspace(1, 10, 20)
     a_bar = np.linspace(0.2, 0.6, 20)
     params = ModelParams(R=150, kappaR=kappaR)
-    ex = expansion.MappedExpansionQuad(a_bar[:, None], kappaR[None, :], 0.001)
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar[:, None], kappaR[None, :], 0.001)
     min_energy, min_angle = sEE_minimum(ex, params, match_expansion_axis_to_params=1, accuracy=0.01)
 
     kR_mesh, a_mesh = np.meshgrid(kappaR, a_bar)
@@ -63,14 +62,14 @@ def kappaR_dependence(kappaR, save_as=None, cmap=cm.jet):
     a_bar = np.linspace(0.12, 0.8, 15)
     params = ModelParams(R=150, kappaR=kappaR)
 
-    ex = expansion.MappedExpansionQuad(a_bar[None, :], kappaR[:, None], 0.001)
+    ex = charge_distributions.create_mapped_quad_expansion(a_bar[None, :], kappaR[:, None], 0.001)
     min_energy, min_angle = sEE_minimum(ex, params, match_expansion_axis_to_params=0, accuracy=0.001,
                                         angle_start=0.5, angle_stop=1.)
 
     kappaR_alt = np.array([0.01, 2, 5, 10, 50])
     params_alt = ModelParams(R=150, kappaR=kappaR_alt)
     a_bar_alt = np.linspace(np.min(a_bar) - 0.05, np.max(a_bar) + 0.05, 20)
-    ex2 = expansion.MappedExpansionQuad(a_bar_alt[:, None], kappaR_alt[None, :], 0.001)
+    ex2 = charge_distributions.create_mapped_quad_expansion(a_bar_alt[:, None], kappaR_alt[None, :], 0.001)
     min_energy_alt, min_angle_alt = sEE_minimum(ex2, params_alt, match_expansion_axis_to_params=1, accuracy=0.001,
                                                 angle_start=0.5, angle_stop=1.)
 
@@ -111,9 +110,9 @@ def kappaR_dependence(kappaR, save_as=None, cmap=cm.jet):
     plt.show()
 
 
-def IC_kappaR_dependence(config_data: dict, which_kappa_lines: list,
+def IC_kappaR_dependence(which_kappa_lines: list,
                          save_as=None, cmap=cm.jet, file_suffix=""):
-    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_9"))
+    em_data_path = ICI_DATA_PATH.joinpath("FIG_SUPP_sEE")
     em_data = np.load(em_data_path.joinpath(f"fig9{file_suffix}.npz"))['arr_0']
     # print(em_data.shape)
 
@@ -208,17 +207,14 @@ def IC_kappaR_dependence(config_data: dict, which_kappa_lines: list,
 
 def main():
 
-    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-        config_data = json.load(config_file)
-
     # contours()
 
     # kappaR_dependence(kappaR=np.linspace(0.1, 30, 20),
     #                   save_as=Path(config_data["figures"]).joinpath('sEE_min_kappaR_abar.png')
     #                   )
 
-    IC_kappaR_dependence(config_data, which_kappa_lines=[0.01, 3., 5., 10., 50.], file_suffix="_1",
-                         save_as=Path(config_data["figures"]).joinpath("Emanuele_data").joinpath('IC_sEE_min_kappaR_abar.png')
+    IC_kappaR_dependence(which_kappa_lines=[0.01, 3., 5., 10., 50.], file_suffix="_1",
+                         # save_as=FIGURES_PATH.joinpath("ICi_data").joinpath('IC_sEE_min_kappaR_abar.png')
                          )
 
 

charged_shells/charge_distributions.py

@@ -0,0 +1,198 @@
+from __future__ import annotations
+import copy
+import numpy as np
+from scipy.special import eval_legendre
+from charged_shells import expansion, functions as fn
+import quaternionic
+
+Array = np.ndarray
+Quaternion = quaternionic.array
+Expansion = expansion.Expansion
+
+
+def create_expansion24(sigma2: float | Array, sigma4: float | Array, sigma0: float | Array = 0.):
+    l_array = np.array([0, 2, 4])
+    try:
+        broadcasted_arrs = np.broadcast_arrays(np.sqrt(4*np.pi) * sigma0, sigma2, sigma4)
+    except ValueError:
+        raise ValueError("Given sigma0, sigma2 and sigma4 arrays cannot be broadcast to a common shape.")
+    coefs = rot_sym_expansion(l_array, np.stack(broadcasted_arrs, axis=-1))
+    return Expansion(l_array, coefs)
+
+
+def create_mapped_rotsym_expansion(a_bar: Array | float,
+                                   kappaR: Array | float,
+                                   sigma_tilde: float,
+                                   l_array: Array,
+                                   sigma0: float | Array = 0) -> Expansion:
+    """
+    Create Expansion that matches the outside potential of an impermeable particle with
+    a rotationally symmetric point charge distribution inside (specifically, either quadrupolar or dipolar).
+    :param a_bar: distance between the center and off center charges
+    :param kappaR: screening parameter
+    :param sigma_tilde: magnitude of off-center charges / 4pi R^2
+    :param l_array:  array of necessary multipole numbers
+    :param sigma0: total (mean) charge density
+    """
+    if not isinstance(sigma0, Array):
+        sigma0 = np.array([sigma0])
+    if sigma0.ndim > 1:
+        raise ValueError(f'Sigma0 parameter cannot be an array of dimensions larger than 1, got dim={sigma0.ndim}')
+    a_bar_bc, kappaR_bc = np.broadcast_arrays(a_bar, kappaR)
+
+    a_bar_bc2, kappaR_bc2, l_array_expanded = np.broadcast_arrays(a_bar_bc[..., None],
+                                                                  kappaR_bc[..., None],
+                                                                  l_array[None, :])
+    coefs = (2 * sigma_tilde * fn.coef_C_diff(l_array_expanded, kappaR_bc2)
+             * np.sqrt(4 * np.pi * (2 * l_array_expanded + 1)) * np.power(a_bar_bc2, l_array_expanded))
+    coefs = np.squeeze(rot_sym_expansion(l_array, coefs))
+    kappaR_bc3, sigma0_bc3 = np.broadcast_arrays(kappaR_bc[..., None], sigma0[None, :])
+    coefs = expansion_total_charge(coefs, net_charge_map(sigma0_bc3, kappaR_bc3), l_array)
+    return Expansion(l_array, np.squeeze(coefs))
+
+
+def create_mapped_quad_expansion(a_bar: Array | float,
+                                 kappaR: Array | float,
+                                 sigma_tilde: float,
+                                 l_max: int = 20,
+                                 sigma0: float | Array = 0) -> Expansion:
+    """Create Expansion mapped from a quadrupolar point charge distribution."""
+    l_array = np.array([l for l in range(l_max + 1) if l % 2 == 0])
+    return create_mapped_rotsym_expansion(a_bar, kappaR, sigma_tilde, l_array, sigma0)
+
+
+def create_mapped_dipolar_expansion(a_bar: Array | float,
+                                    kappaR: Array | float,
+                                    sigma_tilde: float,
+                                    l_max: int = 20,
+                                    sigma0: float | Array = 0) -> Expansion:
+    """Create Expansion mapped from a dipolar point charge distribution."""
+    l_array = np.array([l for l in range(l_max + 1) if l % 2 == 1 or l == 0])
+    return create_mapped_rotsym_expansion(a_bar, kappaR, sigma_tilde, l_array, sigma0)
+
+
+def create_gaussian_charge_expansion(omega_k: Array,
+                                     lambda_k: Array | float,
+                                     sigma1: float,
+                                     l_max: int,
+                                     sigma0: float | Array = 0,
+                                     equal_charges: bool = True)  -> Expansion:
+    """
+    Create Expansion for a collection of smeared charges on the sphere.
+    :param omega_k: array of positions (theta, phi) of all charges
+    :param lambda_k: smear parameter for each charge or smear for different cases (if equal_charges = True)
+    :param sigma1: scaling
+    :param l_max: maximal ell value for the expansion
+    :param sigma0: total (mean) charge density
+    :param equal_charges: if this is False, length of lambda_k should be N. If True, theta0_k array will be treated
+        as different expansion cases
+    """
+    omega_k = omega_k.reshape(-1, 2)
+    if not isinstance(lambda_k, Array):
+        lambda_k = np.array([lambda_k])
+    if equal_charges:
+        if lambda_k.ndim > 1:
+            raise ValueError(f'If equal_charges=True, lambda_k should be a 1D array, got shape {lambda_k.shape}')
+        lambda_k = np.full((omega_k.shape[0], lambda_k.shape[0]), lambda_k).T
+    if lambda_k.shape[-1] != omega_k.shape[0]:
+        raise ValueError("Number of charges (length of omega_k) should match the last dimension of lambda_k array.")
+    lambda_k = lambda_k.reshape(-1, omega_k.shape[0])
+    l_array = np.arange(l_max + 1)
+    full_l_array, full_m_array = expansion.full_lm_arrays(l_array)
+    theta_k = omega_k[:, 0]
+    phi_k = omega_k[:, 1]
+    summands = (lambda_k[:, None, :] / np.sinh(lambda_k[:, None, :])
+                * fn.sph_bessel_i(full_l_array[None, :, None], lambda_k[:, None, :])
+                * np.conj(fn.sph_harm(full_l_array[None, :, None], full_m_array[None, :, None],
+                                      theta_k[None, None, :], phi_k[None, None, :])))
+    coefs = np.squeeze(4 * np.pi * sigma1 * np.sum(summands, axis=-1))
+    coefs = expansion_total_charge(coefs, sigma0, l_array)
+    l_array, coefs = purge_unneeded_l(l_array, coefs)
+    return Expansion(l_array, coefs)
+
+
+def create_spherical_cap_expansion(omega_k: Array,
+                                   theta0_k: Array | float,
+                                   sigma1: float,
+                                   l_max: int,
+                                   sigma0: float | Array = 0,
+                                   equal_sizes: bool = True) -> Expansion:
+    """
+    Create Expansion for a collection of spherical caps.
+    :param omega_k: array of positions (theta, phi) of all spherical caps
+    :param theta0_k: sizes of each spherical caps or cap sizes for different cases (if equal_sizes = True)
+    :param sigma1: charge magnitude for the single cap, currently assumes that this is equal for all caps
+    :param l_max: maximal ell value for the expansion
+    :param sigma0: total (mean) charge density
+    :param equal_sizes: if this is False, length of theta0_k should be N. If True, theta0_k array will be treated as
+        different expansion cases
+    """
+    omega_k = omega_k.reshape(-1, 2)
+    if not isinstance(theta0_k, Array):
+        theta0_k = np.array(theta0_k)
+    if equal_sizes:
+        if theta0_k.ndim == 0:
+            theta0_k = np.full(omega_k.shape[0], theta0_k)
+        elif theta0_k.ndim == 1:
+            theta0_k = np.full((omega_k.shape[0], theta0_k.shape[0]), theta0_k)
+        else:
+            raise ValueError(f'If equal_charges=True, theta0_k should be a 1D array, got shape {theta0_k.shape}')
+    if theta0_k.shape[0] != omega_k.shape[0]:
+        raise ValueError("Number of charges (length of omega_k) should match the last dimension of theta0_k array.")
+    rotations = Quaternion(np.stack((np.cos(omega_k[..., 0] / 2),
+                                     np.sin(omega_k[..., 1]) * np.sin(omega_k[..., 0] / 2),
+                                     np.cos(omega_k[..., 1]) * np.sin(omega_k[..., 0] / 2),
+                                     np.zeros_like(omega_k[..., 0]))).T)
+    l_array = np.arange(l_max + 1)
+    coefs_l0 = -sigma1 * (np.sqrt(np.pi / (2 * l_array[None, :] + 1)) *
+                          (eval_legendre(l_array[None, :] + 1, np.cos(theta0_k[..., None]))
+                           - eval_legendre(l_array[None, :] - 1, np.cos(theta0_k[..., None]))))
+    coefs = rot_sym_expansion(l_array, coefs_l0)
+    coefs_all_single_caps = expansion.expansion_rotation(rotations, coefs, l_array)
+    # Rotating is implemented in such a way that it rotates every patch to every position,
+    # hence the redundancy of out of diagonal elements.
+    coefs_all = np.sum(np.diagonal(coefs_all_single_caps), axis=-1)
+    coefs_all = expansion_total_charge(coefs_all, sigma0, l_array)
+    return Expansion(l_array, np.squeeze(coefs_all))
+
+
+def net_charge_map(sigma0: float | Array, kappaR: float | Array):
+    return sigma0 * np.exp(kappaR) / np.sinh(kappaR) * kappaR / (1 + kappaR)
+
+
+def rot_sym_expansion(l_array: Array, coefs: Array) -> Array:
+    """Create full expansion array for rotationally symmetric distributions with only m=0 terms different form 0."""
+    full_coefs = np.zeros(coefs.shape[:-1] + (np.sum(2 * l_array + 1),))
+    full_coefs[..., np.cumsum(2 * l_array + 1) - l_array - 1] = coefs
+    return full_coefs
+
+
+def expansion_total_charge(coefs: Array, sigma0: float | Array | None, l_vals: Array):
+    """Adds a new axis to the expansion coefficients that modifies expansion based on given net charge density."""
+    if l_vals[0] != 0:
+        raise NotImplementedError("To modify total charge, the first expansion term should be l=0. "
+                                  "Adding l=0 term not yet implemented (TODO).")
+    if sigma0 is None:
+        return coefs
+    if not isinstance(sigma0, Array):
+        x = copy.deepcopy(coefs)
+        x[..., 0] = sigma0 * np.sqrt(4 * np.pi)
+        return x
+
+    # insert new axis in the 2nd-to-last place and repeat the expansion data over this new axis
+    x = np.repeat(np.expand_dims(coefs, -2), sigma0.shape[-1], axis=-2)
+    x[..., 0] = sigma0 * np.sqrt(4 * np.pi)
+    return x
+
+
+def purge_unneeded_l(l_array: Array, coefs: Array) -> (Array, Array):
+    """Remove ell values from expansion for which all (ell, m) coefficients are zero."""
+    def delete_zero_entries(l, l_arr, cfs):
+        l_idx = np.where(l_arr == l)[0][0]
+        m_indices = expansion.m_indices_at_l(l_arr, l_idx)
+        if np.all(cfs[..., m_indices] == 0):
+            return np.delete(l_arr, l_idx), np.delete(cfs, m_indices, axis=-1)
+        return l_arr, cfs
+    for l in l_array:
+        l_array, coefs = delete_zero_entries(l, l_array, coefs)
+    return l_array, coefs

charged_shells/expansion.py

@@ -5,8 +5,6 @@ import charged_shells.functions as fn
 import quaternionic
 import spherical
 import copy
-from scipy.special import eval_legendre
-
 
 Array = np.ndarray
 Quaternion = quaternionic.array
@@ -23,7 +21,7 @@ class Expansion:
     l_array: Array
     coefs: Array
     _starting_coefs: Array = None  # initialized with the __post_init__ method
-    _rotations: Quaternion = Quaternion([1., 0., 0., 0.])
+    _rotations: Quaternion = None
 
     def __post_init__(self):
         if self.coefs.shape[-1] != np.sum(2 * self.l_array + 1):
@@ -32,6 +30,7 @@ class Expansion:
             raise InvalidExpansion('Array of l values should be unique and sorted.')
         self.coefs = self.coefs.astype(np.complex128)
         self._starting_coefs = np.copy(self.coefs)
+        self._rotations = Quaternion([1., 0., 0., 0.])
 
     def __getitem__(self, item):
         return Expansion(self.l_array, self.coefs[item])
@@ -71,8 +70,6 @@ class Expansion:
         self.coefs = expansion_rotation(rotations, self._starting_coefs, self.l_array)
 
     def rotate_euler(self, alpha: Array = 0, beta: Array = 0, gamma: Array = 0, rotate_existing=False):
-        # TODO: additional care required on the convention used to transform euler angles to quaternions
-        # TODO: might be off for a minus sign for each? angle !!
         R_euler = quaternionic.array.from_euler_angles(alpha, beta, gamma)
         self.rotate(R_euler, rotate_existing=rotate_existing)
 
@@ -100,161 +97,12 @@ class Expansion:
         return copy.deepcopy(self)
 
 
-class Expansion24(Expansion):
-
-    def __init__(self, sigma2: float | Array, sigma4: float | Array, sigma0: float | Array = 0.):
-        l_array = np.array([0, 2, 4])
-        try:
-            broadcasted_arrs = np.broadcast_arrays(np.sqrt(4*np.pi) * sigma0, sigma2, sigma4)
-        except ValueError:
-            raise ValueError("Given sigma0, sigma2 and sigma4 arrays cannot be broadcast to a common shape.")
-        coefs = rot_sym_expansion(l_array, np.stack(broadcasted_arrs, axis=-1))
-        super().__init__(l_array, coefs)
-
-
-class MappedExpansionSym(Expansion):
+def m_indices_at_l(l_arr: Array, l_idx: int):
     """
-    Expansion that matches the outside potential of an impermeable particle with
-    a rotationally symmetric point charge distribution inside (specifically, either quadrupolar or dipolar).
+    For a given l_array and index l_idx for some ell in this array, get indices of all (ell, m) coefficients
+    in coefficients array.
     """
-
-    def __init__(self,
-                 a_bar: Array | float,
-                 kappaR: Array | float,
-                 sigma_tilde: float,
-                 l_array: Array,
-                 sigma0: float | Array = 0):
-        """
-        :param a_bar: distance between the center and off center charges
-        :param kappaR: screening parameter
-        :param sigma_tilde: magnitude of off-center charges / 4pi R^2
-        :param l_max: maximal ell value for the expansion
-        :param sigma0: total (mean) charge density
-        """
-        if not isinstance(sigma0, Array):
-            sigma0 = np.array([sigma0])
-        if sigma0.ndim > 1:
-            raise ValueError(f'Sigma0 parameter cannot be an array of dimensions larger than 1, got dim={sigma0.ndim}')
-        a_bar_bc, kappaR_bc = np.broadcast_arrays(a_bar, kappaR)
-
-        a_bar_bc2, kappaR_bc2, l_array_expanded = np.broadcast_arrays(a_bar_bc[..., None],
-                                                                      kappaR_bc[..., None],
-                                                                      l_array[None, :])
-        coefs = (2 * sigma_tilde * fn.coef_C_diff(l_array_expanded, kappaR_bc2)
-                 * np.sqrt(4 * np.pi * (2 * l_array_expanded + 1)) * np.power(a_bar_bc2, l_array_expanded))
-        coefs = np.squeeze(rot_sym_expansion(l_array, coefs))
-        kappaR_bc3, sigma0_bc3 = np.broadcast_arrays(kappaR_bc[..., None], sigma0[None, :])
-        coefs = expansion_total_charge(coefs, net_charge_map(sigma0_bc3, kappaR_bc3), l_array)
-        super().__init__(l_array, np.squeeze(coefs))
-
-
-class MappedExpansionQuad(MappedExpansionSym):
-
-    def __init__(self,
-                 a_bar: Array | float,
-                 kappaR: Array | float,
-                 sigma_tilde: float,
-                 l_max: int = 20,
-                 sigma0: float | Array = 0):
-        l_array = np.array([l for l in range(l_max + 1) if l % 2 == 0])
-        super().__init__(a_bar, kappaR, sigma_tilde, l_array, sigma0)
-
-
-class MappedExpansionDipole(MappedExpansionSym):
-
-    def __init__(self,
-                 a_bar: Array | float,
-                 kappaR: Array | float,
-                 sigma_tilde: float,
-                 l_max: int = 20,
-                 sigma0: float | Array = 0):
-        l_array = np.array([l for l in range(l_max + 1) if l % 2 == 1 or l == 0])
-        super().__init__(a_bar, kappaR, sigma_tilde, l_array, sigma0)
-
-
-class GaussianCharges(Expansion):
-    """Expansion for a collection of smeared charges on the sphere."""
-
-    def __init__(self, omega_k: Array, lambda_k: Array | float, sigma1: float, l_max: int,
-                 sigma0: float | Array = 0, equal_charges: bool = True):
-        """
-        :param omega_k: array of positions (theta, phi) of all charges
-        :param lambda_k: smear parameter for each charge or smear for different cases (if equal_charges = True)
-        :param sigma1: scaling
-        :param l_max: maximal ell value for the expansion
-        :param sigma0: total (mean) charge density
-        :param equal_charges: if this is False, length of lambda_k should be N. If True, theta0_k array will be treated
-            as different expansion cases
-        """
-        omega_k = omega_k.reshape(-1, 2)
-        if not isinstance(lambda_k, Array):
-            lambda_k = np.array([lambda_k])
-        if equal_charges:
-            if lambda_k.ndim > 1:
-                raise ValueError(f'If equal_charges=True, lambda_k should be a 1D array, got shape {lambda_k.shape}')
-            lambda_k = np.full((omega_k.shape[0], lambda_k.shape[0]), lambda_k).T
-        if lambda_k.shape[-1] != omega_k.shape[0]:
-            raise ValueError("Number of charges (length of omega_k) should match the last dimension of lambda_k array.")
-        lambda_k = lambda_k.reshape(-1, omega_k.shape[0])
-        l_array = np.arange(l_max + 1)
-        full_l_array, full_m_array = full_lm_arrays(l_array)
-        theta_k = omega_k[:, 0]
-        phi_k = omega_k[:, 1]
-        summands = (lambda_k[:, None, :] / np.sinh(lambda_k[:, None, :])
-                    * fn.sph_bessel_i(full_l_array[None, :, None], lambda_k[:, None, :])
-                    * np.conj(fn.sph_harm(full_l_array[None, :, None], full_m_array[None, :, None],
-                                          theta_k[None, None, :], phi_k[None, None, :])))
-        coefs = np.squeeze(4 * np.pi * sigma1 * np.sum(summands, axis=-1))
-        coefs = expansion_total_charge(coefs, sigma0, l_array)
-        l_array, coefs = purge_unneeded_l(l_array, coefs)
-        super().__init__(l_array, coefs)
-
-
-class SphericalCap(Expansion):
-    """Expansion for a collection of spherical caps."""
-
-    def __init__(self, omega_k: Array, theta0_k: Array | float, sigma1: float, l_max: int, sigma0: float | Array = 0,
-                 equal_sizes: bool = True):
-        """
-        :param omega_k: array of positions (theta, phi) of all spherical caps
-        :param theta0_k: sizes of each spherical caps or cap sizes for different cases (if equal_sizes = True)
-        :param sigma1: charge magnitude for the single cap, currently assumes that this is equal for all caps
-        :param l_max: maximal ell value for the expansion
-        :param sigma0: total (mean) charge density
-        :param equal_sizes: if this is False, length of theta0_k should be N. If True, theta0_k array will be treated as
-            different expansion cases
-        """
-        omega_k = omega_k.reshape(-1, 2)
-        if not isinstance(theta0_k, Array):
-            theta0_k = np.array(theta0_k)
-        if equal_sizes:
-            if theta0_k.ndim == 0:
-                theta0_k = np.full(omega_k.shape[0], theta0_k)
-            elif theta0_k.ndim == 1:
-                theta0_k = np.full((omega_k.shape[0], theta0_k.shape[0]), theta0_k)
-            else:
-                raise ValueError(f'If equal_charges=True, theta0_k should be a 1D array, got shape {theta0_k.shape}')
-        if theta0_k.shape[0] != omega_k.shape[0]:
-            raise ValueError("Number of charges (length of omega_k) should match the last dimension of theta0_k array.")
-        rotations = Quaternion(np.stack((np.cos(omega_k[..., 0] / 2),
-                                         np.sin(omega_k[..., 1]) * np.sin(omega_k[..., 0] / 2),
-                                         np.cos(omega_k[..., 1]) * np.sin(omega_k[..., 0] / 2),
-                                         np.zeros_like(omega_k[..., 0]))).T)
-        l_array = np.arange(l_max + 1)
-        coefs_l0 = -sigma1 * (np.sqrt(np.pi / (2 * l_array[None, :] + 1)) *
-                              (eval_legendre(l_array[None, :] + 1, np.cos(theta0_k[..., None]))
-                               - eval_legendre(l_array[None, :] - 1, np.cos(theta0_k[..., None]))))
-        coefs = rot_sym_expansion(l_array, coefs_l0)
-        coefs_all_single_caps = expansion_rotation(rotations, coefs, l_array)
-        # Rotating is implemented in such a way that it rotates every patch to every position,
-        # hence the redundancy of out of diagonal elements.
-        coefs_all = np.sum(np.diagonal(coefs_all_single_caps), axis=-1)
-        coefs_all = expansion_total_charge(coefs_all, sigma0, l_array)
-        super().__init__(l_array, np.squeeze(coefs_all))
-
-
-def net_charge_map(sigma0: float | Array, kappaR: float | Array):
-    return sigma0 * np.exp(kappaR) / np.sinh(kappaR) * kappaR / (1 + kappaR)
+    return np.arange(np.sum(2 * l_arr[:l_idx] + 1), np.sum(2 * l_arr[:l_idx+1] + 1))
 
 
 def full_lm_arrays(l_array: Array) -> (Array, Array):
@@ -266,52 +114,6 @@ def full_lm_arrays(l_array: Array) -> (Array, Array):
     return np.repeat(l_array, 2 * l_array + 1), np.array(all_m_list)
 
 
-def rot_sym_expansion(l_array: Array, coefs: Array) -> Array:
-    """Create full expansion array for rotationally symmetric distributions with only m=0 terms different form 0."""
-    full_coefs = np.zeros(coefs.shape[:-1] + (np.sum(2 * l_array + 1),))
-    full_coefs[..., np.cumsum(2 * l_array + 1) - l_array - 1] = coefs
-    return full_coefs
-
-
-def expansion_total_charge(coefs: Array, sigma0: float | Array | None, l_vals: Array):
-    """Adds a new axis to the expansion coefficients that modifies expansion based on given net charge density."""
-    if l_vals[0] != 0:
-        raise NotImplementedError("To modify total charge, the first expansion term should be l=0. "
-                                  "Adding l=0 term not yet implemented (TODO).")
-    if sigma0 is None:
-        return coefs
-    if not isinstance(sigma0, Array):
-        x = copy.deepcopy(coefs)
-        x[..., 0] = sigma0 * np.sqrt(4 * np.pi)
-        return x
-
-    # insert new axis in the 2nd-to-last place and repeat the expansion data over this new axis
-    x = np.repeat(np.expand_dims(coefs, -2), sigma0.shape[-1], axis=-2)
-    x[..., 0] = sigma0 * np.sqrt(4 * np.pi)
-    return x
-
-
-def m_indices_at_l(l_arr: Array, l_idx: int):
-    """
-    For a given l_array and index l_idx for some ell in this array, get indices of all (ell, m) coefficients
-    in coefficients array.
-    """
-    return np.arange(np.sum(2 * l_arr[:l_idx] + 1), np.sum(2 * l_arr[:l_idx+1] + 1))
-
-
-def purge_unneeded_l(l_array: Array, coefs: Array) -> (Array, Array):
-    """Remove ell values from expansion for which all (ell, m) coefficients are zero."""
-    def delete_zero_entries(l, l_arr, cfs):
-        l_idx = np.where(l_arr == l)[0][0]
-        m_indices = m_indices_at_l(l_arr, l_idx)
-        if np.all(cfs[..., m_indices] == 0):
-            return np.delete(l_arr, l_idx), np.delete(cfs, m_indices, axis=-1)
-        return l_arr, cfs
-    for l in l_array:
-        l_array, coefs = delete_zero_entries(l, l_array, coefs)
-    return l_array, coefs
-
-
 def coefs_fill_missing_l(expansion: Expansion, target_l_array: Array) -> Expansion:
     """Explicitly add missing expansion coefficients so that expansion includes all ell values from the target array."""
     missing_l = np.setdiff1d(target_l_array, expansion.l_array, assume_unique=True)
@@ -327,7 +129,7 @@ def coefs_fill_missing_l(expansion: Expansion, target_l_array: Array) -> Expansi
 def expansions_to_common_l(ex1: Expansion, ex2: Expansion) -> (Expansion, Expansion):
     """Explicitly add zero expansion coefficients so that both expansions include coefficients for the same ell."""
     common_l_array = np.union1d(ex1.l_array, ex2.l_array)
-    return coefs_fill_missing_l(ex1, common_l_array),  coefs_fill_missing_l(ex2, common_l_array)
+    return coefs_fill_missing_l(ex1, common_l_array), coefs_fill_missing_l(ex2, common_l_array)
 
 
 def expansion_rotation(rotations: Quaternion, coefs: Array, l_array: Array):

charged_shells/interactions.py

@@ -5,7 +5,6 @@ import numpy as np
 from typing import Literal
 import charged_shells.units_and_constants as uc
 
-
 Array = np.ndarray
 Expansion = expansion.Expansion
 ModelParams = parameters.ModelParams
@@ -23,7 +22,8 @@ def energy_units(units: EnergyUnit, params: ModelParams) -> float:
             return uc.UNITS.energy
 
 
-def charged_shell_energy(ex1: Expansion, ex2: Expansion, params: ModelParams, dist: float = 2, units: EnergyUnit = 'kT'):
+def charged_shell_energy(ex1: Expansion, ex2: Expansion, params: ModelParams, dist: float = 2, units: EnergyUnit = 'kT',
+                         chunk_size: int = 1000):
 
     ex1, ex2 = expansion.expansions_to_common_l(ex1, ex2)
     dist = dist * params.R
@@ -52,36 +52,50 @@ def charged_shell_energy(ex1: Expansion, ex2: Expansion, params: ModelParams, di
     # additional selection that excludes terms where Wigner 3j symbols are 0 by definition
     s_bool1 = np.abs(flat_l[:, None] - all_s_array[None, :]) <= flat_p[:, None]
     s_bool2 = flat_p[:, None] <= (flat_l[:, None] + all_s_array[None, :])
-    indices_lpm, indices_s = np.nonzero(s_bool1 * s_bool2)
+    indices_lpm_all, indices_s_all = np.nonzero(s_bool1 * s_bool2)
+
+    # indices array can get really large (a lot of combinations) so we split the calculation into chunks to preserve RAM
+    # interestingly, this also leads to performance improvements if chunks are still large enough
+    if chunk_size is None:
+        chunk_size = len(indices_lpm_all)
+    num_sections = np.ceil(len(indices_lpm_all) / chunk_size)
+    energy = 0
+    for indices_lpm, indices_s in zip(np.array_split(indices_lpm_all, num_sections),
+                                      np.array_split(indices_s_all, num_sections)):
+
+        l_vals = flat_l[indices_lpm]
+        p_vals = flat_p[indices_lpm]
+        m_vals = flat_m[indices_lpm]
+        s_vals = all_s_array[indices_s]
+        bessel_vals = bessels[indices_s]
+
+        # While all other arrays are 1D, this one can have extra leading axes corresponding to leading dimensions
+        # of expansion coefficients. The last dimension is the same as other arrays.
+        charge_vals = charge_factor[..., indices_lpm]
+
+        lps_terms = (2 * s_vals + 1) * np.sqrt((2 * l_vals + 1) * (2 * p_vals + 1))
+
+        # the same combination of l, s, and p is repeated many times
+        _, unique_indices1, inverse1 = np.unique(np.stack((l_vals, s_vals, p_vals)),
+                                                 axis=1, return_inverse=True, return_index=True)
+        wigner1 = wigner3j(2 * l_vals[unique_indices1], 2 * s_vals[unique_indices1], 2 * p_vals[unique_indices1],
+                           0, 0, 0)[inverse1]
 
-    l_vals = flat_l[indices_lpm]
-    p_vals = flat_p[indices_lpm]
-    m_vals = flat_m[indices_lpm]
-    s_vals = all_s_array[indices_s]
-    bessel_vals = bessels[indices_s]
+        # all the combinations (l, s, p, m) are unique
+        wigner2 = wigner3j(2 * l_vals, 2 * s_vals, 2 * p_vals,
+                           -2 * m_vals, 0, 2 * m_vals)
 
-    # While all other arrays are 1D, this one can have extra leading axes corresponding to leading dimensions
-    # of expansion coefficients. The last dimension is the same as other arrays.
-    charge_vals = charge_factor[..., indices_lpm]
+        constants = params.R ** 2 / (params.kappa * params.epsilon * uc.CONSTANTS.epsilon0) * energy_units(units, params)
 
-    lps_terms = (2 * s_vals + 1) * np.sqrt((2 * l_vals + 1) * (2 * p_vals + 1))
+        C_vals = fn.interaction_coef_C(l_vals, p_vals, params.kappaR)
+        lspm_vals = C_vals * (-1) ** (l_vals + m_vals) * lps_terms * bessel_vals * wigner1 * wigner2
+        broadcasted_lspm_vals = np.broadcast_to(lspm_vals, charge_vals.shape)
 
-    # the same combination of l, s, and p is repeated many times
-    _, unique_indices1, inverse1 = np.unique(np.stack((l_vals, s_vals, p_vals)),
-                                             axis=1, return_inverse=True, return_index=True)
-    wigner1 = wigner3j(2 * l_vals[unique_indices1], 2 * s_vals[unique_indices1], 2 * p_vals[unique_indices1],
-                       0, 0, 0)[inverse1]
+        rescale_at_equal_lp = np.where(l_vals == p_vals, 0.5, 1)
 
-    # all the combinations (l, s, p, m) are unique
-    wigner2 = wigner3j(2 * l_vals, 2 * s_vals, 2 * p_vals,
-                       -2 * m_vals, 0, 2 * m_vals)
+        energy +=  constants * np.sum(rescale_at_equal_lp * broadcasted_lspm_vals * charge_vals, axis=-1)
 
-    constants = params.R ** 2 / (params.kappa * params.epsilon * uc.CONSTANTS.epsilon0) * energy_units(units, params)
+    return energy
 
-    C_vals = fn.interaction_coef_C(l_vals, p_vals, params.kappaR)
-    lspm_vals = C_vals * (-1) ** (l_vals + m_vals) * lps_terms * bessel_vals * wigner1 * wigner2
-    broadcasted_lspm_vals = np.broadcast_to(lspm_vals, charge_vals.shape)
 
-    rescale_at_equal_lp = np.where(l_vals == p_vals, 0.5, 1)
 
-    return constants * np.sum(rescale_at_equal_lp * broadcasted_lspm_vals * charge_vals, axis=-1)

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,6 +0,0 @@
-{
-  "emanuele_data":  "/home/andraz/ChargedShells/data_Emanuele/2024-06-24",
-  "figures": "/home/andraz/ChargedShells/Figures",
-  "figure_data": "/home/andraz/ChargedShells/Figures/Data",
-  "expansion_data": "/home/andraz/ChargedShells/ExpansionData/"
-}

requirements.txt

@@ -0,0 +1,6 @@
+numpy>=2.0.2
+py3nj>=0.2.1
+quaternionic>=1.0.12
+scipy>=1.14.1
+setuptools>=75.1.0
+spherical>=1.0.14

setup.py

@@ -0,0 +1,16 @@
+from setuptools import setup, find_packages
+import os
+
+def read_requirements():
+    reqs_path = os.path.join(os.path.dirname(__file__), 'requirements.txt')
+    with open(reqs_path, 'r') as reqs_file:
+        return reqs_file.read().splitlines()
+
+setup(
+    name="charged_shells",
+    version="0.1.0",
+    description="A description of my package",
+    packages=find_packages(),
+    install_requires=read_requirements(),
+    python_requires='>=3.10',
+)

tests/expansion_mapping_test.py

@@ -1,6 +1,6 @@
 import unittest
 import numpy as np
-from charged_shells import expansion, interactions, parameters, units_and_constants, mapping, potentials
+from charged_shells import interactions, parameters, units_and_constants, mapping, potentials, charge_distributions
 from functools import partial
 
 
@@ -15,7 +15,7 @@ class IsotropicTest(unittest.TestCase):
         self.radius = 150
         self.dist = 2 * self.radius
         self.sigma0 = self.charge / (4 * np.pi * self.radius ** 2)
-        self.ex1 = expansion.MappedExpansionQuad(0, self.kappaR, 0, l_max=10, sigma0=self.sigma0)
+        self.ex1 = charge_distributions.create_mapped_quad_expansion(0, self.kappaR, 0, l_max=10, sigma0=self.sigma0)
         self.ex2 = self.ex1.clone()
         self.params = parameters.ModelParams(kappaR=self.kappaR, R=self.radius)
 

tests/interactions_test.py

@@ -1,4 +1,4 @@
-from charged_shells import expansion, interactions
+from charged_shells import interactions, charge_distributions
 from charged_shells.parameters import ModelParams
 import time
 import numpy as np
@@ -8,7 +8,7 @@ import matplotlib.pyplot as plt
 def v22_distance_test():
 
     params = ModelParams(R=10, kappaR=3.29)
-    ex0 = expansion.Expansion24(1, 0, 0)
+    ex0 = charge_distributions.create_expansion24(1, 0, 0)
     ex1 = ex0.clone()
 
     ex0.rotate_euler(0, np.array([0, 0, np.pi / 2]), 0)
@@ -29,7 +29,7 @@ def quadrupole_variation_test():
 
     params = ModelParams(R=10, kappaR=3.29)
     sigma2 = np.array([0.45, 0.5, 0.55, 0.6, 0.65])
-    ex0 = expansion.Expansion24(sigma2, 0, sigma0=0.1)
+    ex0 = charge_distributions.create_expansion24(sigma2, 0, sigma0=0.1)
     ex1 = ex0.clone()
 
     ex1.rotate_euler(0, np.pi / 2, 0)
@@ -45,7 +45,7 @@ def quadrupole_variation_test():
 
 def timing():
     params = ModelParams(R=150, kappaR=3)
-    ex1 = expansion.MappedExpansionQuad(0.44, params.kappaR, 0.001, l_max=20)
+    ex1 = charge_distributions.create_mapped_quad_expansion(0.44, params.kappaR, 0.001, l_max=20)
     ex2 = ex1.clone()
 
     dist = 2.

tests/janus_interaction_test.py

@@ -1,11 +1,11 @@
-from charged_shells import expansion, interactions, parameters
+from charged_shells import charge_distributions, interactions, parameters
 import numpy as np
 
 
 def charge_energy_test():
 
     params = parameters.ModelParams(R=150, kappaR=3)
-    ex1 = expansion.MappedExpansionDipole(0.5, params.kappaR, 0.001, l_max=30, sigma0=0.0002)
+    ex1 = charge_distributions.create_mapped_dipolar_expansion(0.5, params.kappaR, 0.001, l_max=30, sigma0=0.0002)
     ex2 = ex1.clone()
     ex3 = ex1.clone().inverse_sign(exclude_00=True)
 

tests/potential_tests.py

@@ -1,4 +1,4 @@
-from charged_shells import expansion, potentials
+from charged_shells import charge_distributions, potentials
 from charged_shells.parameters import ModelParams
 import numpy as np
 import matplotlib.pyplot as plt
@@ -6,7 +6,7 @@ import matplotlib.pyplot as plt
 
 def main():
     params = ModelParams(R=150, kappaR=3)
-    ex1 = expansion.MappedExpansionQuad(0.44, params.kappaR, 0.001, l_max=30)
+    ex1 = charge_distributions.create_mapped_quad_expansion(0.44, params.kappaR, 0.001, l_max=30)
     # ex2 = expansion.SphericalCap(np.array([[0, 0], [np.pi, 0]]), 0.5, 0.003, l_max=70)
 
     theta = np.linspace(0, np.pi, 1000)