# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Plotting the atmospheric neutrino angular distributions
# Data taken from the most recent atmospheric model calculation http://www.icrr.u-tokyo.ac.jp/~mhonda/
#
# We use the year averaged Gran Sasso (with mountain) files at Solar minimum 

import sys
sys.path.append('../src')
from numpy import *
from numpy import random
import matplotlib as mpl
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from Params import *
from NeutrinoFuncs import *
from LabFuncs import *
from PlotFuncs import *
import cmocean

# +
fname = 'AtmNu_GranSasso_SolarMin.d'

#### Load high energy data
Phi_tot,E_high,cosZ,phi_Az = GetAtmNuFluxes(fname)
Phi_high = squeeze(sum(sum(Phi_tot,0),0))

# +
import matplotlib.colors as colors
fig,ax1,ax2,ax3 = MyTriplePlot(size_y=6,wspace=0.05,lfs=31)

Phi_tot1 = sum(Phi_tot[:,:,0:20],2)
Phi_tot2 = sum(Phi_tot[:,:,20:40],2)
Phi_tot3 = sum(Phi_tot[:,:,40:60],2)

maxPhi = amax(Phi_tot1)
minPhi = amax(Phi_tot3)

nlevels = 20

#cmap = cmocean.cm.ice
cmap = cm.YlGnBu_r
norm = colors.BoundaryNorm(boundaries=linspace(0,1,nlevels), ncolors=256)

cnt1 = ax1.contourf(phi_Az,cosZ,Phi_tot1.T/amax(Phi_tot1.T),levels=nlevels,cmap=cmap,vmin=0,vmax=1)
cnt2 = ax2.contourf(phi_Az,cosZ,Phi_tot2.T/amax(Phi_tot2.T),levels=nlevels,cmap=cmap,vmin=0,vmax=1)
im = ax3.pcolormesh(phi_Az,cosZ,Phi_tot3.T/amax(Phi_tot3.T),cmap=cmap,norm=norm,vmin=0,vmax=1)
cnt3 = ax3.contourf(phi_Az,cosZ,Phi_tot3.T/amax(Phi_tot3.T),levels=nlevels,cmap=cmap,vmin=0,vmax=1)

# This is the fix for the white lines between contour levels
for c in cnt1.collections:
    c.set_edgecolor("face")
for c in cnt2.collections:
    c.set_edgecolor("face")
for c in cnt3.collections:
    c.set_edgecolor("face")

ax1.set_title(r'$E_\nu = [10^2,10^3]$ MeV',fontsize=30)
ax2.set_title(r'$E_\nu = [10^3,10^4]$ MeV',fontsize=30)
ax3.set_title(r'$E_\nu = [10^4,10^5]$ MeV',fontsize=30)

ax1.set_yticks([cosZ[-1],cosZ[-1]/2.0,0,cosZ[0]/2.0,cosZ[0]])
ax1.set_yticklabels([-1,-0.5,0,0.5,1])
ax2.set_yticklabels([])
ax3.set_yticklabels([])


phi_ticks = [15,90,180,270,360-15]
phi_ticklabels = [r'$\mathcal{S}$',r'$\mathcal{E}$',r'$\mathcal{N}$',r'$\mathcal{W}$',r'$\mathcal{S}$']
ax1.set_ylabel(r'Zenith angle, $\cos{\theta}$')
ax1.set_xlabel(r'Azimuth, $\phi$')
ax2.set_xlabel(r'Azimuth, $\phi$')
ax3.set_xlabel(r'Azimuth, $\phi$')

ax1.set_xticks(phi_ticks)
ax1.set_xticklabels(phi_ticklabels)
ax1.plot([15,360-15],[0,0],'k--',lw=3)
ax1.text(20,0.03,'Horizon',fontsize=25,color='k')

ax2.set_xticks(phi_ticks)
ax2.set_xticklabels(phi_ticklabels)
ax2.plot([15,360-15],[0,0],'k--',lw=3)
ax2.text(20,0.03,'Horizon',fontsize=25,color='k')

ax3.set_xticks(phi_ticks)
ax3.set_xticklabels(phi_ticklabels)
ax3.plot([15,360-15],[0,0],'k--',lw=3)
ax3.text(20,0.03,'Horizon',fontsize=25,color='k')

cb = cbar(im)
cb.set_ticks([0,0.25,0.5,0.75,1])
cb.ax.set_ylabel(r'$\Phi_{\nu}/\max{(\Phi_\nu)}$',rotation=-90,labelpad=40,fontsize=34)

plt.gcf().text(0.5,1.02,r'{\bf Atmospheric neutrino directions} (Gran Sasso)',fontsize=40,horizontalalignment='center')

plt.show()
MySaveFig(fig,"AtmosphericNeutrinos_costhphi")
# -



StarcoderdataJupyterScriptsDedupFiltered