2023-11-26

patch_shape_comparison.py

@@ -48,10 +48,11 @@ if __name__ == '__main__':
     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)
 
-    theta = np.linspace(0, np.pi, 1000)
+    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)
     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)
@@ -60,24 +61,31 @@ if __name__ == '__main__':
 
     # 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()
+    fig, ax = plt.subplots()
+    ax.plot(theta, 1000 * potential1.T, label='CS', linewidth=2)
+    # ax.plot(theta, potential_ic.T, label='IC', ls=':', linewidth=2, marker='o', markevery=50, mfc='none')
+    ax.scatter(theta[::50], 1000 * potential_ic.T[::50], marker='o', label='IC', facecolors='none', edgecolors='tab:red')
+    ax.plot(theta, 1000 * potential2.T, label='Gauss', linewidth=2, ls='--')
+    ax.plot(theta, 1000 * potential3.T, label='cap', linewidth=2, ls='--')
+    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)
+    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"))
+    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

@@ -71,9 +71,9 @@ def inverse_potential_patch_size(target_patch_size: float,
 
 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)
+    a_bar = np.linspace(0.2, 0.8, 101)
+    kappaR = np.array([0.1, 1, 3, 10, 30])
+    params = ModelParams(R=150, kappaR=kappaR)
     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]]))
@@ -82,12 +82,14 @@ if __name__ == '__main__':
     # patch_size = charge_patch_size(ex)
     patch_size = potential_patch_size(ex, params, match_expansion_axis_to_params=1)
 
+    markers = ['o', 's', '^', 'D', 'v', '<', '>', 'p', '*', 'H', '+', 'x']
+
     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}')
+    for patch, kR, marker in zip(patch_size.T, kappaR, markers):
+        ax.plot(a_bar, patch, label=rf'$\kappa$ = {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=13)
-    ax.set_ylabel(r'$\theta$', fontsize=13)
+    ax.set_ylabel(r'$\theta_0$', fontsize=13)
     plt.legend(fontsize=12)
     # plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential.png"), dpi=600)
     plt.show()

potentials.py

@@ -4,6 +4,7 @@ import numpy as np
 import parameters
 import units_and_constants as uc
 import matplotlib.pyplot as plt
+from scipy.special import eval_legendre, kv
 
 
 Array = np.ndarray
@@ -50,21 +51,45 @@ def charged_shell_potential(theta: Array | float,
     return pot.reshape(ex.shape + angles_shape)
 
 
+def inverse_patchy_particle_potential(theta: Array | float,
+                                      dist: float,
+                                      abar: float,
+                                      Zc: float,
+                                      Zp: tuple[float, float],
+                                      params: ModelParams,
+                                      lmax: int = 20):
+    dist = dist * params.R
+
+    out0 = ((Zc + Zp[0] + Zp[1]) * (np.exp(params.kappaR) / (1. + params.kappaR)) *
+            (np.exp(-params.kappa * dist) / dist) / (params.epsilon * uc.CONSTANTS.epsilon0)) * params.R ** 2
+
+    for l in range(2, lmax, 2):
+        out1 = (Zp[0] * np.power(abar, l) * eval_legendre(l, np.cos(theta))
+                + Zp[1] * np.power(abar, l) * eval_legendre(l, np.cos(np.pi - theta)))
+        out0 += (((2 * l + 1.) * (kv(l + 0.5, params.kappa * dist) / kv(l + 1.5, params.kappaR))) * out1
+                 / (params.kappa * params.epsilon * uc.CONSTANTS.epsilon0))
+
+    return out0
+
+
 if __name__ == '__main__':
 
     params = ModelParams(R=150, kappaR=3)
     ex1 = expansion.MappedExpansionQuad(0.44, params.kappaR, 0.001, l_max=30)
-    ex2 = expansion.SphericalCap(np.array([[0, 0], [np.pi, 0]]), 0.5, 0.003, l_max=70)
+    # ex2 = expansion.SphericalCap(np.array([[0, 0], [np.pi, 0]]), 0.5, 0.003, l_max=70)
 
     theta = np.linspace(0, np.pi, 1000)
     phi = 0.
     dist = 1
 
+    potential_ic = inverse_patchy_particle_potential(theta, dist, 0.44, -0.002, (0.001, 0.001), params)
+
     potential1 = charged_shell_potential(theta, phi, dist, ex1, params)
-    potential2 = charged_shell_potential(theta, phi, dist, ex2, params)
+    # potential2 = charged_shell_potential(theta, phi, dist, ex2, params)
     # print(potential.shape)
     # print(potential)
 
+    plt.plot(theta, potential_ic)
     plt.plot(theta, potential1.T)
-    plt.plot(theta, potential2.T)
+    # plt.plot(theta, potential2.T)
     plt.show()