charged-shells/analysis/peak_heigth.py

import matplotlib.pyplot as plt
from charged_shells import expansion, interactions, mapping
from charged_shells.parameters import ModelParams
import numpy as np
from typing import Literal
from pathlib import Path
from functools import partial
import json

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)

    ex1 = expansion.MappedExpansionQuad(a_bar[:, None], params.kappaR[None, :], 0.001, l_max=l_max)
    ex2 = ex1.clone()
    if which == 'ep':
        ex1.rotate_euler(alpha=0, beta=np.pi / 2, gamma=0)

    energy_fn = mapping.parameter_map_two_expansions(partial(interactions.charged_shell_energy, dist=dist), 1)
    energy = energy_fn(ex1, ex2, params)

    fig, ax = plt.subplots()
    for en, ab in zip(energy, a_bar):
        ax.plot(kappaR, en / en[0], label=rf'$\bar a = {ab:.1f}$')
        # ax.plot(kappaR, en, label=rf'$\bar a = {ab:.1f}$')
    ax.legend(fontsize=12)
    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
    ax.set_xlabel(r'$\kappa R$', fontsize=15)
    ax.set_ylabel(rf'$\bar V$', fontsize=15)
    plt.tight_layout()
    if save_as is not None:
        plt.savefig(save_as, dpi=600)
    plt.show()


def IC_peak_energy_plot(config_data: dict,
                        a_bar: list,
        which: Literal['ep', 'pp'],
                     save_as: Path = None):

    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_11"))
    em_data = np.load(em_data_path.joinpath("pair_energy_fig11.npz"))
    data = em_data['fixA']

    if which == 'ep':
        column_idx = 4
    elif which == 'pp':
        column_idx = 3
    else:
        raise ValueError

    abar, inverse, counts = np.unique(data[:, 1], return_counts=True, return_inverse=True)
    # print(indices, counts)
    print(inverse)
    # sort_abar = np.argsort(indices)
    # print(data[:, 1][indices])
    # print(sort_abar)
    # print(data[:, 1][indices[sort_abar]])

    # energies = data[:, column_idx].reshape(-1, counts[0])
    # kappaR = data[:, 2].reshape(-1, counts[0])
    # print(len(kappaR), len(energies), len(abar))

    fig, ax = plt.subplots()
    for i in range(len(abar)):
        if abar[i] in a_bar:
            idx, = np.nonzero(inverse == i)
            kR = data[idx, 2]
            en = data[idx, column_idx]
            sort = np.argsort(kR)
            # ax.plot(kR[sort], en[sort] / en[sort][0], label=rf'$\bar a = {abar[i]:.2f}$')
            ax.plot(kR[sort], en[sort], label=rf'$\bar a = {abar[i]:.2f}$')
    ax.legend(fontsize=12)
    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=15)
    ax.set_xlabel(r'$\kappa R$', fontsize=15)
    ax.set_ylabel(rf'$\bar V_{{{which}}}$', fontsize=15)
    plt.tight_layout()
    ax.set_yscale('log')
    ax.set_xscale('log')
    if save_as is not None:
        plt.savefig(save_as, dpi=600)
    plt.show()


if __name__ == '__main__':
    
    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
        config_data = json.load(config_file)

    kappaR = np.arange(0.5, 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_ep.pdf')
    #                  )
    
    IC_peak_energy_plot(config_data, a_bar=[0.2, 0.24, 0.36, 0.52, 0.8], which='pp',
                     # save_as=Path('/home/andraz/ChargedShells/Figures/Emanuele_data/nonmonotonicity_check_ep.pdf')
                     )