charged_shells is now a package, with a large refactor of expansion construction using factory functions. Analysis now handles config loading in a separate python file.

analysis/ICi_ICp_term_comparison.py

@@ -1,7 +1,6 @@
 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
@@ -122,11 +121,11 @@ def plot_kappaR_dep_ratio(l: int, abar: float, l_diff: int = 2):
 if __name__ == '__main__':
 
     # plot_exp([1, 10, 50], 0.4,
-    #          # save_as='/home/andraz/ChargedShells/Figures/ici_icp_expansion_a08_kr10.png'
+    #          # save_as=FIGURES_PATH.joinpath('ici_icp_expansion_a08_kr10.png')
     #          )
 
     plot_exp2(1, [0.5, 0.8], legend=True,
-             save_as='/home/andraz/ChargedShells/Figures/ici_icp_expansion_kr1.png'
+             # save_as=FIGURES_PATH.joinpath('ici_icp_expansion_kr1.png')
              )
 
     # plot_kappaR_dep_ratio(2, 0.8, l_diff=2)

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,12 +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 expansion, 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
@@ -54,7 +51,7 @@ 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)
@@ -65,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)
@@ -78,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)
@@ -91,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)
@@ -101,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()
 
@@ -138,9 +135,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_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")
+    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):]
@@ -166,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 = []
@@ -181,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()
 
@@ -206,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 = []
@@ -220,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()
 
@@ -245,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):]
@@ -262,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)
 
@@ -287,7 +284,7 @@ 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.00099
     kappaR_list = [1, 3, 10]
@@ -296,7 +293,7 @@ def combined_all(config_data: dict, save_as=None):
     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 = []
@@ -309,7 +306,7 @@ 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(exclude_00=True)
 
@@ -320,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 = []
@@ -334,7 +331,7 @@ 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(exclude_00=True)
 
@@ -345,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']
@@ -361,7 +358,7 @@ def combined_all(config_data: dict, save_as=None):
     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(exclude_00=True)
 
@@ -407,7 +404,7 @@ def combined_all(config_data: dict, save_as=None):
     for ax in axs:
         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()
@@ -415,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,
@@ -432,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("Emanuele_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('Emanuele_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

@@ -5,9 +5,13 @@ Expansion = expansion.Expansion
 RotConfig = Literal['ep', 'pp', 'sep']
 
 
-def get_kappaR_higher_multipole_dicts(peak: Peak, params: ModelParams, emanuele_data: Array,
+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,
@@ -27,8 +31,9 @@ def get_kappaR_higher_multipole_dicts(peak: Peak, params: ModelParams, emanuele_
         d = np.around(distances[i], 5)
         if d in np.around(dist_cs, 5):
             idx, = np.nonzero(inverse == i)
-            kR = emanuele_data[idx, 1][(emanuele_data[idx, 1] <= max_kappaR) * (min_kappa <= emanuele_data[idx, 1])]
-            en = emanuele_data[idx, peak.emanuele_data_column-1][(emanuele_data[idx, 1] <= max_kappaR) * (min_kappa <= emanuele_data[idx, 1])]
+            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.emanuele_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
@@ -41,10 +46,9 @@ def get_kappaR_higher_multipole_dicts(peak: Peak, params: ModelParams, emanuele_
     return data_dict_ic, data_dict_cs
 
 
-def IC_peak_energy_plot(config_data: dict,
-                        dist: list,
+def IC_peak_energy_plot(dist: list,
                         which: RotConfig,
-                        a_bar=0.8,
+                        abar=0.8,
                         min_kappaR: float = 0.01,
                         max_kappaR: float = 50.,
                         R: float = 150,
@@ -53,9 +57,9 @@ 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("TEST_HIGHER_MULTIPOLES"))
-    em_data = np.load(em_data_path.joinpath("higher_multiplot_energy.npz"))
+    # 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']
@@ -65,7 +69,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(a_bar, params.kappaR, 0.00099, l_max=20)
+    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)
@@ -76,15 +80,15 @@ def IC_peak_energy_plot(config_data: dict,
     else:
         raise ValueError
 
-    data_dict_ic2, data_dict_cs = get_kappaR_higher_multipole_dicts(peak, params, data2, kappaR, dist, max_kappaR, min_kappaR)
-    data_dict_ic4, _ = get_kappaR_higher_multipole_dicts(peak, params, data4, kappaR, dist, max_kappaR, min_kappaR)
-    data_dict_ic6, _ = get_kappaR_higher_multipole_dicts(peak, params, data6, kappaR, dist, max_kappaR, min_kappaR)
+    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, 1.8))
+    fig = plt.figure(figsize=(2, 2))
     gs = gridspec.GridSpec(1, 1, figure=fig)
-    gs.update(left=0.25, right=0.95, top=0.97, bottom=0.21)
+    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)
@@ -93,12 +97,13 @@ def IC_peak_energy_plot(config_data: dict,
         ax.plot(data4[:, 0], data4[:, 1], ls='-', c=next(colors), label=r'$l=4$')
         ax.plot(data6[:, 0], data6[:, 1], ls='-', c=next(colors), label=r'$l=6$')
     ax.legend(fontsize=9, ncol=1, frameon=False, handlelength=0.7, loc='upper right',
-              bbox_to_anchor=(0.6, 0.5),
+              bbox_to_anchor=(0.5, 0.5),
               # bbox_to_anchor=(0.42, 0.9)
               )
-    ax.set_title(rf'$\rho = {key:.1f}$', loc='right',
+    ax.set_title(f'$\\bar a = {abar:.1f}$, $\\rho = {key:.1f}$', loc='center',
                  # x=0.95, y=0.8,
-                 x=0.45, y=0.6,
+                 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)
@@ -114,13 +119,10 @@ def IC_peak_energy_plot(config_data: dict,
 
 
 if __name__ == '__main__':
-    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-        config_data = json.load(config_file)
 
-    dist = [2.5]
-    IC_peak_energy_plot(config_data, dist=dist, which='ep', min_kappaR=0.1, max_kappaR=50,
-                        save_as=Path(
-                            '/home/andraz/ChargedShells/Figures/Emanuele_data/peak_ep_kappaR_higher_correction_rho25.png'),
+    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("Emanuele_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,
@@ -248,7 +248,7 @@ def ic_cs_comparison2(save_data=False):
     # plt.tight_layout()
     plt.subplots_adjust(left=0.22, right=0.97, top=0.95, bottom=0.22)
     ax.xaxis.set_label_coords(0.5, -0.17)
-    plt.savefig(Path("/home/andraz/ChargedShells/Figures/final_figures/potential_ic_cs_comparison.png"), dpi=300)
+    plt.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,9 +24,7 @@ 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=(2, 1.7))
@@ -45,14 +42,14 @@ def plot_abar_dependence(*, save=False, save_data=False):
     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)
 
@@ -67,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:
@@ -91,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()
 
 
@@ -102,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)
@@ -121,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()
 
 
@@ -132,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)
@@ -151,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,11 +1,10 @@
 from matplotlib import gridspec
-from charged_shells import expansion, interactions, mapping, functions
+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
 
@@ -156,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.,
@@ -173,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']
 
@@ -182,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.00099, 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)
@@ -235,7 +231,7 @@ def IC_peak_energy_plot(config_data: dict,
     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,
@@ -246,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']
 
@@ -255,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)
@@ -293,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,
@@ -303,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']
 
@@ -312,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)
@@ -350,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']
 
@@ -367,7 +361,7 @@ 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()
+    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)
@@ -418,14 +412,13 @@ def IC_peak_energy_kappaR_combined_plot(config_data: dict,
     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']
 
@@ -433,8 +426,8 @@ def IC_peak_energy_charge_combined_plot(config_data: dict,
     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)
@@ -486,32 +479,29 @@ def higher_multipole_ic_correction(kappaR, abar, higher_l: int = 4):
 
 
 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='pp', min_kappaR=0.01, max_kappaR=50,
-    #                     save_as=Path('/home/andraz/ChargedShells/Figures/Emanuele_data/peak_pp_kappaR_higher_correction_abar02.png'),
+    # IC_peak_energy_plot(a_bar=a_bar, which='pp', min_kappaR=0.01, max_kappaR=50,
+    #                     save_as=FIGURES_PATH.joinpath('Emanuele_data').joinpath('peak_pp_kappaR_higher_correction_abar02.png'),
     #                     save_data=False,
     #                     quad_correction=False,
     #                     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'),
+    # IC_peak_energy_charge_plot(a_bar=a_bar, which='sep',
+    #                            # save_as=FIGURES_PATH.joinpath('Emanuele_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'),
+    # IC_peak_energy_abar_plot(kappaR=kappaR, which='sep', min_abar=0.2, max_abar=0.8,
+    #                     save_as=FIGURES_PATH.joinpath('Emanuele_data').joinpath('peak_sep_abar.png'),
     #                     save_data=False,
     #                     quad_correction=False,
     #                     log_y=False

analysis/peak_heigth_janus.py

@@ -1,13 +1,11 @@
 from matplotlib import gridspec
-from charged_shells import expansion, interactions, mapping, functions
+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
 import peak_heigth
 
 Array = np.ndarray
@@ -95,14 +93,13 @@ def get_charge_energy_dicts(peak: peak_heigth.Peak, params: ModelParams, emanuel
     return data_dict_ic, data_dict_cs
 
 
-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_5_JANUS")).joinpath("FIG_5_JANUS_NEW_CHARGE")
+    # 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']
 
@@ -110,8 +107,8 @@ def IC_peak_energy_charge_combined_plot(config_data: dict,
     sigma0 = np.linspace(-0.00099, 0.00099, 300)
     sigma_tilde = 0.00099
 
-    ex = expansion.MappedExpansionDipole(np.sort(np.array(a_bar)),  # sorting necessary as it is also in energy_dicts()
-                                       kappaR=3, sigma0=sigma0, l_max=20, sigma_tilde=sigma_tilde)
+    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)
@@ -164,19 +161,12 @@ def IC_peak_energy_charge_combined_plot(config_data: dict,
 
 
 def main():
-    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
-        config_data = json.load(config_file)
 
     a_bar = [0.2, 0.5, 0.8]
-    # IC_peak_energy_kappaR_combined_plot(config_data, a_bar,
-    #                                     save_as=Path(
-    #                                         '/home/andraz/ChargedShells/Figures/final_figures/peak_combined_kappaR.png')
-    #                                     )
-
     # a_bar = [0.1, 0.2, 0.3]
-    # IC_peak_energy_charge_combined_plot(config_data, a_bar,
-    #                                     save_as=Path('/home/andraz/ChargedShells/Figures/final_figures/janus_peak_combined_charge.png')
-    #                                     )
+    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,13 +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
@@ -31,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
@@ -58,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']
 
@@ -87,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)
@@ -100,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)
@@ -113,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)
@@ -126,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 = []
@@ -148,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']
@@ -170,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']
@@ -190,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']
@@ -209,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.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 = []
@@ -259,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():
@@ -277,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 = []
@@ -301,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 = []
@@ -314,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)
@@ -333,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:
@@ -345,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)
@@ -364,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)
@@ -392,7 +388,7 @@ 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.00099
     kappaR_list = [1, 3, 10]
@@ -401,8 +397,7 @@ def combined_all(config_data: dict, save_as=None):
     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 = []
@@ -411,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)
@@ -419,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:
@@ -429,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)
@@ -438,14 +432,14 @@ 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)
@@ -491,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,
@@ -511,39 +502,39 @@ 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("Emanuele_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("Emanuele_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("Emanuele_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')
                  )
 
 

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("Emanuele_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_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)
+    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):
@@ -71,8 +69,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 +96,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 +113,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 +128,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)

config.json

@@ -1,6 +0,0 @@
-{
-  "emanuele_data":  "/home/andraz/ChargedShells/data_Emanuele/2024-08-26",
-  "figures": "/home/andraz/ChargedShells/Figures",
-  "figure_data": "/home/andraz/ChargedShells/Figures/Data",
-  "expansion_data": "/home/andraz/ChargedShells/ExpansionData/"
-}

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)