from numpy import log10, histogram
from matplotlib.pyplot import figure, show
from matplotlib import gridspec
from Read import ReadRadial, ReadAngleVsEnergy
from Analytic import Analytic_theta
from Constants import degre

def angle_vs_energy(theta,energy,weight,nbBins=50):
   '''
      Compute arrival angle versus energy
      Input: directory name
      Output: energy, arrival angle
   '''

   energy =log10(energy)
   angle,ener = histogram(energy,nbBins,weights=weight*theta)
   nb,   ener = histogram(energy,nbBins,weights=weight)
   ener = 10**ener
   enercenter = (ener[1:nbBins+1]+ener[0:nbBins])/2
   angle /= nb

   return enercenter, angle

def radial(theta,weight,nbBins=50):
   cond = theta!=0
   theta = 2*log10(theta[cond])
   weight = weight[cond]
   dn,angle2 = histogram(theta,nbBins,range=[2*log10(1e-3),2*log10(25)],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 drawRadial(files,all_source_spectrum=False,all_jets=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:
      if all_source_spectrum:
         i=1
         for powerlaw_index in [1,1.5,2,2.5]: 
            theta,dndtheta2 = ReadRadial(fileId,[0,i])
            ax.plot(theta,dndtheta2,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
            i+=1
      elif all_jets:
         i=5
         for jet_opening in ["isotrop","60 deg","30 deg"]: 
            theta,dndtheta2 = ReadRadial(fileId,[0,i])
            ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=jet_opening)
            i+=1
      else:
         theta,dndtheta2 = ReadRadial(fileId,[0,5])
         ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId)

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

def drawAngle_vs_energy(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=False):
   '''
      Plot arrival angle versus energy + analytic expression
      Input: list of directories
      Output: graph of arrival angle versus energy
   '''
   fig = figure()
   ax1 = fig.add_subplot(111)

   for fileId in files:
      if all_source_spectrum:
         i=1
         for powerlaw_index in [1,1.5,2,2.5]: 
            ener,angle = ReadAngleVsEnergy(fileId,[0,i])
            p=ax1.plot(ener,angle,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
      elif all_jets:
         i=5
         for jet_opening in ["isotrop","60 deg","30 deg"]: 
            ener,angle = ReadAngleVsEnergy(fileId,[0,i])
            p=ax1.plot(ener,angle,drawstyle='steps-mid',label=jet_opening)
      else:
         ener,angle = ReadAngleVsEnergy(fileId,[0,1])
         p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId)

      if PlotAnalytic:
         yfit = Analytic_theta(ener,fileId)
         ax1.plot(ener,yfit,'--',color=p[0].get_color(),linewidth=2,label="analytic")

   ax1.legend(loc="best")
   ax1.set_xscale('log')
   ax1.set_yscale('log')
   ax1.set_ylim([0,180])
   ax1.grid(b=True,which='major')
   ax1.set_ylabel("$\\theta$ [deg]")
   ax1.set_xlabel("energy [GeV]")

   show()