analysis/energy_gap.py

@@ -74,30 +74,34 @@ def abar_kappaR_dependence2(save_as=None):
 def charge_kappaR_dependence(a_bar, min_charge, max_charge, save_as=None, cmap=cm.jet):
 
     kappaR = np.linspace(0.01, 10, 50)
+    sigma_tilde = 0.001
     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, sigma_tilde, 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 = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=np.min(charge) / sigma_tilde,
+                                                         vmax=np.max(charge) / sigma_tilde))
     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)
+    for g, c in zip(gap.T, colors):
+        ax.plot(kappaR, g, 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.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+    ax.set_xlabel(r'$\kappa R$', fontsize=20)
+    ax.set_ylabel(r'$(V_{pp}-V_{ep})/V_{pp}$', fontsize=20)
 
     divider = make_axes_locatable(ax)
     cax = divider.append_axes('right', size='5%', pad=0.05)
+    cax.tick_params(labelsize=15)
     cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
-    cbar.set_label(r'$\tilde \sigma_0$', rotation=90, labelpad=15, fontsize=12)
+    cbar.set_label(r'$\eta$', rotation=0, labelpad=15, fontsize=20)
 
-    plt.tight_layout()
+    # plt.tight_layout()
+    plt.subplots_adjust(left=0.1, right=0.9, top=0.95, bottom=0.12)
     if save_as is not None:
         plt.savefig(save_as, dpi=600)
     plt.show()
@@ -193,22 +197,42 @@ 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, config_data: dict, 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 = np.load(em_data_path.joinpath("relative_gap.npz"))
+    em_data = np.load(em_data_path.joinpath("relative_gap_ZC.npz"))
     data = em_data[which_change]
 
+    sigma_tilde = 0.001
+
     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)
+    # print(f'All charge: {charge}')
 
-    colors = cmap(np.linspace(0, 1, len(charge)))
-    sm = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=np.min(charge), vmax=np.max(charge)))
+    eta = charge / sigma_tilde
+
+    if eta_min is None:
+        eta_min = np.min(eta)
+    if eta_max is None:
+        eta_max = np.max(eta)
+
+    def map_eta_to_unit(x):
+        return (x - eta_min) / (eta_max - eta_min)
+
+    # print(eta[0], eta[1])
+    # print(map_eta_to_unit(eta[0]), map_eta_to_unit(eta[-1]))
+
+    colors_linspace = np.linspace(map_eta_to_unit(eta[0]), map_eta_to_unit(eta[-1]), len(charge))
+    colors_linspace[colors_linspace > 1] = 1
+    colors_linspace[colors_linspace < 0] = 0
+
+    colors = cmap(colors_linspace)
+    sm = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=eta_min, vmax=eta_max))
     sm.set_array([])
 
     fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
@@ -218,33 +242,44 @@ def IC_gap_charge_at_abar(a_bar, config_data: dict, save_as=None, cmap=cm.coolwa
         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.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+    ax.set_xlabel(r'$\kappa R$', fontsize=20)
+    ax.set_ylabel(r'$(V_{pp}-V_{ep})/V_{pp}$', fontsize=20)
     ax.set_xlim(-0.25, 10.25)
-    # ax.set_ylim(-4.25, 6.25)
+    ax.set_ylim(-4, 5)
 
     divider = make_axes_locatable(ax)
     cax = divider.append_axes('right', size='5%', pad=0.05)
+    cax.tick_params(labelsize=15)
     cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
-    cbar.set_label(r'$\tilde \sigma_0$', rotation=90, labelpad=15, fontsize=12)
+    cbar.set_label(r'$\eta$', rotation=0, labelpad=15, fontsize=20)
 
-    plt.tight_layout()
+    # plt.tight_layout()
+    plt.subplots_adjust(left=0.1, right=0.9, top=0.95, bottom=0.12)
     if save_as is not None:
         plt.savefig(save_as, dpi=600)
     plt.show()
 
 
+def test_gap(a_bar, kappaR, charge):
+    params = ModelParams(R=150, kappaR=kappaR)
+    ex = expansion.MappedExpansionQuad(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)
+
     # 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"),
+    # charge_kappaR_dependence(a_bar=0.8, min_charge=-0.002, max_charge=0.002,
+    #                          save_as=Path("/home/andraz/ChargedShells/Figures/full_amplitude_charge_abar08.png"),
     #                          cmap=cm.coolwarm)
 
     # charge_kappaR_dependence_heatmap(a_bar=0.5, min_charge=-0.003, max_charge=0.003,
@@ -256,8 +291,8 @@ def main():
     # 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")
+    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")
                           )
 
 

analysis/expansion_plot.py

@@ -1,8 +1,10 @@
-from charged_shells import expansion
+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
+import quadrupole_model_mappings
 
 Expansion = expansion.Expansion
 
@@ -30,7 +32,7 @@ def plot_theta_profile_multiple(ex_list: list[Expansion], label_list, phi: float
     plt.show()
 
 
-def plot_charge_3d(ex: Expansion, num_theta=100, num_phi=100):
+def plot_charge_3d(ex: Expansion, num_theta=100, num_phi=100, save_as: Path = None):
     theta = np.linspace(0, np.pi, num_theta)
     phi = np.linspace(0, 2 * np.pi, num_phi)
     theta, phi = np.meshgrid(theta, phi)
@@ -42,21 +44,40 @@ def plot_charge_3d(ex: Expansion, num_theta=100, num_phi=100):
     z = np.cos(theta)
 
     # Create a heatmap on the sphere
-    fig = go.Figure(data=go.Surface(x=x, y=y, z=z, surfacecolor=r, colorscale='Jet'))
+    fig = go.Figure(data=go.Surface(x=x, y=y, z=z, surfacecolor=r,
+                                    colorscale='RdBu', reversescale=True))
     fig.update_layout(scene=dict(aspectmode='data'))
-    fig.update_layout(scene=dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z'))
+    fig.update_layout(scene=dict(xaxis_title='', yaxis_title='', zaxis_title=''))
+
+    # Remove axes planes, background, ticks, and labels
+    fig.update_layout(scene=dict(xaxis=dict(showbackground=False, gridcolor='white', showticklabels=False, ticks=''),
+                                 yaxis=dict(showbackground=False, gridcolor='white', showticklabels=False, ticks=''),
+                                 zaxis=dict(showbackground=False, gridcolor='white', showticklabels=False, ticks='')))
+
+    # Adjust the width and height for higher resolution
+    fig.update_layout(width=1200, height=1200)
+
+    # Save as PNG with higher resolution
+    if save_as is not None:
+        fig.write_image(save_as, scale=3)  # Adjust the scale as needed
 
     fig.show()
 
 
 def main():
-    # ex = expansion.MappedExpansionQuad(0.44, 3, 1, 10)
+    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)
     # ex = expansion.Expansion(np.arange(3), np.array([1, -1, 0, 1, 2, 0, 3, 0, 2]))
-    # ex = expansion.GaussianCharges(omega_k=np.array([[0, 0], [np.pi, 0]]), lambda_k=10, sigma1=0.001, l_max=10)
-    ex = expansion.SphericalCap(np.array([[0, 0], [np.pi, 0]]), 0.5, 0.1, 30)
+    # ex = expansion.GaussianCharges(omega_k=np.array([[0, 0], [np.pi, 0]]), lambda_k=2.676, sigma1=0.00044, l_max=30)
+    # ex = expansion.SphericalCap(np.array([[0, 0], [np.pi, 0]]), 0.894, 0.00132, 50)
+    # ex = quadrupole_model_mappings.ic_to_gauss(0.001, 0.328, params, l_max=30)
+    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)
+    plot_charge_3d(ex, save_as=Path(config_data["figures"]).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/patch_shape_comparison.py

@@ -1,21 +1,11 @@
-from charged_shells import expansion, functions as fn, potentials, patch_size
+from charged_shells import expansion, potentials, patch_size
 import numpy as np
 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):
-    return (sigma_m * fn.coefficient_Cim(2, params.kappaR) / fn.coefficient_Cpm(2, params.kappaR)
-            * np.sinh(lbd) / (lbd * fn.sph_bessel_i(2, lbd)) * a_bar ** 2)
-
-
-def point_to_cap_map(sigma_m, a_bar, theta0, params: ModelParams):
-    return (sigma_m * 10 * fn.coefficient_Cim(2, params.kappaR) / fn.coefficient_Cpm(2, params.kappaR) *
-            a_bar ** 2 / (sps.eval_legendre(1, np.cos(theta0)) - sps.eval_legendre(3, np.cos(theta0))))
+import json
+import quadrupole_model_mappings
 
 
 def ic_cs_comparison():
@@ -54,36 +44,21 @@ def ic_cs_comparison():
     plt.show()
 
 
-def models_comparison():
+def models_comparison(save_data=False):
     target_patch_size = 0.92
     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):
+    def fn(x):
         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, 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]]), 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)
-
+    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)
 
-    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, 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, sigma0=sigma0_mapped)
+    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)
 
     theta = np.linspace(0, np.pi, 1001)
     phi = 0.
@@ -97,6 +72,15 @@ def models_comparison():
     # print(potential.shape)
     # print(potential)
 
+    if save_data:
+        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_shape_models.npz"),
+                 ICi=np.stack((theta, potential_ic)).T,
+                 CSp_mapped=np.stack((theta, potential1)).T,
+                 CSp_gauss=np.stack((theta, potential2)).T,
+                 CSp_caps=np.stack((theta, potential3)).T)
+
     # expansion.plot_theta_profile_multiple([ex_point, ex_gauss, ex_cap], ['IC', 'Gauss', 'cap'], num=1000)
 
     fig, ax = plt.subplots()
@@ -188,7 +172,7 @@ def charge_comparsion():
 # TODO: comparison of patch shape at different kappaR and a neutral particle, with fixed patch size
 
 
-def IC_CS_comparison():
+def ic_cs_comparison2(save_data=False):
     target_patch_size = 0.92
     kappaR = 3
     params = ModelParams(R=150, kappaR=kappaR)
@@ -215,17 +199,30 @@ def IC_CS_comparison():
                                                                     (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:
+            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,
+                 ICi_eta_0=np.stack((theta, potential_ic[1])).T,
+                 ICi_eta_pos=np.stack((theta, potential_ic[2])).T,
+                 CSp_eta_neg=np.stack((theta, potential_cs[0])).T,
+                 CSp_eta_0=np.stack((theta, potential_cs[1])).T,
+                 CSp_eta_pos=np.stack((theta, potential_cs[2])).T)
+
     fig, ax = plt.subplots()
     lines = []
+    lines.append(plt.scatter([], [], marker='x', c='k', label='CSp'))
+    lines.append(plt.scatter([], [], marker='o', facecolor='none', edgecolor='k', label='ICi'))
     for p_ic, p_cs, charge in zip(potential_ic, potential_cs, sigma0_array):
-        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)
+        l, = ax.plot(theta, 1000 * p_cs.T, label=rf'$\eta = {charge/sigma_m}$', linewidth=3)
         lines.append(l)
+        ax.scatter(theta[::50], 1000 * p_ic.T[::50], marker='o', facecolors='none', edgecolors='k', s=75)
+        ax.scatter(theta[::50], 1000 * p_cs.T[::50], marker='x', c='k', s=75)
     # 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)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+    ax.set_xlabel(r'$\theta$', fontsize=20)
+    ax.set_ylabel(r'$\phi$ [mV]', fontsize=20)
     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=':')
@@ -242,9 +239,9 @@ def IC_CS_comparison():
 
 
 def main():
-    # models_comparison()
+    models_comparison(save_data=False)
     # ic_cs_comparison()
-    IC_CS_comparison()
+    # ic_cs_comparison2()
 
     # abar_comparison()
     # charge_comparsion()

analysis/patch_size_dependence.py

@@ -4,9 +4,10 @@ import numpy as np
 from charged_shells import expansion, patch_size
 from charged_shells.parameters import ModelParams
 import time
+import json
 
 
-def plot_abar_dependence(*, save=False):
+def plot_abar_dependence(*, save=False, save_data=False):
     a_bar = np.linspace(0.2, 0.8, 101)
     kappaR = np.array([0.1, 1, 3, 10, 30])
     params = ModelParams(R=150, kappaR=kappaR)
@@ -16,14 +17,22 @@ def plot_abar_dependence(*, save=False):
 
     markers = ['o', 's', '^', 'D', 'v', '<', '>', 'p', '*', 'H', '+', 'x']
 
+    if save_data:
+        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)
+
     # 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)
-    ax.set_xlabel(r'$\bar a$', fontsize=15)
-    ax.set_ylabel(r'$\theta_0$', fontsize=15)
-    plt.legend(fontsize=14)
+        ax.plot(a_bar, patch, label=rf'$\kappa R$ = {kR}', marker=marker, markerfacecolor='none', markevery=5, linewidth=2.5, markersize=8)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+    ax.set_xlabel(r'$\bar a$', fontsize=20)
+    ax.set_ylabel(r'$\theta_0$', fontsize=20)
+    plt.legend(fontsize=18)
     plt.tight_layout()
     if save:
         plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential_new.pdf"), dpi=300)

analysis/path_plot.py

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

analysis/peak_heigth.py

@@ -1,11 +1,12 @@
 import matplotlib.pyplot as plt
-from charged_shells import expansion, interactions, mapping
+from charged_shells import expansion, interactions, mapping, functions
 from charged_shells.parameters import ModelParams
 import numpy as np
 from typing import Literal
 from pathlib import Path
 from functools import partial
 import json
+from itertools import cycle
 
 Array = np.ndarray
 Expansion = expansion.Expansion
@@ -22,16 +23,6 @@ def peak_energy_plot(kappaR: Array,
 
     params = ModelParams(R=R, kappaR=kappaR)
 
-    # energy = []
-    # for params in params.unravel():
-    #     ex1 = expansion.MappedExpansionQuad(a_bar, params.kappaR, 0.001, l_max=l_max)
-    #     ex2 = ex1.clone()
-    #     if which == 'ep':
-    #         ex1.rotate_euler(alpha=0, beta=np.pi / 2, gamma=0)
-    #     energy.append(interactions.charged_shell_energy(ex1, ex2, dist, params))
-    #
-    # energy = np.array(energy)
-
     ex1 = expansion.MappedExpansionQuad(a_bar[:, None], params.kappaR[None, :], 0.001, l_max=l_max)
     ex2 = ex1.clone()
     if which == 'ep':
@@ -56,8 +47,12 @@ def peak_energy_plot(kappaR: Array,
 
 def IC_peak_energy_plot(config_data: dict,
                         a_bar: list,
-        which: Literal['ep', 'pp'],
-                     save_as: Path = None):
+                        which: RotConfig,
+                        max_kappaR: float = 30.,
+                        R: float = 150,
+                        save_as: Path = None,
+                        save_data: bool = False,
+                        quad_correction: bool = False):
 
     em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_11"))
     em_data = np.load(em_data_path.joinpath("pair_energy_fig11.npz"))
@@ -71,50 +66,100 @@ def IC_peak_energy_plot(config_data: dict,
         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]])
+    all_kappaR = np.unique(data[:, 2])
 
-    # energies = data[:, column_idx].reshape(-1, counts[0])
-    # kappaR = data[:, 2].reshape(-1, counts[0])
-    # print(len(kappaR), len(energies), len(abar))
+    kappaR = np.geomspace(np.min(all_kappaR), max_kappaR, 30)
+    params = ModelParams(R=R, kappaR=kappaR)
 
-    fig, ax = plt.subplots()
+    ex1 = expansion.MappedExpansionQuad(np.array(a_bar)[:, None], params.kappaR[None, :], 0.001, l_max=20)
+    ex2 = ex1.clone()
+    if which == 'ep':
+        ex1.rotate_euler(alpha=0, beta=np.pi / 2, gamma=0)
+
+    energy_fn = mapping.parameter_map_two_expansions(partial(interactions.charged_shell_energy, dist=2), 1)
+    energy = energy_fn(ex1, ex2, params)
+
+    data_dict_ic = {}
+    data_dict_cs = {}
+    k = 0
     for i in range(len(abar)):
-        if abar[i] in a_bar:
+        ab = np.around(abar[i], 5)
+        if ab in np.around(a_bar, 5):
             idx, = np.nonzero(inverse == i)
-            kR = data[idx, 2]
-            en = data[idx, column_idx]
+            kR = data[idx, 2][data[idx, 2] <= max_kappaR]
+            en = data[idx, column_idx][data[idx, 2] <= max_kappaR]
+            if quad_correction:
+                en *= extra_correction(ab, kR)
             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)
+            data_dict_ic[ab] = np.stack((kR[sort], np.abs(en)[sort])).T
+            data_dict_cs[ab] = np.stack((kappaR, np.abs(energy[k]))).T
+            k += 1
+
+    if save_data:
+        ic_save = {str(key): val for key, val in data_dict_ic.items()}
+        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)
+
+    colors = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
+
+    fig, ax = plt.subplots()
+    for (key, data1), data2 in zip(data_dict_ic.items(), data_dict_cs.values()):
+        current_color = next(colors)
+        ax.plot(data1[:, 0], data1[:, 1], label=rf'$\bar a = {key:.2f}$', c=current_color)
+        ax.plot(data2[:, 0], data2[:, 1], ls='--', c=current_color)
+    ax.legend(fontsize=14, ncol=2)
     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()
+
+    y_label = r'$U_{pp}$' if which == 'pp' else r'$|U_{ep}|$'
+    ax.set_ylabel(y_label, fontsize=15)
     ax.set_yscale('log')
     ax.set_xscale('log')
+    plt.tight_layout()
     if save_as is not None:
         plt.savefig(save_as, dpi=600)
     plt.show()
 
 
-if __name__ == '__main__':
-    
+def extra_correction(a_bar, kappaR):
+    y = a_bar * kappaR
+    # return (kappaR ** 2 * a_bar ** 4 * functions.sph_bessel_k(1, kappaR) / functions.sph_bessel_k(3, kappaR) *
+    #         np.sinh(y) / (-3 * y * np.cosh(y) + (3 + y ** 2) * np.sinh(y)))
+    return (a_bar ** 2 * functions.sph_bessel_k(1, kappaR) / functions.sph_bessel_k(3, kappaR)
+            * functions.sph_bessel_i(0, y) / functions.sph_bessel_i(2, y))
+
+
+def main():
     with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
         config_data = json.load(config_file)
 
-    kappaR = np.arange(0.5, 10, 0.1)
-    a_bar = np.arange(0.2, 0.8, 0.2)
-
+    # kappaR = np.arange(0.5, 50, 0.1)
+    # a_bar = np.arange(0.2, 0.8, 0.2)
     # peak_energy_plot(kappaR, a_bar, which='pp',
     #                  # save_as=Path('/home/andraz/ChargedShells/Figures/nonmonotonicity_check_ep.pdf')
     #                  )
+
+    a_bar = [0.2, 0.32, 0.52, 0.8]
+    # a_bar = list(np.arange(0.2, 0.81, 0.04))
+    IC_peak_energy_plot(config_data, a_bar=a_bar, which='pp', max_kappaR=50,
+                        save_as=Path('/home/andraz/ChargedShells/Figures/Emanuele_data/both_nonmonotonicity_check_pp.png'),
+                        save_data=False,
+                        quad_correction=False
+                        )
+
+    # kr = np.linspace(0.1, 100, 1000)
+    # for ab in [0.2, 0.32, 0.52, 0.8]:
+    #     plt.plot(kr, extra_correction(ab, kr))
+    # plt.xscale('log')
+    # plt.yscale('log')
+    # plt.show()
+
+
+if __name__ == '__main__':
     
-    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')
-                     )
+    main()

analysis/quadrupole_model_mappings.py

@@ -0,0 +1,53 @@
+from charged_shells import patch_size, expansion, parameters
+from charged_shells import functions as fn
+from scipy.special import eval_legendre
+import numpy as np
+
+
+ModelParams = parameters.ModelParams
+Array = np.ndarray
+
+
+def point_to_gauss_magnitude(sigma_tilde: Array, a_bar: Array, lbd: Array, kappaR: Array):
+    return (sigma_tilde * fn.coefficient_Cim(2, kappaR) / fn.coefficient_Cpm(2, kappaR)
+            * np.sinh(lbd) / (lbd * fn.sph_bessel_i(2, lbd)) * a_bar ** 2)
+
+
+def point_to_cap_magnitude(sigma_tilde: Array, a_bar: Array, theta0: Array, kappaR: Array):
+    return (sigma_tilde * 10 * fn.coefficient_Cim(2, kappaR) / fn.coefficient_Cpm(2, kappaR) *
+            a_bar ** 2 / (eval_legendre(1, np.cos(theta0)) - eval_legendre(3, np.cos(theta0))))
+
+
+def ic_to_gauss(sigma_tilde, a_bar, params: ModelParams, l_max: int = 30,
+                sigma0: float = 0) -> expansion.GaussianCharges:
+
+    ex_mapped = expansion.MappedExpansionQuad(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)
+
+    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)
+
+    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)
+
+
+def ic_to_cap(sigma_tilde, a_bar, params: ModelParams, l_max: int = 30, sigma0: float = 0) -> expansion.SphericalCap:
+
+    ex_mapped = expansion.MappedExpansionQuad(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)
+
+    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)
+
+    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)

analysis/sEE_minimum.py

@@ -87,22 +87,23 @@ def kappaR_dependence(kappaR, save_as=None, cmap=cm.jet):
     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))
+    labelLines(ax.get_lines(), align=False, fontsize=15,  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.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+    ax.set_xlabel('pry angle', fontsize=20)
     ax.set_ylabel('U', fontsize=20)
     ax.set_xlim(2.2, 2.65)
-    ax.set_ylim(-20, 1)
+    ax.set_ylim(-31, 1)
 
     divider = make_axes_locatable(ax)
     cax = divider.append_axes('right', size='5%', pad=0.05)
+    cax.tick_params(labelsize=15)
     cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
-    cbar.set_label(r'$\bar a$', rotation=90, labelpad=15, fontsize=15)
+    cbar.set_label(r'$\bar a$', rotation=0, labelpad=15, fontsize=20)
 
     plt.tight_layout()
     if save_as is not None:
@@ -147,22 +148,23 @@ def IC_kappaR_dependence(config_data: dict, which_kappa_lines: list,
         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))
+    labelLines(ax.get_lines(), align=False, fontsize=15, 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.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+    ax.set_xlabel('pry angle', fontsize=20)
     ax.set_ylabel('U', fontsize=20)
     ax.set_xlim(2.2, 2.65)
-    ax.set_ylim(-40, 1)
+    ax.set_ylim(-41, 1)
 
     divider = make_axes_locatable(ax)
     cax = divider.append_axes('right', size='5%', pad=0.05)
+    cax.tick_params(labelsize=15)
     cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
-    cbar.set_label(r'$\bar a$', rotation=90, labelpad=15, fontsize=12)
+    cbar.set_label(r'$\bar a$', rotation=0, labelpad=15, fontsize=20)
 
     plt.tight_layout()
     if save_as is not None:
@@ -211,13 +213,13 @@ def main():
 
     # contours()
 
-    kappaR_dependence(kappaR=np.linspace(0.1, 30, 20),
-                      save_as=Path(config_data["figures"]).joinpath('sEE_min_kappaR_abar.png')
-                      )
+    # 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')
-    #                      )
+    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')
+                         )
 
 
 if __name__ == '__main__':

charged_shells/functions.py

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

charged_shells/rotational_path.py

@@ -19,18 +19,20 @@ class PairRotationalPath:
     rotations1: list[Quaternion] = field(default_factory=list)
     rotations2: list[Quaternion] = field(default_factory=list)
     x_axis: list[Array] = field(default_factory=list)
+    ticks: dict[float, str] = field(default_factory=dict)
     overlapping_last: bool = True
     _default_x_axis: Array = None
 
-    def add(self, rotation1: Quaternion, rotation2: Quaternion, x_axis: Array = None):
+    def add(self, rotation1: Quaternion, rotation2: Quaternion, x_axis: Array = None,
+            start_name: str | float = None, end_name: str | float = None):
         rotation1, rotation2 = np.broadcast_arrays(rotation1, rotation2)
         self.rotations1.append(Quaternion(rotation1))
         self.rotations2.append(Quaternion(rotation2))
         if x_axis is None:
             x_axis = np.arange(len(rotation1)) if self._default_x_axis is None else self._default_x_axis
-        self.add_x_axis(x_axis)
+        self.add_x_axis(x_axis, start_name, end_name)
 
-    def add_x_axis(self, x_axis: Array):
+    def add_x_axis(self, x_axis: Array, start_name: str | float = None, end_name: str | float = None):
         try:
             last_x_val = self.x_axis[-1][-1]
         except IndexError:
@@ -39,16 +41,23 @@ class PairRotationalPath:
             self.x_axis.append(x_axis + last_x_val)
         else:
             raise NotImplementedError('Currently only overlapping end points for x-axes are supported.')
+        # adding ticks to x_axis
+        s_n = x_axis[0] + last_x_val if start_name is None else start_name  # defaults to just numbers
+        e_n = x_axis[-1] + last_x_val if end_name is None else end_name
+        start_position = float(self.x_axis[-1][0])
+        end_position = float(self.x_axis[-1][-1])
+        if start_position not in self.ticks or start_name is not None:  # allows overwriting previous end with new start
+            self.ticks[start_position] = s_n
+        self.ticks[end_position] = e_n
 
     def set_default_x_axis(self, default_x_axis: Array):
         self._default_x_axis = default_x_axis
 
     def add_euler(self, *, alpha1: Array = 0, beta1: Array = 0, gamma1: Array = 0,
-                  alpha2: Array = 0, beta2: Array = 0, gamma2: Array = 0,
-                  x_axis: Array = None):
+                  alpha2: Array = 0, beta2: Array = 0, gamma2: Array = 0, **add_kwargs):
         R1_euler = quaternionic.array.from_euler_angles(alpha1, beta1, gamma1)
         R2_euler = quaternionic.array.from_euler_angles(alpha2, beta2, gamma2)
-        self.add(Quaternion(R1_euler), Quaternion(R2_euler), x_axis)
+        self.add(Quaternion(R1_euler), Quaternion(R2_euler), **add_kwargs)
 
     def stack_rotations(self) -> (Quaternion, Quaternion):
         return Quaternion(np.vstack(self.rotations1)), Quaternion(np.vstack(self.rotations2))
@@ -56,6 +65,16 @@ class PairRotationalPath:
     def stack_x_axes(self) -> Array:
         return np.hstack(self.x_axis)
 
+    def plot_style(self, fig, ax, energy_units: interactions.EnergyUnit = 'kT'):
+        ax.axhline(y=0, c='k', linestyle=':')
+        ax.legend(fontsize=18)
+        ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
+        # ax.set_xlabel('angle', fontsize=15)
+        ax.set_ylabel(f'U [{energy_units}]', fontsize=20)
+        ax.set_xticks(list(self.ticks.keys()), list(self.ticks.values()), fontsize=18)
+        fig.set_size_inches(plt.figaspect(0.5))
+        plt.tight_layout()
+
 
 @dataclass
 class PathEnergyPlot:
@@ -76,15 +95,8 @@ class PathEnergyPlot:
         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 plot_style(self, fig, ax, energy_units: interactions.EnergyUnit = 'kT'):
+        self.rot_path.plot_style(fig, ax, energy_units)
 
     def evaluate_energy(self):
         energy = []

config.json

@@ -1,4 +1,5 @@
 {
-  "emanuele_data":  "/home/andraz/ChargedShells/data_Emanuele/2024-01-14",
-  "figures": "/home/andraz/ChargedShells/Figures"
+  "emanuele_data":  "/home/andraz/ChargedShells/data_Emanuele/2024-02-15",
+  "figures": "/home/andraz/ChargedShells/Figures",
+  "figure_data": "/home/andraz/ChargedShells/Figures/Data"
 }

tests/expansion_mapping_test.py

@@ -46,9 +46,9 @@ class IsotropicTest(unittest.TestCase):
                                                          0)
         int_comp = energy_fn(self.ex1, self.ex2, self.params)
 
-        # print(int_analytic)
-        # print(int_comp)
-        # print(int_comp / int_analytic)
+        print(int_analytic)
+        print(int_comp)
+        print(int_comp / int_analytic)
 
         np.testing.assert_almost_equal(int_comp, int_analytic)