charged-shells/patch_size.py

import numpy as np
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
ModelParams = parameters.ModelParams


@expansion.map_over_expansion
def charge_patch_size(ex: Expansion, phi: float = 0, theta0: Array | float = 0, theta1: Array | float = np.pi / 2):
    return bisect(lambda theta: ex.charge_value(theta, phi), theta0, theta1)


def potential_patch_size(ex: Expansion, params: ModelParams,
                         phi: float = 0, theta0: Array | float = 0, theta1: Array | float = np.pi / 2,
                         match_expansion_axis_to_params: int = None):

    # this is more complicate to map over leading axes of the expansion as potential also depends on model parameters,
    # with some, such as kappaR, also being the parameters of the expansion in the first place. When mapping,
    # we must therefore provide the expansion axis that should match the collection of parameters in params.

    meap = match_expansion_axis_to_params  # just a shorter variable name

    @expansion.map_over_expansion
    def potential_zero(exp: Expansion, prms: ModelParams):
        return bisect(lambda theta: potentials.charged_shell_potential(theta, phi, 1, exp, prms), theta0, theta1)

    params_list = params.unravel()
    if meap is not None:
        expansion_list = [Expansion(ex.l_array, np.take(ex.coefs, i, axis=meap)) for i in range(ex.shape[meap])]
    else:
        expansion_list = [ex for _ in params_list]

    if not len(expansion_list) == len(params_list):
        raise ValueError(f'Axis of expansion that is supposed to match params does not have the same length, got '
                         f'len(params.unravel()) = {len(params_list)} and '
                         f'expansion.shape[{meap}] = {len(expansion_list)}')

    results = []
    for exp, prms in zip(expansion_list, params_list):
        results.append(potential_zero(exp, prms))

    if meap is not None:
        return np.array(results).swapaxes(0, meap)  # return the params-matched axis to where it belongs
    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, 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]]))
    # 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)

    markers = ['o', 's', '^', 'D', 'v', '<', '>', 'p', '*', 'H', '+', 'x']

    fig, ax = plt.subplots()
    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_0$', fontsize=13)
    plt.legend(fontsize=12)
    # plt.savefig(Path("/home/andraz/ChargedShells/Figures/patch_size_potential.png"), dpi=600)
    plt.show()