from numpy import log10, histogram, arange
from matplotlib.pyplot import figure, show
from matplotlib import gridspec
from read import ReadRadial, ReadAngleVsEnergy, ReadGeneration
from constants import degre

def arrivalAngle(theta,weight,nbBins=50,theta_range=[1e-6,180]):
   cond = theta>0
   theta = log10(theta[cond])
   weight = weight[cond]
   dn,angle2 = histogram(theta,nbBins,range=log10(theta_range),weights=weight)
   angle2=10**angle2
   #theta = theta**2
   #dn,angle2 = histogram(theta,nbBins,range=[0,25],weights=weight)
   thetacenter = (angle2[1:nbBins+1]+angle2[0:nbBins])/2
   dtheta = angle2[1:nbBins+1]-angle2[0:nbBins]
   dndtheta = dn/dtheta
   dndtheta /= max(dndtheta)

   return thetacenter, dndtheta
 
def drawArrivalAngle(files,plot_gen=False,plot_energy_range=False,plot_time_range=False):
   '''
      Plot flux versus arrival angle
      Input: list of directories
      Output: graph of flux versus angle
   '''
   fig = figure()
   ax = fig.add_subplot(111)

   for fileId in files:
      theta,dndtheta = ReadRadial(fileId,[0,1])
      ax.plot(theta,theta*dndtheta,drawstyle='steps-mid',label=fileId)

      if plot_energy_range:
         theta,dNdtheta1,dNdtheta2,dNdtheta3 = ReadRadial(fileId,[0,2,3,4])
         ax.plot(theta,theta*dNdtheta1,drawstyle='steps-mid',linestyle='--',label="MeV band")
         ax.plot(theta,theta*dNdtheta2,drawstyle='steps-mid',linestyle='--',label="GeV band")
         ax.plot(theta,theta*dNdtheta3,drawstyle='steps-mid',linestyle='--',label="TeV band")

      if plot_time_range:
         theta,dNdtheta1,dNdtheta2,dNdtheta3 = ReadRadial(fileId,[0,5,6,7])
         ax.plot(theta,theta*dNdtheta1,drawstyle='steps-mid',linestyle='--',label="$t_{obs}$=30days")
         ax.plot(theta,theta*dNdtheta2,drawstyle='steps-mid',linestyle='--',label="$t_{obs}$=1yr")
         ax.plot(theta,theta*dNdtheta3,drawstyle='steps-mid',linestyle='--',label="$t_{obs}$=5yrs")

      if plot_gen:
         theta,dNdtheta1,dNdtheta2,dNdtheta3 = ReadRadial(fileId,[0,8,9,10])
         ax.plot(theta,theta*dNdtheta1,drawstyle='steps-mid',linestyle='--',label="gen=1")
         ax.plot(theta,theta*dNdtheta2,drawstyle='steps-mid',linestyle='--',label="gen=2")
         ax.plot(theta,theta*dNdtheta3,drawstyle='steps-mid',linestyle='--',label="gen=3")
         #i=8
         #for gen in ReadGeneration(fileId,[0]):
         #   theta,dndtheta2 = ReadRadial(fileId,[0,i])
         #   ax.plot(theta,dndtheta2,drawstyle='steps-mid',linestyle='--',label="gen=%.0f"%gen)
         #   i+=1

   ax.legend(loc="best")
   ax.set_xscale('log')
   ax.set_yscale('log')
   ax.grid(b=True,which='major')
   ax.set_xlabel("$\\theta$ [deg]")
   ax.set_ylabel("$\\theta*dN/d\\theta$ [arbitrary]")
   
   show()