2023-11-02

expansion.py

@@ -117,7 +117,7 @@ class MappedExpansionQuad(Expansion):
         coefs = (2 * sigma_m * 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))
         coefs[..., 0] = sigma0 / np.sqrt(4 * np.pi)
-        coefs = rot_sym_expansion(l_array, coefs)
+        coefs = np.squeeze(rot_sym_expansion(l_array, coefs))
         super().__init__(l_array, coefs)
 
 
@@ -300,10 +300,10 @@ if __name__ == '__main__':
     # ex = MappedExpansionQuad(0.44, 3, 1, 10)
     # 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 / 2, 0]]), 0.5, 0.1, 30)
-    # print(ex.coefs)
-    # plot_theta_profile(ex, num=1000, theta_end=2 * np.pi, phi=0)
-    plot_charge_3d(ex)
+    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)
 
     # new_coeffs = expansion_rotation(Quaternion(np.arange(20).reshape(5, 4)).normalized, ex.coeffs, ex.l_array)
     # print(new_coeffs.shape)

functions.py

@@ -22,5 +22,9 @@ def coefficient_Cpm(l, x):
     return x * sps.kv(l + 1 / 2, x) * sps.iv(l + 1 / 2, x)
 
 
+def coefficient_Cim(l, x):
+    return sps.kv(l + 1 / 2, x) / sps.kv(l + 3 / 2, x)
+
+
 def coef_C_diff(l, x):
     return 1 / (x * sps.iv(l + 1 / 2, x) * sps.kv(l + 3 / 2, x))

patch_shape_comparison.py

@@ -0,0 +1,71 @@
+import expansion
+import patch_size
+import potentials
+import numpy as np
+from parameters import ModelParams
+import functions as fn
+import matplotlib.pyplot as plt
+import scipy.special as sps
+
+
+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, 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
+
+
+if __name__ == '__main__':
+
+    target_patch_size = 0.9
+    params = ModelParams(R=150, kappaR=3)
+    sigma_m = 0.001
+
+    def fn1(x):
+        return expansion.MappedExpansionQuad(a_bar=x, kappaR=params.kappaR, sigma_m=sigma_m, l_max=30)
+
+    def fn2(x):
+        return expansion.GaussianCharges(lambda_k=x, omega_k=np.array([[0, 0], [np.pi, 0]]), sigma1=0.001, l_max=30)
+
+    def fn3(x):
+        return expansion.SphericalCap(theta0_k=x, sigma1=0.001, l_max=50, omega_k=np.array([[0, 0], [np.pi, 0]]))
+
+    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_m=sigma_m, l_max=30)
+
+    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)
+
+    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)
+
+    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)
+    # print(potential.shape)
+    # print(potential)
+
+    plt.plot(theta, potential1.T)
+    plt.plot(theta, potential2.T)
+    plt.plot(theta, potential3.T)
+    plt.show()
+
+
+
+

patch_size.py

@@ -1,10 +1,11 @@
 import numpy as np
-from scipy.optimize import bisect
+from scipy.optimize import bisect, root_scalar
 import expansion
 from matplotlib import pyplot as plt
 import parameters
 import potentials
 from pathlib import Path
+from typing import Callable
 
 Expansion = expansion.Expansion
 Array = np.ndarray
@@ -50,23 +51,44 @@ def potential_patch_size(ex: Expansion, params: ModelParams,
     return np.squeeze(np.array(results))
 
 
+def inverse_potential_patch_size(target_patch_size: float,
+                                 ex_generator: Callable[[float], Expansion],
+                                 x0: float,
+                                 params: ModelParams, **ps_kwargs):
+
+    def patch_size_dif(x):
+        ex = ex_generator(x)
+        return potential_patch_size(ex, params, **ps_kwargs) - target_patch_size
+
+    root_result = root_scalar(patch_size_dif, x0=x0)
+
+    if not root_result.converged:
+        raise ValueError('No convergence. Patches of desired size might not be achievable in the given model. '
+                         'Conversely, a common mistake might be target patch size input in degrees.')
+
+    return root_result.root
+
+
 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=kappaR)
-    ex = expansion.MappedExpansionQuad(a_bar=a_bar[:, None], sigma_m=0.001, l_max=30, kappaR=kappaR[None, :])
-
-    # 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()
+    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()
 

potentials.py

@@ -15,7 +15,8 @@ def charged_shell_potential(theta: Array | float,
                             phi: Array | float,
                             dist: float,
                             ex: Expansion,
-                            params: ModelParams) -> Array:
+                            params: ModelParams,
+                            print_idx: int = None) -> Array:
     """
     Electrostatic potential around a charged shell with patches given by expansion over spherical harmonics.
 
@@ -24,6 +25,7 @@ def charged_shell_potential(theta: Array | float,
     :param dist: distance between the particles in units of radius R
     :param ex: Expansion object detailing patch distribution
     :param params: ModelParams object specifying parameter values for the model
+    :param print_idx: if not None, print a single term for debugging purposes
     """
     theta, phi = np.broadcast_arrays(theta, phi)
     angles_shape = theta.shape
@@ -39,9 +41,11 @@ def charged_shell_potential(theta: Array | float,
                  / fn.sph_bessel_k(ex.l_array, params.kappaR))
     l_factors = ex.repeat_over_m(l_factors)
 
-    pot = (1 / (params.kappa * params.epsilon * uc.CONSTANTS.epsilon0)
-            * np.real(np.sum(l_factors[:, None] * ex.coefs[..., None]
-                             * fn.sph_harm(l_array[:, None], m_array[:, None], theta[None, :], phi[None, :]), axis=-2)))
+    factors = l_factors[:, None] * ex.coefs[..., None] * fn.sph_harm(l_array[:, None], m_array[:, None], theta[None, :], phi[None, :])
+    if print_idx is not None:
+        print(l_array[print_idx], m_array[print_idx], np.real(factors[print_idx]))
+
+    pot = (1 / (params.kappa * params.epsilon * uc.CONSTANTS.epsilon0) * np.real(np.sum(factors, axis=-2)))
 
     return pot.reshape(ex.shape + angles_shape)
 
@@ -49,15 +53,18 @@ def charged_shell_potential(theta: Array | float,
 if __name__ == '__main__':
 
     params = ModelParams(R=150, kappaR=3)
-    ex = expansion.MappedExpansionQuad(np.array([0.44, 0.5]), params.kappaR, 0.001, l_max=10)
+    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)
 
     theta = np.linspace(0, np.pi, 1000)
     phi = 0.
     dist = 1
 
-    potential = charged_shell_potential(theta, phi, dist, ex, params)
-    print(potential.shape)
+    potential1 = charged_shell_potential(theta, phi, dist, ex1, params)
+    potential2 = charged_shell_potential(theta, phi, dist, ex2, params)
+    # print(potential.shape)
     # print(potential)
 
-    plt.plot(theta, potential.T)
+    plt.plot(theta, potential1.T)
+    plt.plot(theta, potential2.T)
     plt.show()