charged-shells/patch_shape_comparison.py

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()