import numpy as np
from matplotlib.ticker import NullFormatter
import matplotlib as mpl
import matplotlib.pyplot as plt

# ----------
# Formatting 
# ----------
mpl.rcParams.update({
    'mathtext.fontset'            : 'cm',
    'axes.unicode_minus'          : False,
    'axes.formatter.use_mathtext' : True,
})

plt.rcParams.update({
    'font.size': 14,
    'figure.figsize'    : [8.0, 5.0],
    'font.family'       : "serif",
    'font.serif'        : "cmr10",
    'xtick.major.size'  : 6,
    'xtick.minor.size'  : 3,
    'ytick.major.size'  : 6,
    'ytick.minor.size'  : 3,
})

LabelSize=14
MarkerSize=75
MarkerSize=95
TickSize=20
LegendSize=14
nullfmt=NullFormatter() # no labels

SMALL_SIZE_2 = 19
SMALL_SIZE = 20
MEDIUM_SIZE = 24
BIGGER_SIZE = 24
BIGGEST_SIZE = 30

plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE_2)    # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title

# ----------
# Import data
# ----------
nsosminmaxgapLO = np.loadtxt('cs2_LOgap_D0p140.dat')
nsosminmaxgapNLO = np.loadtxt('cs2_NLOgap_D0p140.dat')
nsosminmax = np.loadtxt('cs2_NLOgap_D0p0.dat')
muB2p6Delta0p136 = np.loadtxt('cs2_NLOgap_D0p140.dat')
muB2p6Delta0p136LO = np.loadtxt('cs2_NLOgap_D0p140.dat')
muB2p6Delta0p0 = np.loadtxt('cs2_NLOgap_D0p0.dat')

# ----------
# Plot
# ----------
fig, ax = plt.subplots()
plt.plot(nsosminmaxgapNLO.transpose()[0],nsosminmaxgapNLO.transpose()[1], color = 'tab:green')
plt.plot(nsosminmaxgapNLO.transpose()[0],nsosminmaxgapNLO.transpose()[2], color = 'tab:green')
plt.fill_between(nsosminmaxgapNLO.transpose()[0], nsosminmaxgapNLO.transpose()[1], nsosminmaxgapNLO.transpose()[2], 
                 color='tab:green', alpha=0.4, label=r'$\Delta^{\!*}_{\rm CFL} = 140 \, {\rm MeV}$ (NLO)')

plt.plot(nsosminmaxgapLO.transpose()[0],nsosminmaxgapLO.transpose()[1], color = 'tab:blue')
plt.plot(nsosminmaxgapLO.transpose()[0],nsosminmaxgapLO.transpose()[2], color = 'tab:blue')
plt.fill_between(nsosminmaxgapLO.transpose()[0], nsosminmaxgapLO.transpose()[1], nsosminmaxgapLO.transpose()[2], 
                 color='tab:blue', alpha=0.4, label=r'$\Delta^{\!*}_{\rm CFL} = 140 \, {\rm MeV}$ (LO)')

plt.plot(nsosminmax.transpose()[0],nsosminmax.transpose()[1], color = 'black', alpha=0.2)
plt.plot(nsosminmax.transpose()[0],nsosminmax.transpose()[2], color = 'black', alpha=0.2)
plt.fill_between(nsosminmax.transpose()[0], nsosminmax.transpose()[1], nsosminmax.transpose()[2], color='black', label=r'NQM', hatch="X", alpha=0.2)
plt.plot(muB2p6Delta0p136.transpose()[0],muB2p6Delta0p136.transpose()[1], color = 'black', label=r'$\mu_{\rm B} = 2.6 \, {\rm GeV}$', alpha=1, linestyle=':')
plt.plot(muB2p6Delta0p0.transpose()[0][0:10],muB2p6Delta0p0.transpose()[1][0:10], color = 'black', alpha=1, linestyle=':')
plt.plot(muB2p6Delta0p136LO.transpose()[0][0:12],muB2p6Delta0p136LO.transpose()[1][0:12], color = 'black', alpha=1, linestyle=':')

plt.axhline(1/3, c='black', ls='--', alpha=0.5, dashes=(5, 3))
plt.xscale('log')
plt.xlim(25,240)
plt.yticks([0.2,0.3,0.4,0.5])
plt.ylim(0.26,.54)
plt.xlabel(r'$n_{\rm B}/n_0$')
plt.ylabel(r'$c_{\rm s}^2$')

ax.set_xticks([30,40,50,60,70, 100, 150,200])
ax.set_xticklabels([30,r"",50,r"",70, 100,r"", 200]) 
plt.legend(loc='upper right')
plt.savefig('speed_of_sound.pdf', format='pdf', bbox_inches='tight')
plt.close()