2023-11-17

expansion.py

@@ -11,6 +11,7 @@ import matplotlib.pyplot as plt
 from typing import Callable, TypeVar
 from scipy.special import eval_legendre
 import plotly.graph_objects as go
+from pathlib import Path
 
 
 Array = np.ndarray
@@ -249,6 +250,21 @@ def plot_theta_profile(ex: Expansion, phi: float = 0, num: int = 100, theta_star
     plt.show()
 
 
+def plot_theta_profile_multiple(ex_list: list[Expansion], label_list, phi: float = 0, num: int = 100, theta_start: float = 0, theta_end: float = np.pi):
+    theta_vals = np.linspace(theta_start, theta_end, num)
+
+    fig, ax = plt.subplots()
+    for ex, label in zip(ex_list, label_list):
+        ax.plot(theta_vals, ex.charge_value(theta_vals, phi).T, label=label)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\theta$', fontsize=13)
+    ax.set_ylabel(r'$\sigma$', fontsize=13)
+    plt.legend(fontsize=12)
+    plt.tight_layout()
+    plt.savefig(Path("/home/andraz/ChargedShells/Figures/charge_shape_comparison.png"), dpi=600)
+    plt.show()
+
+
 def plot_charge_3d(ex: Expansion, num_theta=100, num_phi=100):
     theta = np.linspace(0, np.pi, num_theta)
     phi = np.linspace(0, 2 * np.pi, num_phi)
@@ -301,9 +317,9 @@ if __name__ == '__main__':
     # ex = Expansion(np.arange(3), np.array([1, -1, 0, 1, 2, 0, 3, 0, 2]))
     # ex = GaussianCharges(omega_k=np.array([[0, 0], [np.pi, 0]]), lambda_k=10, sigma1=0.001, l_max=10)
     ex = SphericalCap(np.array([[0, 0], [np.pi, 0]]), 0.5, 0.1, 70)
-    print(np.real(ex.coefs))
-    plot_theta_profile(ex, num=1000, theta_end=2 * np.pi, phi=0)
-    # plot_charge_3d(ex)
+    # print(np.real(ex.coefs))
+    # plot_theta_profile(ex, num=1000, theta_end=2 * np.pi, phi=0)
+    plot_charge_3d(ex)
 
     # new_coeffs = expansion_rotation(Quaternion(np.arange(20).reshape(5, 4)).normalized, ex.coeffs, ex.l_array)
     # print(new_coeffs.shape)

patch_shape_comparison.py

@@ -6,6 +6,9 @@ from parameters import ModelParams
 import functions as fn
 import matplotlib.pyplot as plt
 import scipy.special as sps
+from pathlib import Path
+import rotational_path
+import path_plot
 
 
 def point_to_gauss_map(sigma_m, a_bar, lbd, params: ModelParams):
@@ -15,18 +18,13 @@ def point_to_gauss_map(sigma_m, a_bar, lbd, params: ModelParams):
 
 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, theta0) - sps.eval_legendre(3, theta0))) / 1.7537
-
-
-def point_to_cap_map2(sigma_1, lbd, theta0):
-    return (10 * sigma_1 * lbd * fn.sph_bessel_i(2, lbd) /
-            (np.sinh(lbd) * (sps.eval_legendre(1, theta0) - sps.eval_legendre(3, theta0)))) / 1.7537
+            a_bar ** 2 / (sps.eval_legendre(1, np.cos(theta0)) - sps.eval_legendre(3, np.cos(theta0))))
 
 
 if __name__ == '__main__':
 
     target_patch_size = 0.9
-    params = ModelParams(R=150, kappaR=3)
+    params = ModelParams(R=150, kappaR=1)
     sigma_m = 0.001
 
     def fn1(x):
@@ -48,23 +46,38 @@ if __name__ == '__main__':
     ex_gauss = expansion.GaussianCharges(lambda_k=lbd, omega_k=np.array([[0, 0], [np.pi, 0]]), sigma1=gauss_sigma, l_max=30)
 
     cap_sigma = point_to_cap_map(sigma_m, a_bar, theta0, params)
-    # cap_sigma = point_to_cap_map2(gauss_sigma, lbd, theta0)
-    ex_cap = expansion.SphericalCap(theta0_k=theta0, sigma1=cap_sigma, omega_k=np.array([[0, 0], [np.pi, 0]]), l_max=50)
+    ex_cap = expansion.SphericalCap(theta0_k=theta0, sigma1=cap_sigma, omega_k=np.array([[0, 0], [np.pi, 0]]), l_max=30)
 
     theta = np.linspace(0, np.pi, 1000)
     phi = 0.
     dist = 1
 
-    potential1 = potentials.charged_shell_potential(theta, phi, dist, ex_point, params, print_idx=3)
-    potential2 = potentials.charged_shell_potential(theta, phi, dist, ex_gauss, params, print_idx=2)
-    potential3 = potentials.charged_shell_potential(theta, phi, dist, ex_cap, params, print_idx=6)
+    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)
     # print(potential.shape)
     # print(potential)
 
-    plt.plot(theta, potential1.T)
-    plt.plot(theta, potential2.T)
-    plt.plot(theta, potential3.T)
-    plt.show()
+    # expansion.plot_theta_profile_multiple([ex_point, ex_gauss, ex_cap], ['IC', 'Gauss', 'cap'], num=1000)
+
+    # fig, ax = plt.subplots()
+    # ax.plot(theta, potential1.T, label='IC')
+    # ax.plot(theta, potential2.T, label='Gauss')
+    # ax.plot(theta, potential3.T, label='cap')
+    # ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    # ax.set_xlabel(r'$\theta$', fontsize=13)
+    # ax.set_ylabel(r'$\phi$', fontsize=13)
+    # plt.legend(fontsize=12)
+    # plt.tight_layout()
+    # plt.savefig(Path("/home/andraz/ChargedShells/Figures/potential_shape_comparison.png"), dpi=600)
+    # 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"))
 
 
 

patch_size.py

@@ -74,21 +74,21 @@ if __name__ == '__main__':
     a_bar = np.linspace(0.2, 0.8, 100)
     kappaR = np.array([0.26, 1, 3, 10, 26])
     params = ModelParams(R=150, kappaR=10)
-    # ex = expansion.MappedExpansionQuad(a_bar=a_bar[:, None], sigma_m=0.001, l_max=30, kappaR=kappaR[None, :])
-
-    fn = lambda x: expansion.SphericalCap(theta0_k=x, sigma1=0.001, l_max=30, omega_k=np.array([[0, 0], [np.pi, 0]]))
-    print(inverse_potential_patch_size(48.8 * np.pi / 180, fn, 0.4, params))
-
-    # # patch_size = charge_patch_size(ex)
-    # patch_size = potential_patch_size(ex, params, match_expansion_axis_to_params=1)
-    #
-    # fig, ax = plt.subplots()
-    # for patch, kR in zip(patch_size.T, kappaR):
-    #     ax.plot(a_bar, patch * 180 / np.pi, label=rf'$\kappa$ = {kR}')
-    # ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
-    # ax.set_xlabel(r'$\bar a$', fontsize=13)
-    # ax.set_ylabel(r'$\theta$', fontsize=13)
-    # plt.legend(fontsize=12)
-    # # plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential.png"), dpi=600)
-    # plt.show()
+    ex = expansion.MappedExpansionQuad(a_bar=a_bar[:, None], sigma_m=0.001, l_max=30, kappaR=kappaR[None, :])
+
+    # fn = lambda x: expansion.SphericalCap(theta0_k=x, sigma1=0.001, l_max=30, omega_k=np.array([[0, 0], [np.pi, 0]]))
+    # print(inverse_potential_patch_size(48.8 * np.pi / 180, fn, 0.4, params))
+
+    # patch_size = charge_patch_size(ex)
+    patch_size = potential_patch_size(ex, params, match_expansion_axis_to_params=1)
+
+    fig, ax = plt.subplots()
+    for patch, kR in zip(patch_size.T, kappaR):
+        ax.plot(a_bar, patch * 180 / np.pi, label=rf'$\kappa$ = {kR}')
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\bar a$', fontsize=13)
+    ax.set_ylabel(r'$\theta$', fontsize=13)
+    plt.legend(fontsize=12)
+    # plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential.png"), dpi=600)
+    plt.show()
 

path_plot.py

@@ -52,26 +52,26 @@ if __name__ == '__main__':
     # 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, 15]), sigma1=0.001, l_max=20)
-    ex1 = expansion.SphericalCap(np.array([[0, 0], [np.pi, 0]]), np.array([0.3, 0.5, 1]), 0.001, 30)
+    ex1 = expansion.SphericalCap(np.array([[0, 0], [np.pi, 0]]), np.array([0.5]), 0.001, 20)
     ex2 = ex1.clone()
 
     # charge_profile = ex1.charge_value(theta=np.linspace(0, np.pi, 100), phi=0)
     # plt.plot(charge_profile.T)
     # plt.show()
 
-    expansion.plot_theta_profile(ex1, num=1000, theta_end=np.pi, phi=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)
 
-    ps = patch_size.potential_patch_size(ex1, params, theta1=np.pi / 2,
-                                         # match_expansion_axis_to_params=1
-                                         )
-    print(ps)
+    path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params)
 
-    # path_plot = PathEnergyPlot(ex1, ex2, QuadPath, dist=2., params=params)
-
-    # path_plot.plot(labels=[rf'$\theta$={180 * patch / np.pi:.1f}' for patch in np.array([0.3, 0.5, 1.0])],
-    #                # norm_euler_angles={'beta2': np.pi},
-    #                # save_as=Path("/home/andraz/ChargedShells/Figures/dipole_path_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")
+                   )
 
     # path_plot.plot(labels=[rf'$\kappa R$={kR}' for kR in kappaR],
     #                # norm_euler_angles={'beta2': np.pi},

peak_heigth.py

@@ -0,0 +1,54 @@
+import matplotlib.pyplot as plt
+import expansion
+import interactions
+from parameters import ModelParams
+import numpy as np
+from typing import Literal
+from pathlib import Path
+
+Array = np.ndarray
+Expansion = expansion.Expansion
+RotConfig = Literal['ep', 'pp']
+
+
+def peak_energy_plot(kappaR: Array,
+                     a_bar: Array,
+                     which: Literal['ep', 'pp'],
+                     R: float = 150,
+                     dist: float = 2.,
+                     l_max=20,
+                     save_as: Path = None):
+
+    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)
+
+    fig, ax = plt.subplots()
+    for en, ab in zip(energy.T, a_bar):
+        ax.plot(kappaR, en / en[0], label=rf'$\bar a = {ab:.1f}$')
+    ax.legend(fontsize=12)
+    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
+    ax.set_xlabel(r'$\kappa R$', fontsize=13)
+    ax.set_ylabel(rf'$\bar V$', fontsize=13)
+    plt.tight_layout()
+    if save_as is not None:
+        plt.savefig(save_as, dpi=600)
+    plt.show()
+
+
+if __name__ == '__main__':
+
+    kappaR = np.arange(1, 10, 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_pp.pdf')
+                     )

rotational_path.py

@@ -59,6 +59,7 @@ class PairRotationalPath:
 
 @dataclass
 class PathEnergyPlot:
+    """Path comparison for a pair of charge distributions, possibly at different parameter values."""
 
     expansion1: Expansion
     expansion2: Expansion
@@ -102,11 +103,15 @@ class PathEnergyPlot:
                                      gamma=norm_euler_angles.get('gamma2', 0))
         return np.abs(self.evaluate_energy())
 
-    def plot(self, labels: list[str] = None, norm_euler_angles: dict = None, save_as: Path = None):
+    def evaluate_path(self, norm_euler_angles: dict = None):
         energy = self.path_energy()
         normalization = self.normalization(norm_euler_angles)
         energy = energy / normalization[None, ...]
         energy = energy.reshape(energy.shape[0], -1)
+        return np.squeeze(energy)
+
+    def plot(self, labels: list[str] = None, norm_euler_angles: dict = None, save_as: Path = None):
+        energy = self.evaluate_path(norm_euler_angles=norm_euler_angles)
         x_axis = self.rot_path.stack_x_axes()
 
         fig, ax = plt.subplots()
@@ -137,3 +142,41 @@ class PathEnergyPlot:
         if save_as is not None:
             plt.savefig(save_as, dpi=600)
         plt.show()
+
+
+@dataclass
+class PathExpansionComparison:
+    """Path comparison for different charge distribution models."""
+
+    ex_list: list[Expansion]
+    rot_path: PairRotationalPath
+    dist: float
+    params: ModelParams
+    path_list: list[PathEnergyPlot] = field(default_factory=list)
+
+    def __post_init__(self):
+        for ex in self.ex_list:
+            self.path_list.append(PathEnergyPlot(ex, ex.clone(), self.rot_path, self.dist, self.params))
+
+    def evaluate_path(self, norm_euler_angles: dict = None):
+        path_vals = []
+        for path in self.path_list:
+            path_vals.append(path.evaluate_path(norm_euler_angles))
+        return np.stack(path_vals)
+
+    def plot(self, labels: list[str] = None, norm_euler_angles: dict = None, save_as: Path = None):
+        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=':')
+        ax.plot(x_axis, np.squeeze(energy).T, 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()
+        if save_as is not None:
+            plt.savefig(save_as, dpi=600)
+        plt.show()
+