charged-shells/analysis/energy_gap.py

import matplotlib.pyplot as plt
import numpy as np
from charged_shells import expansion, interactions, mapping
from charged_shells.parameters import ModelParams
from functools import partial
from pathlib import Path
from matplotlib import cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.colors import TwoSlopeNorm
from matplotlib.ticker import FuncFormatter
import json

Expansion = expansion.Expansion


def energy_gap(ex1: Expansion, params: ModelParams, dist=2., match_expansion_axis_to_params=None):
    ex2 = ex1.clone()
    ex3 = ex1.clone()

    ex2.rotate_euler(alpha=0, beta=np.pi/2, gamma=0)  # to get EP config between ex1 and ex2

    energy_fn = mapping.parameter_map_two_expansions(partial(interactions.charged_shell_energy, dist=dist),
                                                     match_expansion_axis_to_params)
    energy_ep = energy_fn(ex1, ex2, params)
    energy_pp = energy_fn(ex1, ex3, params)

    return (energy_pp - energy_ep) / energy_pp


def abar_kappaR_dependence(save_as=None):

    kappaR = np.linspace(0.01, 25, 25)
    a_bar = np.array([0.2, 0.4, 0.6, 0.8])
    params = ModelParams(R=150, kappaR=kappaR)
    ex = expansion.MappedExpansionQuad(a_bar[:, None], kappaR[None, :], 0.001)
    # ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
    gap = energy_gap(ex, params, match_expansion_axis_to_params=1)

    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
    for g, lbl in zip(gap, [rf'$\bar a={a}$' for a in a_bar]):
        ax.plot(kappaR, g, label=lbl)
    ax.legend(fontsize=17)
    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
    ax.set_xlabel(r'$\kappa R$', fontsize=15)
    ax.set_ylabel(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', fontsize=20)
    plt.tight_layout()
    if save_as is not None:
        plt.savefig(save_as, dpi=600)
    plt.show()


def abar_kappaR_dependence2(save_as=None):

    kappaR = np.array([1, 3, 10, 30])
    a_bar = np.linspace(0.2, 0.8, 30)
    params = ModelParams(R=150, kappaR=kappaR)
    ex = expansion.MappedExpansionQuad(a_bar[:, None], kappaR[None, :], 0.001)
    # ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
    gap = energy_gap(ex, params, match_expansion_axis_to_params=1)

    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
    for g, lbl in zip(gap.T, [rf'$\kappa R={kR}$' for kR in kappaR]):
        ax.plot(a_bar, g, label=lbl)
    ax.legend(fontsize=17)
    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
    ax.set_xlabel(r'$\bar a$', fontsize=15)
    ax.set_ylabel(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', fontsize=20)
    plt.tight_layout()
    if save_as is not None:
        plt.savefig(save_as, dpi=600)
    plt.show()


def charge_kappaR_dependence(a_bar, min_charge, max_charge, save_as=None, cmap=cm.jet):

    kappaR = np.linspace(0.01, 10, 50)
    params = ModelParams(R=150, kappaR=kappaR)
    charge = np.linspace(min_charge, max_charge, 100)
    ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001, sigma0=charge)
    # ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
    gap = energy_gap(ex, params, match_expansion_axis_to_params=0)

    colors = cmap(np.linspace(0, 1, len(charge)))
    sm = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=np.min(charge), vmax=np.max(charge)))
    sm.set_array([])

    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
    for g, lbl, c in zip(gap.T, [rf'$\sigma_0 ={1000 * c:.2f}$' for c in charge], colors):
        ax.plot(kappaR, g, label=lbl, c=c)
    # ax.legend(fontsize=17)
    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
    ax.set_xlabel(r'$\kappa R$', fontsize=15)
    ax.set_ylabel(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', fontsize=20)

    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='5%', pad=0.05)
    cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
    cbar.set_label(r'$\tilde \sigma_0$', rotation=90, labelpad=15, fontsize=12)

    plt.tight_layout()
    if save_as is not None:
        plt.savefig(save_as, dpi=600)
    plt.show()


def charge_kappaR_dependence_heatmap(a_bar, min_charge, max_charge, save_as=None, cmap=cm.jet):

    kappaR = np.linspace(0.01, 10, 50)
    params = ModelParams(R=150, kappaR=kappaR)
    charge = np.linspace(min_charge, max_charge, 100)
    ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001, sigma0=charge)
    # ex = expansion.MappedExpansionQuad(a_bar, kappaR, 0.001)
    gap = energy_gap(ex, params, match_expansion_axis_to_params=0)

    norm = TwoSlopeNorm(vmin=np.min(gap), vcenter=0, vmax=np.max(gap))
    sm = cm.ScalarMappable(cmap=cmap, norm=norm)
    sm.set_array([])

    def y_formatter(x, pos):
        return f"{charge[int(x)-1]:.2f}"

    def x_formatter(x, pos):
        return f"{kappaR[int(x)-1]:.2f}"

    fig, ax = plt.subplots(figsize=plt.figaspect(1))
    ax.imshow(gap.T, cmap=cmap, origin='lower',
              # extent=[kappaR.min(), kappaR.max(), charge.min(), charge.max()]
              )
    # ax.legend(fontsize=17)
    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
    ax.set_xlabel(r'$\kappa R$', fontsize=15)
    ax.set_ylabel(r'$\tilde \sigma_0$', fontsize=15)

    plt.gca().xaxis.set_major_formatter(FuncFormatter(x_formatter))
    plt.gca().yaxis.set_major_formatter(FuncFormatter(y_formatter))

    # ax.set_xticks(kappaR)
    # ax.set_yticks(charge)

    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='5%', pad=0.05)
    cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
    cbar.set_label(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', rotation=90, labelpad=20, fontsize=12)

    plt.tight_layout()
    if save_as is not None:
        plt.savefig(save_as, dpi=600)
    plt.show()


def IC_gap_plot(config_data: dict, save_as=None):
    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
    em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
    for k in list(em_data.keys()):
        data = em_data[k]
        print(k, data.shape)
        for i in range(3):
            print(np.unique(data[:, i]))
        print('\n')


def IC_gap_kappaR(config_data: dict, save_as=None):
    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
    em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
    data = em_data['fixA']

    print(data)

    sort = np.argsort(data[:, 2])
    xdata = data[:, 2][sort]
    ydata = data[:, 3][sort]

    plt.plot(xdata, ydata)
    plt.xlabel('kappaR')
    plt.ylabel('gap')
    plt.show()


def IC_gap_abar(config_data: dict, save_as=None):
    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
    em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
    data = em_data['fixM']

    print(data)

    sort = np.argsort(data[:, 1])
    xdata = data[:, 1][sort]
    ydata = data[:, 3][sort]

    plt.plot(xdata, ydata)
    plt.xlabel('abar')
    plt.ylabel('gap')
    plt.show()


def IC_gap_charge_at_abar(a_bar, config_data: dict, save_as=None, cmap=cm.coolwarm, which_change='changezp'):
    em_data_path = (Path(config_data["emanuele_data"]).joinpath("FIG_10"))
    em_data = np.load(em_data_path.joinpath("relative_gap.npz"))
    data = em_data[which_change]

    relevant_indices = data[:, 1] == a_bar
    if not np.any(relevant_indices):
        raise ValueError(f'No results for given a_bar = {a_bar}. Possible values: {np.unique(data[:, 1])}')

    data = data[relevant_indices]
    print('Unique gap sizes', np.unique(data[:, 3]))

    charge, inverse, counts = np.unique(data[:, 0], return_counts=True, return_inverse=True)

    colors = cmap(np.linspace(0, 1, len(charge)))
    sm = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=np.min(charge), vmax=np.max(charge)))
    sm.set_array([])

    fig, ax = plt.subplots(figsize=plt.figaspect(0.5))
    for i, c in enumerate(colors):
        idx, = np.nonzero(inverse == i)
        kR = data[idx, 2]
        gap = data[idx, 3]
        sort = np.argsort(kR)
        ax.plot(kR[sort], gap[sort], c=c)
    # ax.legend(fontsize=17)
    ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12)
    ax.set_xlabel(r'$\kappa R$', fontsize=15)
    ax.set_ylabel(r'$\frac{V_{pp}-V_{ep}}{V_{pp}}$', fontsize=20)
    ax.set_xlim(-0.25, 10.25)
    # ax.set_ylim(-4.25, 6.25)

    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='5%', pad=0.05)
    cbar = fig.colorbar(sm, cax=cax, orientation='vertical')
    cbar.set_label(r'$\tilde \sigma_0$', rotation=90, labelpad=15, fontsize=12)

    plt.tight_layout()
    if save_as is not None:
        plt.savefig(save_as, dpi=600)
    plt.show()


def main():
    with open(Path("/home/andraz/ChargedShells/charged-shells/config.json")) as config_file:
        config_data = json.load(config_file)

    # abar_kappaR_dependence(Path("/home/andraz/ChargedShells/Figures/full_amplitude_kappaR_dep.png"))
    # abar_kappaR_dependence2(Path("/home/andraz/ChargedShells/Figures/full_amplitude_abar_dep.png"))

    # charge_kappaR_dependence(a_bar=0.5, min_charge=-0.003, max_charge=0.003,
    #                          save_as=Path("/home/andraz/ChargedShells/Figures/full_amplitude_charge_abar05.png"),
    #                          cmap=cm.coolwarm)

    # charge_kappaR_dependence_heatmap(a_bar=0.5, min_charge=-0.003, max_charge=0.003,
    #                          save_as=Path("/home/andraz/ChargedShells/Figures/full_amplitude_heatmap_abar05.png"),
    #                          cmap=cm.bwr)

    # IC_gap_plot(config_data)

    # IC_gap_kappaR(config_data)
    # IC_gap_abar(config_data)

    IC_gap_charge_at_abar(0.8, config_data, which_change='changezc',
                          # save_as=Path("/home/andraz/ChargedShells/Figures/Emanuele_data/IC_full_amplitude_charge_abar08.png")
                          )


if __name__ == '__main__':

    main()