2024-01-16

analysis/energy_gap.py

@@ -0,0 +1,266 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from charged_shells import expansion, interactions, mapping
+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
+
+Expansion = expansion.Expansion
+
+
+def energy_gap(ex1: Expansion, params: ModelParams, dist=2., match_expansion_axis_to_params=None):
+    ex2 = ex1.clone()
+    ex3 = ex1.clone()
+
+    ex2.rotate_euler(alpha=0, beta=np.pi/2, gamma=0)  # to get EP config between ex1 and ex2
+
+    energy_fn = mapping.parameter_map_two_expansions(partial(interactions.charged_shell_energy, dist=dist),
+                                                     match_expansion_axis_to_params)
+    energy_ep = energy_fn(ex1, ex2, params)
+    energy_pp = energy_fn(ex1, ex3, params)
+
+    return (energy_pp - energy_ep) / energy_pp
+
+
+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 = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
+    gap = energy_gap(ex, params, match_expansion_axis_to_params=1)
+
+    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
+    for g, lbl in zip(gap, [rf'$\bar a={a}$' for a in a_bar]):
+        ax.plot(kappaR, g, label=lbl)
+    ax.legend(fontsize=17)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\kappa R$', fontsize=15)
+    ax.set_ylabel(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', fontsize=20)
+    plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+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 = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
+    gap = energy_gap(ex, params, match_expansion_axis_to_params=1)
+
+    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
+    for g, lbl in zip(gap.T, [rf'$\kappa R={kR}$' for kR in kappaR]):
+        ax.plot(a_bar, g, label=lbl)
+    ax.legend(fontsize=17)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\bar a$', fontsize=15)
+    ax.set_ylabel(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', fontsize=20)
+    plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+def charge_kappaR_dependence(a_bar, min_charge, max_charge, save_as=None, cmap=cm.jet):
+
+    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 = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
+    gap = energy_gap(ex, params, match_expansion_axis_to_params=0)
+
+    colors = cmap(np.linspace(0, 1, len(charge)))
+    sm = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=np.min(charge), vmax=np.max(charge)))
+    sm.set_array([])
+
+    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
+    for g, lbl, c in zip(gap.T, [rf'$\sigma_0 ={1000 * c:.2f}$' for c in charge], colors):
+        ax.plot(kappaR, g, label=lbl, c=c)
+    # ax.legend(fontsize=17)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\kappa R$', fontsize=15)
+    ax.set_ylabel(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', fontsize=20)
+
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes('right', size='5%', pad=0.05)
+    cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
+    cbar.set_label(r'$\tilde \sigma_0$', rotation=90, labelpad=15, fontsize=12)
+
+    plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+def charge_kappaR_dependence_heatmap(a_bar, min_charge, max_charge, save_as=None, cmap=cm.jet):
+
+    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 = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
+    gap = energy_gap(ex, params, match_expansion_axis_to_params=0)
+
+    norm = TwoSlopeNorm(vmin=np.min(gap), vcenter=0, vmax=np.max(gap))
+    sm = cm.ScalarMappable(cmap=cmap, norm=norm)
+    sm.set_array([])
+
+    def y_formatter(x, pos):
+        return f"{charge[int(x)-1]:.2f}"
+
+    def x_formatter(x, pos):
+        return f"{kappaR[int(x)-1]:.2f}"
+
+    fig, ax = plt.subplots(figsize=plt.figaspect(1))
+    ax.imshow(gap.T, cmap=cmap, origin='lower',
+              # extent=[kappaR.min(), kappaR.max(), charge.min(), charge.max()]
+              )
+    # ax.legend(fontsize=17)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\kappa R$', fontsize=15)
+    ax.set_ylabel(r'$\tilde \sigma_0$', fontsize=15)
+
+    plt.gca().xaxis.set_major_formatter(FuncFormatter(x_formatter))
+    plt.gca().yaxis.set_major_formatter(FuncFormatter(y_formatter))
+
+    # ax.set_xticks(kappaR)
+    # ax.set_yticks(charge)
+
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes('right', size='5%', pad=0.05)
+    cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
+    cbar.set_label(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', rotation=90, labelpad=20, fontsize=12)
+
+    plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+def IC_gap_plot(config_data: dict, save_as=None):
+    em_data_path = (Path(config_data["emanuele_data"]).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]
+        print(k, data.shape)
+        for i in range(3):
+            print(np.unique(data[:, i]))
+        print('\n')
+
+
+def IC_gap_kappaR(config_data: dict, save_as=None):
+    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
+    em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
+    data = em_data['fixA']
+
+    print(data)
+
+    sort = np.argsort(data[:, 2])
+    xdata = data[:, 2][sort]
+    ydata = data[:, 3][sort]
+
+    plt.plot(xdata, ydata)
+    plt.xlabel('kappaR')
+    plt.ylabel('gap')
+    plt.show()
+
+
+def IC_gap_abar(config_data: dict, save_as=None):
+    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
+    em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
+    data = em_data['fixM']
+
+    print(data)
+
+    sort = np.argsort(data[:, 1])
+    xdata = data[:, 1][sort]
+    ydata = data[:, 3][sort]
+
+    plt.plot(xdata, ydata)
+    plt.xlabel('abar')
+    plt.ylabel('gap')
+    plt.show()
+
+
+def IC_gap_charge_at_abar(a_bar, config_data: dict, save_as=None, cmap=cm.coolwarm, which_change='changezp'):
+    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
+    em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
+    data = em_data[which_change]
+
+    relevant_indices = data[:, 1] == a_bar
+    if not np.any(relevant_indices):
+        raise ValueError(f'No results for given a_bar = {a_bar}. Possible values: {np.unique(data[:, 1])}')
+
+    data = data[relevant_indices]
+    print('Unique gap sizes', np.unique(data[:, 3]))
+
+    charge, inverse, counts = np.unique(data[:, 0], return_counts=True, return_inverse=True)
+
+    colors = cmap(np.linspace(0, 1, len(charge)))
+    sm = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=np.min(charge), vmax=np.max(charge)))
+    sm.set_array([])
+
+    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
+    for i, c in enumerate(colors):
+        idx, = np.nonzero(inverse == i)
+        kR = data[idx, 2]
+        gap = data[idx, 3]
+        sort = np.argsort(kR)
+        ax.plot(kR[sort], gap[sort], c=c)
+    # ax.legend(fontsize=17)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\kappa R$', fontsize=15)
+    ax.set_ylabel(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', fontsize=20)
+    ax.set_xlim(-0.25, 10.25)
+    # ax.set_ylim(-4.25, 6.25)
+
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes('right', size='5%', pad=0.05)
+    cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
+    cbar.set_label(r'$\tilde \sigma_0$', rotation=90, labelpad=15, fontsize=12)
+
+    plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+def main():
+    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
+        config_data = json.load(config_file)
+
+    # abar_kappaR_dependence(Path("/home/andraz/ChargedShells/Figures/full_amplitude_kappaR_dep.png"))
+    # abar_kappaR_dependence2(Path("/home/andraz/ChargedShells/Figures/full_amplitude_abar_dep.png"))
+
+    # charge_kappaR_dependence(a_bar=0.5, min_charge=-0.003, max_charge=0.003,
+    #                          save_as=Path("/home/andraz/ChargedShells/Figures/full_amplitude_charge_abar05.png"),
+    #                          cmap=cm.coolwarm)
+
+    # charge_kappaR_dependence_heatmap(a_bar=0.5, min_charge=-0.003, max_charge=0.003,
+    #                          save_as=Path("/home/andraz/ChargedShells/Figures/full_amplitude_heatmap_abar05.png"),
+    #                          cmap=cm.bwr)
+
+    # IC_gap_plot(config_data)
+
+    # IC_gap_kappaR(config_data)
+    # IC_gap_abar(config_data)
+
+    IC_gap_charge_at_abar(0.8, config_data, which_change='changezc',
+                          # save_as=Path("/home/andraz/ChargedShells/Figures/Emanuele_data/IC_full_amplitude_charge_abar08.png")
+                          )
+
+
+if __name__ == '__main__':
+
+    main()

analysis/patch_shape_comparison.py

@@ -4,6 +4,8 @@ from charged_shells.parameters import ModelParams
 import matplotlib.pyplot as plt
 import scipy.special as sps
 from pathlib import Path
+from functools import partial
+from matplotlib.lines import Line2D
 
 
 def point_to_gauss_map(sigma_m, a_bar, lbd, params: ModelParams):
@@ -16,38 +18,79 @@ def point_to_cap_map(sigma_m, a_bar, theta0, params: ModelParams):
             a_bar ** 2 / (sps.eval_legendre(1, np.cos(theta0)) - sps.eval_legendre(3, np.cos(theta0))))
 
 
-if __name__ == '__main__':
+def ic_cs_comparison():
+    kappaR = 15
+    params = ModelParams(R=150, kappaR=kappaR)
+    sigma_m = 0.001
+    sigma0 = 0.0001
+    a_bar = 0.3
+
+    theta = np.linspace(0, np.pi, 1001)
+    phi = 0.
+    dist = 1
+
+    ex_cp = expansion.MappedExpansionQuad(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,
+                                                                (sigma_m, sigma_m), params, 30)
+
+    fig, ax = plt.subplots()
+    ax.scatter(theta[::50], 1000 * potential_ic.T[::50], marker='o', label='ICi', facecolors='none',
+               edgecolors='tab:red')
+    ax.plot(theta, 1000 * potential_cp.T, label='CSp - mapped', linewidth=2)
+    # ax.plot(1000 * potential_ic.T - 1000 * potential_cp.T)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\theta$', fontsize=15)
+    ax.set_ylabel(r'$\phi$ [mV]', fontsize=15)
+    custom_ticks = [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4, np.pi]
+    custom_labels = ['0', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$', r'$\pi$']
+    plt.axhline(y=0, color='black', linestyle=':')
+    plt.xticks(custom_ticks, custom_labels, fontsize=13)
+    plt.legend(fontsize=12)
+    plt.tight_layout()
+    # plt.savefig(Path("/home/andraz/ChargedShells/Figures/potential_shape_comparison_overcharged05.pdf"), dpi=300)
+    plt.show()
+
+
+def models_comparison():
     target_patch_size = 0.92
-    params = ModelParams(R=150, kappaR=3)
+    kappaR = 3
+    params = ModelParams(R=150, kappaR=kappaR)
     sigma_m = 0.001
+    sigma0 = 0  # taking this total charge parameter nonzero causes cap model to fail in finding the appropriate theta0
+
+    sigma0_mapped = expansion.net_charge_map(sigma0, kappaR)
 
     def fn1(x):
-        return expansion.MappedExpansionQuad(a_bar=x, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30)
+        return expansion.MappedExpansionQuad(a_bar=x, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30, sigma0=sigma0)
 
     def fn2(x):
-        return expansion.GaussianCharges(lambda_k=x, omega_k=np.array([[0, 0], [np.pi, 0]]), sigma1=0.001, l_max=30)
+        return expansion.GaussianCharges(lambda_k=x, omega_k=np.array([[0, 0], [np.pi, 0]]), sigma1=0.001, l_max=30, sigma0=sigma0_mapped)
 
     def fn3(x):
-        return expansion.SphericalCap(theta0_k=x, sigma1=0.001, l_max=50, omega_k=np.array([[0, 0], [np.pi, 0]]))
+        return expansion.SphericalCap(theta0_k=x, sigma1=0.001, l_max=50, omega_k=np.array([[0, 0], [np.pi, 0]]), sigma0=sigma0_mapped)
 
     a_bar = patch_size.inverse_potential_patch_size(target_patch_size, fn1, 0.5, params)
     lbd = patch_size.inverse_potential_patch_size(target_patch_size, fn2, 5, params)
     theta0 = patch_size.inverse_potential_patch_size(target_patch_size, fn3, 0.5, params)
 
-    ex_point = expansion.MappedExpansionQuad(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30)
+    ex_point = expansion.MappedExpansionQuad(a_bar=a_bar, kappaR=params.kappaR, sigma_tilde=sigma_m, l_max=30, sigma0=sigma0)
 
     gauss_sigma = point_to_gauss_map(sigma_m, a_bar, lbd, params)
-    ex_gauss = expansion.GaussianCharges(lambda_k=lbd, omega_k=np.array([[0, 0], [np.pi, 0]]), sigma1=gauss_sigma, l_max=30)
+    ex_gauss = expansion.GaussianCharges(lambda_k=lbd, omega_k=np.array([[0, 0], [np.pi, 0]]), sigma1=gauss_sigma,
+                                         l_max=30, sigma0=sigma0_mapped)
 
     cap_sigma = point_to_cap_map(sigma_m, a_bar, theta0, params)
-    ex_cap = expansion.SphericalCap(theta0_k=theta0, sigma1=cap_sigma, omega_k=np.array([[0, 0], [np.pi, 0]]), l_max=30)
+    ex_cap = expansion.SphericalCap(theta0_k=theta0, sigma1=cap_sigma, omega_k=np.array([[0, 0], [np.pi, 0]]), l_max=30, sigma0=sigma0_mapped)
 
     theta = np.linspace(0, np.pi, 1001)
     phi = 0.
     dist = 1
 
-    potential_ic = potentials.inverse_patchy_particle_potential(theta, dist, a_bar, -2 * sigma_m, (sigma_m, sigma_m), params, 30)
+    potential_ic = potentials.inverse_patchy_particle_potential(theta, dist, a_bar, -2 * sigma_m + sigma0,
+                                                                (sigma_m, sigma_m), params, 30)
     potential1 = potentials.charged_shell_potential(theta, phi, dist, ex_point, params)
     potential2 = potentials.charged_shell_potential(theta, phi, dist, ex_gauss, params)
     potential3 = potentials.charged_shell_potential(theta, phi, dist, ex_cap, params)
@@ -57,7 +100,8 @@ if __name__ == '__main__':
     # expansion.plot_theta_profile_multiple([ex_point, ex_gauss, ex_cap], ['IC', 'Gauss', 'cap'], num=1000)
 
     fig, ax = plt.subplots()
-    ax.scatter(theta[::50], 1000 * potential_ic.T[::50], marker='o', label='ICi', facecolors='none', edgecolors='tab:red')
+    ax.scatter(theta[::50], 1000 * potential_ic.T[::50], marker='o', label='ICi', facecolors='none',
+               edgecolors='tab:red')
     ax.plot(theta, 1000 * potential1.T, label='CSp - mapped', linewidth=2)
     # ax.plot(theta, potential_ic.T, label='IC', ls=':', linewidth=2, marker='o', markevery=50, mfc='none')
     ax.plot(theta, 1000 * potential2.T, label='CSp - Gauss', linewidth=2, ls='--')
@@ -75,10 +119,138 @@ if __name__ == '__main__':
     plt.savefig(Path("/home/andraz/ChargedShells/Figures/potential_shape_comparison.pdf"), dpi=300)
     plt.show()
 
-    # path_comparison = rotational_path.PathExpansionComparison([ex_point, ex_gauss, ex_cap],
-    #                                                           path_plot.QuadPath,
-    #                                                           dist=2,
-    #                                                           params=params)
-    # path_comparison.plot(['IC', 'Gauss', 'cap'],
-    #                      save_as=Path("/home/andraz/ChargedShells/Figures/energy_shape_comparison_kR1.png"))
+
+def abar_comparison():
+    target_patch_sizes = [0.8, 0.85, 0.92]
+    params = ModelParams(R=150, kappaR=3)
+    sigma_m = 0.001
+
+    theta = np.linspace(0, np.pi, 1001)
+    phi = 0.
+    dist = 1
+
+    def fn1(x):
+        return expansion.MappedExpansionQuad(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)
+        pots.append(potentials.charged_shell_potential(theta, phi, dist, ex_point, params))
+
+    fig, ax = plt.subplots()
+    for pot, ps in zip(pots, target_patch_sizes):
+        ax.plot(theta, 1000 * pot, label=fr'$\theta_0={180 / np.pi * ps:.1f}^o$', linewidth=2)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\theta$', fontsize=15)
+    ax.set_ylabel(r'$\phi$ [mV]', fontsize=15)
+    custom_ticks = [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4, np.pi]
+    custom_labels = ['0', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$', r'$\pi$']
+    plt.axhline(y=0, color='black', linestyle=':')
+    plt.xticks(custom_ticks, custom_labels, fontsize=13)
+    plt.legend(fontsize=12)
+    plt.tight_layout()
+    plt.savefig(Path("/home/andraz/ChargedShells/Figures/potential_amplitude_comparison.pdf"), dpi=300)
+    plt.show()
+
+
+def charge_comparsion():
+
+    charges = np.array([-0.001, 0, 0.001, 0.002])
+    a_bar = 0.6
+    params = ModelParams(R=150, kappaR=10)
+    sigma_m = 0.001
+
+    theta = np.linspace(0, np.pi, 1001)
+    phi = 0.
+    dist = 1
+
+    ex_point = expansion.MappedExpansionQuad(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()
+    for p, lbl in zip(pot, [fr'$\sigma_0={c}$' for c in charges]):
+        ax.plot(theta, 1000 * p, linewidth=2, label=lbl)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\theta$', fontsize=15)
+    ax.set_ylabel(r'$\phi$ [mV]', fontsize=15)
+    custom_ticks = [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4, np.pi]
+    custom_labels = ['0', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$', r'$\pi$']
+    plt.axhline(y=0, color='black', linestyle=':')
+    plt.xticks(custom_ticks, custom_labels, fontsize=13)
+    plt.legend(fontsize=12)
+    plt.tight_layout()
+    # plt.savefig(Path("/home/andraz/ChargedShells/Figures/potential_charge_comparison.pdf"), dpi=300)
+    plt.show()
+
+
+# TODO: comparison of patch shape at different kappaR and a neutral particle, with fixed patch size
+
+
+def IC_CS_comparison():
+    target_patch_size = 0.92
+    kappaR = 3
+    params = ModelParams(R=150, kappaR=kappaR)
+    sigma_m = 0.001
+    sigma0_array = np.array([-0.0002, 0, 0.0002])
+
+    theta = np.linspace(0, np.pi, 1001)
+    phi = 0.
+    dist = 1
+
+    def fn1(x):
+        return expansion.MappedExpansionQuad(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)
+
+    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)
+
+        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))
+
+    fig, ax = plt.subplots()
+    lines = []
+    for p_ic, p_cs, charge in zip(potential_ic, potential_cs, sigma0_array):
+        ax.scatter(theta[::50], 1000 * p_ic.T[::50], marker='o', facecolors='none',
+                   edgecolors='tab:red')
+        l, = ax.plot(theta, 1000 * p_cs.T, label=rf'$\tilde \sigma_0 = {charge}$', linewidth=2)
+        lines.append(l)
+    # ax.plot(theta, potential_ic.T, label='IC', ls=':', linewidth=2, marker='o', markevery=50, mfc='none')
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\theta$', fontsize=15)
+    ax.set_ylabel(r'$\phi$ [mV]', fontsize=15)
+    custom_ticks = [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4, np.pi]
+    custom_labels = ['0', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$', r'$\pi$']
+    plt.axhline(y=0, color='black', linestyle=':')
+    plt.axvline(x=target_patch_size, color='black', linestyle=':')
+    plt.xticks(custom_ticks, custom_labels, fontsize=13)
+
+    lines.append(Line2D([0], [0], color='k', label='CSp'))
+    lines.append(plt.scatter([], [], marker='o', facecolor='none', edgecolor='red', label='ICi'))
+    plt.legend(handles=lines, fontsize=12)
+
+    plt.tight_layout()
+    plt.savefig(Path("/home/andraz/ChargedShells/Figures/potential_ic_cs_comparison_new.pdf"), dpi=300)
+    plt.show()
+
+
+def main():
+    # models_comparison()
+    # ic_cs_comparison()
+    IC_CS_comparison()
+
+    # abar_comparison()
+    # charge_comparsion()
+
+
+if __name__ == '__main__':
+
+    main()
 

analysis/patch_size_dependence.py

@@ -3,6 +3,7 @@ from pathlib import Path
 import numpy as np
 from charged_shells import expansion, patch_size
 from charged_shells.parameters import ModelParams
+import time
 
 
 def plot_abar_dependence(*, save=False):
@@ -15,7 +16,8 @@ def plot_abar_dependence(*, save=False):
 
     markers = ['o', 's', '^', 'D', 'v', '<', '>', 'p', '*', 'H', '+', 'x']
 
-    fig, ax = plt.subplots(figsize=plt.figaspect(0.5)/1.5)
+    # fig, ax = plt.subplots(figsize=plt.figaspect(0.5)/1.5)
+    fig, ax = plt.subplots()
     for patch, kR, marker in zip(ps.T, kappaR, markers):
         ax.plot(a_bar, patch, label=rf'$\kappa R$ = {kR}', marker=marker, markerfacecolor='none', markevery=5)
     ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
@@ -24,7 +26,7 @@ def plot_abar_dependence(*, save=False):
     plt.legend(fontsize=14)
     plt.tight_layout()
     if save:
-        plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential.pdf"), dpi=300)
+        plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential_new.pdf"), dpi=300)
     plt.show()
 
 
@@ -77,47 +79,53 @@ def plot_kappaR_charge_dependence(*, normalized=False, save=False):
 def plot_sigma0_dependence(*, save=False):
 
     # kappaR = np.array([0.1, 1, 3, 10, 30])
-    kappaR = np.linspace(0.1, 15, 101)
+    kappaR = np.linspace(0.1, 15, 201)
     # sigma0 = np.linspace(-0.001, 0.0015, 6)
-    sigma0 = np.array([-0.001, 0, 0.001])
+    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 = expansion.MappedExpansionQuad(a_bar=0.3, sigma_tilde=0.001, l_max=30, kappaR=kappaR, sigma0=sigma0)
 
-    ps = patch_size.potential_patch_size(ex, params, match_expansion_axis_to_params=0)
+    t0 = time.perf_counter()
+    ps = patch_size.potential_patch_size(ex, params, match_expansion_axis_to_params=0, skip_nozero_cases=True)
+    t1 = time.perf_counter()
+    print(f'Time: {t1 - t0}')
 
     markers = ['o', 's', '^', 'D', 'p', 'v', '<', '>', 'x', '*', 'H', '+']
 
     fig, ax = plt.subplots(figsize=plt.figaspect(0.5)/1.5)
     for patch, sgm, marker in zip(ps.T, sigma0, markers):
-        ax.plot(kappaR, patch, label=rf'$\tilde \sigma_0$ = {sgm}', marker=marker, markerfacecolor='none', markevery=5)
+        nonzero_ps = np.nonzero(patch)
+        ax.plot(kappaR[nonzero_ps], patch[nonzero_ps], label=rf'$\tilde \sigma_0$ = {sgm}', marker=marker, markerfacecolor='none', markevery=10)
     ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
     ax.set_xlabel(r'$\kappa R$', fontsize=15)
     ax.set_ylabel(r'$\theta_0$', fontsize=15)
-    plt.legend(fontsize=14)
+    plt.legend(fontsize=14, loc='upper right')
     plt.tight_layout()
     if save:
-        plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential_charge_abar08.pdf"), dpi=300)
+        plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential_charge_abar03.pdf"), dpi=300)
     plt.show()
 
 
 def plot_sigma0_dependence_relative(*, save=False):
 
     # kappaR = np.array([0.1, 1, 3, 10, 30])
-    kappaR = np.linspace(0.1, 15, 101)
+    kappaR = np.linspace(0.1, 15, 201)
     # sigma0 = np.linspace(-0.001, 0.0015, 6)
-    sigma0 = np.array([-0.001, 0, 0.001])
+    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 = 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)
 
-    ps = patch_size.potential_patch_size(ex, params, match_expansion_axis_to_params=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)
 
     markers = ['o', 's', '^', 'D', 'p', 'v', '<', '>', 'x', '*', 'H', '+']
 
     fig, ax = plt.subplots(figsize=plt.figaspect(0.5)/1.5)
     for patch, sgm, marker in zip(ps.T, sigma0, markers):
-        ax.plot(kappaR, patch / ps_neutral, label=rf'$\tilde \sigma_0$ = {sgm}', marker=marker, markerfacecolor='none', markevery=5)
+        nonzero_ps = np.nonzero(patch)
+        ax.plot(kappaR[nonzero_ps], patch[nonzero_ps] / ps_neutral[nonzero_ps],
+                label=rf'$\tilde \sigma_0$ = {sgm}', marker=marker, markerfacecolor='none', markevery=10)
     ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
     ax.set_xlabel(r'$\kappa R$', fontsize=15)
     ax.set_ylabel(r'$\theta_0$', fontsize=15)
@@ -128,11 +136,18 @@ def plot_sigma0_dependence_relative(*, save=False):
     plt.show()
 
 
-if __name__ == '__main__':
+def main():
+
+    plot_abar_dependence(save=True)
 
-    # plot_abar_dependence(save=False)
     # plot_sigma0_dependence(save=True)
+
     # plot_sigma0_dependence_relative(save=True)
 
     # plot_abar_charge_dependence(save=True)
-    plot_kappaR_charge_dependence(normalized=True, save=True)
+    # plot_kappaR_charge_dependence(normalized=True, save=True)
+
+
+if __name__ == '__main__':
+
+    main()

analysis/path_plot.py

@@ -1,7 +1,10 @@
 import numpy as np
 from charged_shells.rotational_path import PairRotationalPath, PathEnergyPlot
-from charged_shells import expansion, patch_size
+from charged_shells import expansion, patch_size, interactions
 from charged_shells.parameters import ModelParams
+import matplotlib.pyplot as plt
+from pathlib import Path
+import json
 
 Array = np.ndarray
 
@@ -43,72 +46,196 @@ DipolePath2.add_euler(beta2=zero_to_pi_half, beta1=pi_half_to_pi[::-1], alpha1=n
 DipolePath2.add_euler(beta2=pi_half_to_pi, beta1=zero_to_pi_half[::-1], alpha1=np.pi)
 
 
-def kappaR_dependence(kappaR: Array, abar: float, sigma_tilde=0.001, save=False):
+def model_comparison(save_as=None):
+    params = ModelParams(R=150, kappaR=3)
+    abar = 0.5
+    sigma_tilde = 0.001
+
+    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex2 = ex1.clone()
+
+    path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params)
+    energy = path_plot.evaluate_path()
+    x_axis = path_plot.rot_path.stack_x_axes()
+
+    # peak_energy_sanity_check
+    ex1new = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex2new = ex1new.clone()
+    pp_energy = interactions.charged_shell_energy(ex1new, ex2new, params)
+    print(f'PP energy: {pp_energy}')
+
+    # Emanuele data
+    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
+        config_data = json.load(config_file)
+    em_data = np.load(Path(config_data["emanuele_data"]).joinpath("FIG_7").joinpath("pathway.npz"))['arr_0']
+
+    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
+    ax.plot(em_data[:, 0], em_data[:, 1], label='ICi')
+    ax.plot(x_axis, np.squeeze(energy), label='CSp')
+    # ax.plot(x_axis, em_data[:, 1] / np.squeeze(energy), label='CSp')
+    PathEnergyPlot.plot_style(fig, ax)
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+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, 0.001, l_max=30)
+    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
     ex2 = ex1.clone()
 
-    # charge_profile = ex1.charge_value(theta=np.linspace(0, np.pi, 100), phi=0)
-    # plt.plot(charge_profile.T)
-    # plt.show()
+    path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=0)
 
-    # expansion.plot_theta_profile(ex1, num=1000, theta_end=np.pi, phi=0)
-    #
-    # ps = patch_size.potential_patch_size(ex1, params, theta1=np.pi / 2,
-    #                                      # match_expansion_axis_to_params=1
-    #                                      )
-    # print(ps)
+    path_plot.plot(labels=[rf'$\kappa R$={kR}' for kR in kappaR],
+                   # norm_euler_angles={'beta2': np.pi},
+                   save_as=save_as)
 
-    path_plot = PathEnergyPlot(ex1, ex2, DipolePath2, dist=2., params=params, match_expansion_axis_to_params=None)
 
-    path_plot.plot(  # labels=[rf'$\theta$={180 * patch / np.pi:.1f}' for patch in np.array([0.5])],
-        # norm_euler_angles={'beta2': np.pi},
-        # save_as=Path("/home/andraz/ChargedShells/Figures/dipole_path_kappaR3.png")
-    )
+def abar_dependence(abar: Array, kappaR: float, sigma_tilde=0.001, save_as=None):
+    params = ModelParams(R=150, kappaR=kappaR)
 
-    # path_plot.plot(labels=[rf'$\kappa R$={kR}' for kR in kappaR],
-    #                # norm_euler_angles={'beta2': np.pi},
-    #                # save_as=Path("/home/andraz/ChargedShells/Figures/dipole_path_lambda5_norm2.png")
-    #                )
+    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30)
+    ex2 = ex1.clone()
 
-    # path_plot.plot_sections(save_as=Path('/home/andraz/ChargedShells/Figures/dipole_path2.png'))
+    path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=None)
 
+    path_plot.plot(labels=[rf'$\bar a$={a}' for a in abar],
+                   # norm_euler_angles={'beta2': np.pi},
+                   save_as=save_as)
 
 
-if __name__ == '__main__':
+def sigma0_dependence(sigma0: Array, kappaR: float, abar: float, sigma_tilde=0.001, save_as=None):
+    params = ModelParams(R=150, kappaR=kappaR)
 
-    # kappaR = np.array([1, 3, 10])
-    params = ModelParams(R=150, kappaR=3)
-    # ex1 = expansion.MappedExpansionQuad(0.4, params.kappaR, 0.001, l_max=20)
-    # ex1 = expansion.Expansion24(sigma2=0.001, sigma4=0)
-    # ex1 = expansion.Expansion(l_array=np.array([1]), coefs=expansion.rot_sym_expansion(np.array([1]), np.array([0.001])))
-    # ex1 = expansion.GaussianCharges(omega_k=np.array([[0, 0], [np.pi, 0]]), lambda_k=np.array([5, 10]), sigma1=0.001, l_max=10)
-    ex1 = expansion.GaussianCharges(omega_k=np.array([0, 0]), lambda_k=np.array([1, 5, 10]), sigma1=0.001, l_max=20)
-    # ex1 = expansion.SphericalCap(np.array([[0, 0], [np.pi, 0]]), np.array([0.5]), 0.001, 20)
+    ex1 = expansion.MappedExpansionQuad(abar, params.kappaR, sigma_tilde, l_max=30, sigma0=sigma0)
     ex2 = ex1.clone()
 
-    # charge_profile = ex1.charge_value(theta=np.linspace(0, np.pi, 100), phi=0)
-    # plt.plot(charge_profile.T)
-    # plt.show()
+    path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params, match_expansion_axis_to_params=None)
+
+    path_plot.plot(labels=[rf'$\tilde \sigma_0={s0}$' for s0 in sigma0],
+                   # norm_euler_angles={'beta2': np.pi},
+                   save_as=save_as)
+
+
+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)
+    ex2 = ex1.clone()
+
+    plots = []
+    for d in dist:
+        path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=d, params=params)
+        x = d * kappaR
+        plots.append(path_plot.evaluate_path() * np.exp(x) * x)
+
+    x_axis = path_plot.rot_path.stack_x_axes()
+    labels = [rf'$\rho/R ={d}$' for d in dist]
 
-    # expansion.plot_theta_profile(ex1, num=1000, theta_end=np.pi, phi=0)
+    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
+    for pl, lbl in zip(plots, labels):
+        ax.plot(x_axis, pl, label=lbl)
+    PathEnergyPlot.plot_style(fig, ax)
+    ax.set_ylabel(r'$U \kappa\rho e^{\kappa\rho}$', fontsize=15)
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+def IC_kappaR_dependence(config_data: dict, save_as=None):
+
+    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_8").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']
+    kR10 = np.load(em_data_path.joinpath("EMME_10.").joinpath("pathway.npz"))['arr_0']
+
+    labels = [rf'$\kappa R$={kR}' for kR in [1, 3, 10]]
+
+    fig, ax = plt.subplots()
+    ax.plot(kR1[:, 0], kR1[:, 1], label=labels[0])
+    ax.plot(kR3[:, 0], kR3[:, 1], label=labels[1])
+    ax.plot(kR10[:, 0], kR10[:, 1], label=labels[2])
+    PathEnergyPlot.plot_style(fig, ax)
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+def IC_abar_dependence(config_data: dict, save_as=None):
+
+    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_8").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']
+
+    labels =[rf'$\bar a$={a}' for a in [0.3, 0.4, 0.5]]
+
+    fig, ax = plt.subplots()
+    ax.plot(a03[:, 0], a03[:, 1], label=labels[0])
+    ax.plot(a04[:, 0], a04[:, 1], label=labels[1])
+    ax.plot(a05[:, 0], a05[:, 1], label=labels[2])
+    PathEnergyPlot.plot_style(fig, ax)
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    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"))
+    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']
+
+    labels = [rf'$\tilde \sigma_0={s0}$' for s0 in [-0.001, 0, 0.001]]
+
+    fig, ax = plt.subplots()
+    ax.plot(overchargerd[:, 0], overchargerd[:, 1], label=labels[0])
+    ax.plot(neutral[:, 0], neutral[:, 1], label=labels[1])
+    ax.plot(undercharged[:, 0], undercharged[:, 1], label=labels[2])
+    PathEnergyPlot.plot_style(fig, ax)
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+def main():
+
+    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
+        config_data = json.load(config_file)
+
+    # model_comparison(
+    #     # save_as=Path(config_data["figures"]).joinpath("Emanuele_data").joinpath('IC_CS_quadrupole_path.pdf')
+    # )
+
+    # kappaR_dependence(np.array([1, 3, 10]), 0.5,
+    #                   # save_as=Path("/home/andraz/ChargedShells/Figures/quadrupole_kappaR_dep.png")
+    #                   )
     #
-    # ps = patch_size.potential_patch_size(ex1, params, theta1=np.pi / 2,
-    #                                      # match_expansion_axis_to_params=1
-    #                                      )
-    # print(ps)
+    # abar_dependence(np.array([0.3, 0.4, 0.5]), 3,
+    #                 save_as=Path("/home/andraz/ChargedShells/Figures/quadrupole_abar_dep.png"))
 
-    path_plot = PathEnergyPlot(ex1, ex2, DipolePath2, dist=2., params=params, match_expansion_axis_to_params=None)
+    sigma0_dependence(np.array([0.001, 0.002, 0.003]), 3, 0.5,
+                      # save_as=Path("/home/andraz/ChargedShells/Figures/quadrupole_charge_dep_abar05_kappaR3.png")
+    )
 
-    path_plot.plot(# labels=[rf'$\theta$={180 * patch / np.pi:.1f}' for patch in np.array([0.5])],
-                   # norm_euler_angles={'beta2': np.pi},
-                   # save_as=Path("/home/andraz/ChargedShells/Figures/dipole_path_kappaR3.png")
-                   )
+    # distance_dependence(dist=np.array([2, 3, 4, 6, 10, 20]), kappaR=3, abar=0.5,
+    #                     # save_as=Path(config_data["figures"]).joinpath('quadrupole_distance_dep.png')
+    #                     )
+
+    # IC_kappaR_dependence(config_data,
+    #                      # save_as=Path(config_data["figures"]).joinpath("Emanuele_data"). joinpath('quadrupole_kappaR_dep.png')
+    #                      )
+    #
+    # IC_abar_dependence(config_data, save_as=Path(config_data["figures"]).joinpath("Emanuele_data").
+    #                    joinpath('quadrupole_abar_dep.png'))
+    #
+    # IC_sigma0_dependence(config_data, save_as=Path(config_data["figures"]).joinpath("Emanuele_data").
+    #                      joinpath('quadrupole_charge_dep_abar05_kappaR3.png'))
+
+
+if __name__ == '__main__':
 
-    # path_plot.plot(labels=[rf'$\kappa R$={kR}' for kR in kappaR],
-    #                # norm_euler_angles={'beta2': np.pi},
-    #                # save_as=Path("/home/andraz/ChargedShells/Figures/dipole_path_lambda5_norm2.png")
-    #                )
+    main()
 
-    # path_plot.plot_sections(save_as=Path('/home/andraz/ChargedShells/Figures/dipole_path2.png'))

analysis/peak_heigth.py

@@ -5,6 +5,7 @@ import numpy as np
 from typing import Literal
 from pathlib import Path
 from functools import partial
+import json
 
 Array = np.ndarray
 Expansion = expansion.Expansion
@@ -42,6 +43,7 @@ def peak_energy_plot(kappaR: Array,
     fig, ax = plt.subplots()
     for en, ab in zip(energy, a_bar):
         ax.plot(kappaR, en / en[0], label=rf'$\bar a = {ab:.1f}$')
+        # ax.plot(kappaR, en, label=rf'$\bar a = {ab:.1f}$')
     ax.legend(fontsize=12)
     ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
     ax.set_xlabel(r'$\kappa R$', fontsize=15)
@@ -52,11 +54,67 @@ def peak_energy_plot(kappaR: Array,
     plt.show()
 
 
+def IC_peak_energy_plot(config_data: dict,
+                        a_bar: list,
+        which: Literal['ep', 'pp'],
+                     save_as: Path = None):
+
+    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_11"))
+    em_data = np.load(em_data_path.joinpath("pair_energy_fig11.npz"))
+    data = em_data['fixA']
+
+    if which == 'ep':
+        column_idx = 4
+    elif which == 'pp':
+        column_idx = 3
+    else:
+        raise ValueError
+
+    abar, inverse, counts = np.unique(data[:, 1], return_counts=True, return_inverse=True)
+    # print(indices, counts)
+    print(inverse)
+    # sort_abar = np.argsort(indices)
+    # print(data[:, 1][indices])
+    # print(sort_abar)
+    # print(data[:, 1][indices[sort_abar]])
+
+    # energies = data[:, column_idx].reshape(-1, counts[0])
+    # kappaR = data[:, 2].reshape(-1, counts[0])
+    # print(len(kappaR), len(energies), len(abar))
+
+    fig, ax = plt.subplots()
+    for i in range(len(abar)):
+        if abar[i] in a_bar:
+            idx, = np.nonzero(inverse == i)
+            kR = data[idx, 2]
+            en = data[idx, column_idx]
+            sort = np.argsort(kR)
+            # ax.plot(kR[sort], en[sort] / en[sort][0], label=rf'$\bar a = {abar[i]:.2f}$')
+            ax.plot(kR[sort], en[sort], label=rf'$\bar a = {abar[i]:.2f}$')
+    ax.legend(fontsize=12)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+    ax.set_xlabel(r'$\kappa R$', fontsize=15)
+    ax.set_ylabel(rf'$\bar V_{{{which}}}$', fontsize=15)
+    plt.tight_layout()
+    ax.set_yscale('log')
+    ax.set_xscale('log')
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
 if __name__ == '__main__':
+    
+    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
+        config_data = json.load(config_file)
 
-    kappaR = np.arange(1, 10, 0.1)
+    kappaR = np.arange(0.5, 10, 0.1)
     a_bar = np.arange(0.2, 0.8, 0.2)
 
-    peak_energy_plot(kappaR, a_bar, which='ep',
-                     # save_as=Path('/home/andraz/ChargedShells/Figures/nonmonotonicity_check_ep.pdf')
+    # peak_energy_plot(kappaR, a_bar, which='pp',
+    #                  # save_as=Path('/home/andraz/ChargedShells/Figures/nonmonotonicity_check_ep.pdf')
+    #                  )
+    
+    IC_peak_energy_plot(config_data, a_bar=[0.2, 0.24, 0.36, 0.52, 0.8], which='pp',
+                     # save_as=Path('/home/andraz/ChargedShells/Figures/Emanuele_data/nonmonotonicity_check_ep.pdf')
                      )

analysis/sEE_minimum.py

@@ -1,17 +1,30 @@
+import matplotlib.pyplot as plt
 import numpy as np
 from charged_shells import expansion, interactions, mapping
 from charged_shells.parameters import ModelParams
 from functools import partial
+from pathlib import Path
+import json
+from matplotlib import cm
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+from labellines import labelLine, labelLines
+from matplotlib.colors import TwoSlopeNorm
+from matplotlib.ticker import FuncFormatter
 
 Expansion = expansion.Expansion
 
 
-def sEE_minimum(ex: Expansion, params: ModelParams, accuracy=1e-3, dist=2., match_expansion_axis_to_params=None):
+def sEE_minimum(ex: Expansion, params: ModelParams, accuracy=1e-2, dist=2., match_expansion_axis_to_params=None,
+                angle_start: float = 0, angle_stop: float = np.pi / 2):
     ex2 = ex.clone()
-    zero_to_pi_half = np.linspace(0, np.pi / 2, int(1 / accuracy), endpoint=True)
+    angle_range = np.linspace(angle_start, angle_stop, int((angle_stop - angle_start) / accuracy), endpoint=True)
+    # angle_range = np.linspace(0.1, 0.7, int(1 / accuracy), endpoint=True)
 
-    ex.rotate_euler(alpha=0, beta=zero_to_pi_half, gamma=0)
-    ex2.rotate_euler(alpha=0, beta=zero_to_pi_half, gamma=0)
+    ex.rotate_euler(alpha=0, beta=angle_range, gamma=0)
+    ex2.rotate_euler(alpha=0, beta=angle_range, gamma=0)
+
+    if match_expansion_axis_to_params is not None:
+        match_expansion_axis_to_params += 1
 
     energy_fn = mapping.parameter_map_two_expansions(partial(interactions.charged_shell_energy, dist=dist),
                                                      match_expansion_axis_to_params)
@@ -20,17 +33,191 @@ def sEE_minimum(ex: Expansion, params: ModelParams, accuracy=1e-3, dist=2., matc
     min_idx = np.argmin(energy, axis=0)
     min_energy = np.min(energy, axis=0)
 
-    return min_energy, zero_to_pi_half[min_idx]
+    return min_energy, angle_range[min_idx]
 
 
-def main():
+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)
+    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)
+
+    print(np.min(min_angle), np.max(min_angle))
+    print(np.min(min_energy), np.max(min_energy))
+
+    plt.contourf(kR_mesh, a_mesh, min_angle)
+    # plt.imshow(min_angle)
+    plt.show()
+
+    plt.contourf(kR_mesh, a_mesh, min_energy, np.array([-9, -5, -2.5, -1, -0.5, -0.2, 0]))
+    # plt.imshow(min_energy)
+    plt.show()
+
+
+def kappaR_dependence(kappaR, save_as=None, cmap=cm.jet):
 
-    kappaR = np.array([1, 3, 10])
+    a_bar = np.linspace(0.12, 0.8, 15)
     params = ModelParams(R=150, kappaR=kappaR)
-    ex = expansion.MappedExpansionQuad(0.4, kappaR, 0.001)
-    min_energy, min_angle = sEE_minimum(ex, params, match_expansion_axis_to_params=1)
-    print(min_energy, min_angle + np.pi/2)
 
+    ex = expansion.MappedExpansionQuad(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)
+    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.)
+
+    colors = cmap(np.linspace(0, 1, len(a_bar)))
+    sm = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=np.min(a_bar), vmax=np.max(a_bar)))
+    sm.set_array([])
+
+    # Transformation to the position along our quadrupolar rotational path
+    min_angle = np.pi - min_angle
+    min_angle_alt = np.pi - min_angle_alt
+
+    kappaR_labels = [rf'$\kappa R={kR:.1f}$' for kR in kappaR_alt]
+
+    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
+    for me, ma, lbl in zip(min_energy_alt.T, min_angle_alt.T, kappaR_labels):
+        ax.plot(ma, me, ls=':', c='k', label=lbl)
+    labelLines(ax.get_lines(), align=False, fontsize=14,  xvals=(2.35, 2.65))
+
+    for me, ma, lbl, c in zip(min_energy.T, min_angle.T, [rf'$\bar a={a:.2f}$' for a in a_bar], colors):
+        ax.plot(ma, me, label=lbl, c=c)
+    # ax.plot(min_angle, min_energy)
+    # ax.legend(fontsize=17)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel('angle', fontsize=15)
+    ax.set_ylabel('U', fontsize=20)
+    ax.set_xlim(2.2, 2.65)
+    ax.set_ylim(-20, 1)
+
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes('right', size='5%', pad=0.05)
+    cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
+    cbar.set_label(r'$\bar a$', rotation=90, labelpad=15, fontsize=15)
+
+    plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+def IC_kappaR_dependence(config_data: dict, 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 = np.load(em_data_path.joinpath(f"fig9{file_suffix}.npz"))['arr_0']
+    # print(em_data.shape)
+
+    a_bar, indices, counts = np.unique(em_data[:, 0], return_counts=True, return_inverse=True)
+    print(a_bar)
+
+    if not np.all(counts == counts[0]):
+        raise ValueError("Data not reshapable.")
+
+    colors = cmap(np.linspace(0, 1, len(a_bar)))
+    sm = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=np.min(a_bar), vmax=np.max(a_bar)))
+    sm.set_array([])
+
+    min_energy = em_data[:, 3].reshape(-1, counts[0])
+    min_angle = em_data[:, 2].reshape(-1, counts[0])
+
+    kappa_line_energy = []
+    kappa_line_angle = []
+    kappaR_labels = []
+    for kappa in which_kappa_lines:
+        ind = em_data[:, 1] == kappa
+        if np.sum(ind) == 0:
+            print(f'No lines with kR={kappa} in data. Available values: {np.unique(em_data[:, 1])}')
+            continue
+        kappa_line_energy.append(em_data[:, 3][ind])
+        kappa_line_angle.append(em_data[:, 2][ind])
+        kappaR_labels.append(rf'$\kappa R={kappa:.1f}$')
+
+
+    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
+    for me, ma, lbl in zip(kappa_line_energy, kappa_line_angle, kappaR_labels):
+        sort = np.argsort(ma)
+        ax.plot(ma[sort], me[sort], ls=':', c='k', label=lbl)
+
+    labelLines(ax.get_lines(), align=False, fontsize=14, xvals=(2.4, 2.65))
+
+    for me, ma, lbl, c in zip(min_energy, min_angle, [rf'$\bar a={a:.2f}$' for a in a_bar], colors):
+        ax.plot(ma, me, label=lbl, c=c)
+    # ax.plot(min_angle, min_energy)
+    # ax.legend(fontsize=17)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel('angle', fontsize=15)
+    ax.set_ylabel('U', fontsize=20)
+    ax.set_xlim(2.2, 2.65)
+    ax.set_ylim(-40, 1)
+
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes('right', size='5%', pad=0.05)
+    cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
+    cbar.set_label(r'$\bar a$', rotation=90, labelpad=15, fontsize=12)
+
+    plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+# def charge_dependence(charge, save_as=None):
+#
+#     a_bar = np.array([0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8])
+#     params = ModelParams(R=150, kappaR=kappaR)
+#
+#     ex = expansion.MappedExpansionQuad(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)
+#     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.)
+#
+#     fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
+#     for me, ma, lbl in zip(min_energy_alt.T, min_angle_alt.T, [rf'$\bar a={a:.2f}$' for a in a_bar]):
+#         ax.plot(ma, me, ls=':', c='k')
+#     for me, ma, lbl in zip(min_energy.T, min_angle.T, [rf'$\bar a={a:.2f}$' for a in a_bar]):
+#         ax.plot(ma, me, label=lbl)
+#     # ax.plot(min_angle, min_energy)
+#     ax.legend(fontsize=17)
+#     ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+#     ax.set_xlabel('angle', fontsize=15)
+#     ax.set_ylabel('U', fontsize=20)
+#     ax.set_xlim(0.56, 0.93)
+#     ax.set_ylim(-20, 1)
+#     plt.tight_layout()
+#     if save_as is not None:
+#         plt.savefig(save_as, dpi=600)
+#     plt.show()
+
+
+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('sEE_min_kappaR_abar.png')
+    #                      )
 
 
 if __name__ == '__main__':

charged_shells/expansion.py

@@ -1,7 +1,6 @@
 from __future__ import annotations
 import numpy as np
 from dataclasses import dataclass
-from functools import cached_property
 import charged_shells.functions as fn
 import quaternionic
 import spherical
@@ -55,7 +54,7 @@ class Expansion:
         self.coefs = self.coefs.reshape(shape + (self.coefs.shape[-1],))
         self._rotations = self._rotations.reshape(shape + (4,))
 
-    @cached_property
+    @property
     def lm_arrays(self) -> (Array, Array):
         """Return l and m arrays containing all (l, m) pairs."""
         return full_lm_arrays(self.l_array)
@@ -94,9 +93,13 @@ class Expansion:
 
 class Expansion24(Expansion):
 
-    def __init__(self, sigma2: float, sigma4: float, sigma0: float = 0.):
+    def __init__(self, sigma2: float | Array, sigma4: float | Array, sigma0: float | Array = 0.):
         l_array = np.array([0, 2, 4])
-        coefs = rot_sym_expansion(l_array, np.array([sigma0, sigma2, sigma4]))
+        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)
 
 
@@ -116,18 +119,22 @@ class MappedExpansionQuad(Expansion):
         :param l_max: maximal ell value for the expansion
         :param sigma0: total (mean) charge density
         """
-        a_bar, kappaR = np.broadcast_arrays(a_bar, kappaR)
+        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)
 
         l_array = np.array([l for l in range(l_max + 1) if l % 2 == 0])
-        a_bar, kappaR, l_array_expanded = np.broadcast_arrays(a_bar[..., None],
-                                                              kappaR[..., None],
-                                                              l_array[None, :])
-
-        coefs = (2 * sigma_tilde * fn.coef_C_diff(l_array_expanded, kappaR)
-                 * np.sqrt(4 * np.pi * (2 * l_array_expanded + 1)) * np.power(a_bar, l_array_expanded))
+        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))
-        coefs = expansion_total_charge(coefs, sigma0)
-        super().__init__(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))
+        super().__init__(l_array, np.squeeze(coefs))
 
 
 class GaussianCharges(Expansion):
@@ -211,6 +218,10 @@ class SphericalCap(Expansion):
         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)
+
+
 def full_lm_arrays(l_array: Array) -> (Array, Array):
     """From an array of l_values get arrays of ell and m that give you all pairs (ell, m)."""
     all_m_list = []
@@ -227,17 +238,18 @@ def rot_sym_expansion(l_array: Array, coefs: Array) -> Array:
     return full_coefs
 
 
-def expansion_total_charge(coefs: Array, sigma0: float | Array):
+def expansion_total_charge(coefs: Array, sigma0: float | Array | None):
     """Adds a new axis to the expansion coefficients that modifies expansion based on given net charge density."""
     if sigma0 is None:
         return coefs
     if not isinstance(sigma0, Array):
         x = copy.deepcopy(coefs)
-        x[..., 0] = sigma0 / np.sqrt(4 * np.pi)
+        x[..., 0] = sigma0 * np.sqrt(4 * np.pi)
         return x
-    sigma0 = sigma0.flatten()
-    x = np.repeat(np.expand_dims(coefs, -2), len(sigma0), axis=-2)
-    x[..., 0] = sigma0 / np.sqrt(4 * np.pi)
+
+    # 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
 
 

charged_shells/functions.py

@@ -7,7 +7,7 @@ def sph_bessel_i(n, x, **kwargs):
 
 
 def sph_bessel_k(n, x, **kwargs):
-    return 2 / np.pi * sps.spherical_kn(n, x, **kwargs)
+    return sps.spherical_kn(n, x, **kwargs)
 
 
 def sph_harm(l, m, theta, phi, **kwargs):

charged_shells/interactions.py

@@ -2,7 +2,6 @@ from charged_shells import expansion, parameters
 import charged_shells.functions as fn
 from py3nj import wigner3j
 import numpy as np
-import time
 from typing import Literal
 import charged_shells.units_and_constants as uc
 
@@ -16,7 +15,8 @@ EnergyUnit = Literal['kT', 'eV', 'J']
 def energy_units(units: EnergyUnit, params: ModelParams) -> float:
     match units:
         case 'eV':
-            return 1 / (uc.CONSTANTS.e0 * uc.UNITS.voltage)
+            # return 1 / (uc.CONSTANTS.e0 * uc.UNITS.voltage)
+            return 1.
         case 'kT':
             return 1 / (params.temperature * uc.CONSTANTS.Boltzmann)
         case 'J':
@@ -36,6 +36,13 @@ def charged_shell_energy(ex1: Expansion, ex2: Expansion, params: ModelParams, di
     flat_p = full_l_array[indices_p]
     flat_m = full_m_array[indices_l]  # the same as full_m_array[indices_p]
 
+    relevant_pairs, = np.nonzero(flat_l >= flat_p)
+    flat_l = flat_l[relevant_pairs]
+    flat_p = flat_p[relevant_pairs]
+    flat_m = flat_m[relevant_pairs]
+    indices_l = indices_l[relevant_pairs]
+    indices_p = indices_p[relevant_pairs]
+
     charge_factor = np.real(ex1.coefs[..., indices_l] * np.conj(ex2.coefs[..., indices_p]) +
                             (-1) ** (flat_l + flat_p) * ex1.coefs[..., indices_p] * np.conj(ex2.coefs[..., indices_l]))
 
@@ -75,4 +82,6 @@ def charged_shell_energy(ex1: Expansion, ex2: Expansion, params: ModelParams, di
     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)
 
-    return 0.5 * constants * np.sum(broadcasted_lspm_vals * charge_vals, axis=-1)
+    rescale_at_equal_lp = np.where(l_vals == p_vals, 0.5, 1)
+
+    return constants * np.sum(rescale_at_equal_lp * broadcasted_lspm_vals * charge_vals, axis=-1)

charged_shells/mapping.py

@@ -36,6 +36,17 @@ def unravel_params(params: ModelParams) -> list[ModelParams]:
         return [ModelParams(R=params.R, kappa=kappa) for kappa in params.kappa]
     if not (isinstance(params.R, Array) or isinstance(params.kappa, Array)):
         return [params]
+    raise NotImplementedError
+
+
+def unravel_expansion_over_axis(ex: Expansion, axis: int | None, param_list_len: int) -> list[Expansion]:
+    if axis is None:
+        return [ex for _ in range(param_list_len)]
+    axis_len = ex.shape[axis]
+    if axis_len != param_list_len:
+        raise ValueError(f'Parameter list has different length than the provided expansion axis, '
+                         f'got param_list_len={param_list_len} and axis_len={axis_len}.')
+    return [Expansion(ex.l_array, np.take(ex.coefs, i, axis=axis)) for i in range(axis_len)]
 
 
 SingleExpansionFn = Callable[[Expansion, ModelParams], T]
@@ -45,58 +56,38 @@ TwoExpansionsFn = Callable[[Expansion, Expansion, ModelParams], T]
 def parameter_map_single_expansion(f: SingleExpansionFn,
                                    match_expansion_axis_to_params: int = None) -> SingleExpansionFn:
 
-    meap = match_expansion_axis_to_params  # just a shorter variable name
-
     def mapped_f(ex: Expansion, params: ModelParams):
         params_list = unravel_params(params)
-        if meap is not None:
-            expansion_list = [Expansion(ex.l_array, np.take(ex.coefs, i, axis=meap)) for i in range(ex.shape[meap])]
-        else:
-            expansion_list = [ex for _ in params_list]
-
-        if not len(expansion_list) == len(params_list):
-            raise ValueError(f'Axis of expansion that is supposed to match params does not have the same length, got '
-                             f'len(unraveled params) = {len(params_list)} and '
-                             f'expansion.shape[{meap}] = {len(expansion_list)}')
+        expansion_list = unravel_expansion_over_axis(ex, match_expansion_axis_to_params, len(params_list))
 
         results = []
         for exp, prms in zip(expansion_list, params_list):
             results.append(f(exp, prms))
 
-        if meap is not None:
-            return np.array(results).swapaxes(0, meap)  # return the params-matched axis to where it belongs
+        if match_expansion_axis_to_params is not None:
+            # return the params-matched axis to where it belongs
+            return np.moveaxis(np.array(results), 0, match_expansion_axis_to_params)
         return np.squeeze(np.array(results))
 
     return mapped_f
 
 
 def parameter_map_two_expansions(f: TwoExpansionsFn,
-                                 match_expansion_axis_to_params: int = None) -> TwoExpansionsFn:
-
-    meap = match_expansion_axis_to_params  # just a shorter variable name
+                                 match_expansion_axis_to_params: int = None,
+                                 ) -> TwoExpansionsFn:
 
     def mapped_f(ex1: Expansion, ex2: Expansion, params: ModelParams):
         params_list = unravel_params(params)
-        if meap is not None:
-            expansion_list1 = [Expansion(ex1.l_array, np.take(ex1.coefs, i, axis=meap)) for i in range(ex1.shape[meap])]
-            expansion_list2 = [Expansion(ex2.l_array, np.take(ex2.coefs, i, axis=meap)) for i in range(ex2.shape[meap])]
-        else:
-            expansion_list1 = [ex1 for _ in params_list]
-            expansion_list2 = [ex2 for _ in params_list]
-
-        if not (len(expansion_list1) == len(params_list) and len(expansion_list2) == len(params_list)):
-            raise ValueError(f'Axis of at least one of the expansions that is supposed to match params '
-                             f'does not have the same length, got '
-                             f'len(unraveled params) = {len(params_list)} and '
-                             f'expansion1.shape[{meap}] = {len(expansion_list1)} and '
-                             f'expansion2.shape[{meap}] = {len(expansion_list2)}')
+        expansion_list1 = unravel_expansion_over_axis(ex1, match_expansion_axis_to_params, len(params_list))
+        expansion_list2 = unravel_expansion_over_axis(ex2, match_expansion_axis_to_params, len(params_list))
 
         results = []
         for exp1, exp2, prms in zip(expansion_list1, expansion_list2, params_list):
             results.append(f(exp1, exp2, prms))
 
-        if meap is not None:
-            return np.array(results).swapaxes(0, meap)  # return the params-matched axis to where it belongs
+        if match_expansion_axis_to_params is not None:
+            # return the params-matched axis to where it belongs
+            return np.moveaxis(np.array(results), 0, match_expansion_axis_to_params)
         return np.squeeze(np.array(results))
 
     return mapped_f

charged_shells/patch_size.py

@@ -15,7 +15,8 @@ def charge_patch_size(ex: Expansion, phi: float = 0, theta0: Array | float = 0,
 
 def potential_patch_size(ex: Expansion, params: ModelParams,
                          phi: float = 0, theta0: Array | float = 0, theta1: Array | float = np.pi / 2,
-                         match_expansion_axis_to_params: int = None):
+                         match_expansion_axis_to_params: int = None,
+                         skip_nozero_cases: bool = False):
 
     # this is more complicate to map over leading axes of the expansion as potential also depends on model parameters,
     # with some, such as kappaR, also being the parameters of the expansion in the first place. When mapping,
@@ -23,7 +24,13 @@ def potential_patch_size(ex: Expansion, params: ModelParams,
 
     @mapping.map_over_expansion
     def potential_zero(exp: Expansion, prms: ModelParams):
-        return bisect(lambda theta: potentials.charged_shell_potential(theta, phi, 1, exp, prms), theta0, theta1)
+        try:
+            return bisect(lambda theta: potentials.charged_shell_potential(theta, phi, 1, exp, prms), theta0, theta1)
+        except ValueError:
+            if skip_nozero_cases:
+                return 0.
+            else:
+                raise ValueError("Potential has no zero on the given interval.")
 
     return mapping.parameter_map_single_expansion(potential_zero, match_expansion_axis_to_params)(ex, params)
 

charged_shells/rotational_path.py

@@ -72,6 +72,19 @@ class PathEnergyPlot:
     def __post_init__(self):
         if not isinstance(self.dist, Array):
             self.dist = np.array([self.dist])
+        # we add 1 to match_expansion_axis as rotations take the new leading axis
+        if self.match_expansion_axis_to_params is not None:
+            self.match_expansion_axis_to_params += 1
+
+    @staticmethod
+    def plot_style(fig, ax, energy_units: interactions.EnergyUnit = 'kT'):
+        ax.axhline(y=0, c='k', linestyle=':')
+        ax.legend(fontsize=15)
+        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+        ax.set_xlabel('angle', fontsize=15)
+        ax.set_ylabel(f'U [{energy_units}]', fontsize=15)
+        fig.set_size_inches(plt.figaspect(0.5))
+        plt.tight_layout()
 
     def evaluate_energy(self):
         energy = []
@@ -118,14 +131,9 @@ class PathEnergyPlot:
         energy = self.evaluate_path(norm_euler_angles=norm_euler_angles)
         x_axis = self.rot_path.stack_x_axes()
 
-        fig, ax = plt.subplots()
-        ax.axhline(y=0, c='k', linestyle=':')
+        fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
         ax.plot(x_axis, np.squeeze(energy), label=labels)
-        ax.legend(fontsize=12)
-        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
-        ax.set_xlabel('angle', fontsize=13)
-        ax.set_ylabel('U', fontsize=13)
-        plt.tight_layout()
+        self.plot_style(fig, ax, energy_units=self.units)
         if save_as is not None:
             plt.savefig(save_as, dpi=600)
         plt.show()
@@ -135,14 +143,10 @@ class PathEnergyPlot:
         normalization = self.normalization(norm_euler_angles)
 
         fig, ax = plt.subplots()
-        ax.axhline(y=0, c='k', linestyle=':')
         for x_axis, energy in zip(self.rot_path.x_axis, energy_list):
             energy, norm = np.broadcast_arrays(energy, normalization)
             ax.plot(x_axis, energy / norm)
-        # ax.legend(fontsize=12)
-        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
-        ax.set_xlabel('angle', fontsize=13)
-        ax.set_ylabel('U', fontsize=13)
+        self.plot_style(fig, ax, energy_units=self.units)
         if save_as is not None:
             plt.savefig(save_as, dpi=600)
         plt.show()

charged_shells/units.json

@@ -0,0 +1,6 @@
+{
+  "distance": 1e-09,
+  "charge": 1.6021766340000001e-19,
+  "voltage": 1,
+  "concentrationM": 0.001
+}

charged_shells/units_and_constants.py

@@ -1,4 +1,6 @@
 from dataclasses import dataclass
+from pathlib import Path
+import json
 
 
 @dataclass
@@ -28,10 +30,8 @@ class Constants:
     e0: float
 
 
-base_units = {'distance': 1e-9,
-              'charge': 1.602176634 * 1e-19,
-              'voltage': 1,
-              'concentrationM': 1e-3}
+with open(Path(__file__).resolve().parent.joinpath('units.json')) as f:
+    base_units = json.load(f)
 
 
 UNITS = BaseUnits(**base_units)

config.json

@@ -0,0 +1,4 @@
+{
+  "emanuele_data":  "/home/andraz/ChargedShells/data_Emanuele/2024-01-14",
+  "figures": "/home/andraz/ChargedShells/Figures"
+}

tests/expansion_mapping_test.py

@@ -0,0 +1,54 @@
+import unittest
+import numpy as np
+from charged_shells import expansion, interactions, parameters, units_and_constants, mapping, potentials
+from functools import partial
+
+
+EPSILON_0 = units_and_constants.CONSTANTS.epsilon0
+
+
+class IsotropicTest(unittest.TestCase):
+
+    def setUp(self):
+        self.charge = 560
+        self.kappaR = np.array([0.1, 1, 3, 5, 10, 20])
+        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.ex2 = self.ex1.clone()
+        self.params = parameters.ModelParams(kappaR=self.kappaR, R=self.radius)
+
+    def test_potential(self):
+        theta = np.linspace(0, np.pi, 100)
+
+        def cs_potential_fn(ex, params):
+            return potentials.charged_shell_potential(theta, 0., self.dist / self.radius, ex, params)
+
+        mapped_cs_potential_fn = mapping.parameter_map_single_expansion(cs_potential_fn, 0)
+        cs_potential = mapped_cs_potential_fn(self.ex1, self.params)
+
+        ic_potential = []
+        for p in mapping.unravel_params(self.params):
+            ic_potential.append(potentials.inverse_patchy_particle_potential(theta, 2., 0, self.sigma0,
+                                                                             (0., 0.), p, lmax=10))
+        ic_potential = np.array(ic_potential)
+
+        np.testing.assert_almost_equal(cs_potential, ic_potential)
+
+    def test_interaction(self):
+
+        int_analytic = (self.charge ** 2 / (4 * np.pi * EPSILON_0 * self.params.epsilon) *
+                        (np.exp(self.params.kappaR) / (1 + self.params.kappaR)) ** 2 * np.exp(-self.params.kappa * self.dist) / self.dist)
+
+        energy_fn = mapping.parameter_map_two_expansions(partial(interactions.charged_shell_energy,
+                                                                 dist=self.dist/self.radius, units='eV'),
+                                                         0)
+        int_comp = energy_fn(self.ex1, self.ex2, self.params)
+
+        # print(int_analytic)
+        # print(int_comp)
+        # print(int_comp / int_analytic)
+
+        np.testing.assert_almost_equal(int_comp, int_analytic)
+

tests/interactions_test.py

@@ -25,6 +25,24 @@ def v22_distance_test():
     plt.show()
 
 
+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)
+    ex1 = ex0.clone()
+
+    ex1.rotate_euler(0, np.pi / 2, 0)
+
+    dist = np.linspace(2, 3.2, 100)
+    energy_array = np.zeros((dist.shape[0], len(sigma2)))
+    for i, d in enumerate(dist):
+        energy_array[i, ...] = interactions.charged_shell_energy(ex0, ex1, params, d)
+
+    plt.plot(dist, energy_array)
+    plt.show()
+
+
 def timing():
     params = ModelParams(R=150, kappaR=3)
     ex1 = expansion.MappedExpansionQuad(0.44, params.kappaR, 0.001, l_max=20)
@@ -49,5 +67,7 @@ def timing():
 
 if __name__ == '__main__':
 
-    v22_distance_test()
-    timing()
+    # v22_distance_test()
+    # timing()
+
+    quadrupole_variation_test()