#!/bin/python

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

from Modules.Constants import *

fig = plt.figure()
ax = fig.add_subplot(111)

z=2

# ==== Dominguez ====
lamb,lambdaI = np.loadtxt("EBL_files/lambdaI_Dominguez.dat",unpack=True,usecols=[0,15])

#lamb,lambdaI =np.loadtxt("EBL_files/lambdaI_Dominguez.dat",unpack=True,usecols=[0,1])

hv = h*c/(lamb*1e-4)  # erg
density = 1e-6*(4*pi/c)*lambdaI/(hv**2) /(erg_to_GeV*1e9) *(1+z)**3
hv = hv *(erg_to_GeV*1e9)

ax.plot(hv,density*hv**2,"--b",label="Dominguez and Al")

# ==== Kneiske and Doll - "best fit" ====
hv,density = np.loadtxt("EBL_files/n_bestfit10.dat",unpack=True,usecols=[0,1])

ax.plot(hv,density*hv**2,"--r",label="Kneiske and Doll - 'best fit'")

# ==== Kneiske and Doll - "lower limit" ====

hv,density = np.loadtxt("EBL_files/n_lowerlimit10.dat",unpack=True,usecols=[0,179])
#hv,density =np.loadtxt("EBL_files/n_lowerlimit10.dat",unpack=True,usecols=[0,1])

density=density*(1+z)**3

ax.plot(hv,density*hv**2,"--g",label="Kneiske and Doll - 'lower limit'")

# ==== Fransceschini ====

hv,density = np.loadtxt("EBL_files/n_Fra.dat",unpack=True,usecols=[0,11])
#hv,density = np.loadtxt("EBL_files/n_Fra.dat",unpack=True,usecols=[0,1])

density=density*(1+z)**3

ax.plot(hv,density*hv**2,"--c",label="Fraceschini")

# ==== Finke ====

hv,density = np.loadtxt("EBL_files/n_Finke.dat",unpack=True,usecols=[0,201])
#hv,density = np.loadtxt("EBL_files/n_Finke.dat",unpack=True,usecols=[0,1])

density=density*(1+z)**3

ax.plot(hv,density*hv**2,"--m",label="Finke and Al")

# ==== Gilmore ====

hv,density = np.loadtxt("EBL_files/n_Gil.dat",unpack=True,usecols=[0,14])
#hv,density = np.loadtxt("EBL_files/n_Gil.dat",unpack=True,usecols=[0,1])

ax.plot(hv,density*hv**2,"--y",label="Gilmore and Al")

#==== CMB ====
def nCMB(E):
   kTcmb = k*Tcmb*erg_to_GeV*1e9*(1+z)
   theta = E/kTcmb
   nCMB=(hb*c*erg_to_GeV*1e9)**(-3) *(E/np.pi)**2 /(np.exp(theta)-1)    
   return nCMB

hv = np.logspace(-4,-1,1000)
ax.plot(hv,nCMB(hv)*hv**2,"-k",label="CMB")

ax.set_xscale('log')
ax.set_yscale('log')
ax.grid(b=True,which='major')
ax.legend(loc="best",title="z = %.0f"%z)
ax.set_xlabel("energy [eV]")

ax.set_ylabel("$n_{EBL}$ [photon.eV.cm$^{-3}$]")

plt.show()