#!/bin/python

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

fig1 = plt.figure()
ax = fig1.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 et 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 et 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 et 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 et Al")

#==== CMB ====
def nCMB(E,z):
   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,z)*hv**2,"-k",label="CMB")

ax.set_xscale('log')
ax.set_yscale('log')
ax.grid(b=True,which='major')
ax.legend(loc="best",frameon=False,framealpha=0.5)
ax.set_xlabel("energy [eV]")
ax.set_ylabel("$n$ [photon.eV.cm$^{-3}$]")

#==== Dominguez vs redshift ==================================================================
fig2 = plt.figure()
ax2 = fig2.add_subplot(111)

z = np.loadtxt("EBL_files/z_Dominguez.dat",unpack=True,usecols=[0])
for z_index in [1,9,12,15,17]:
   lamb,lambdaI = np.loadtxt("EBL_files/lambdaI_Dominguez.dat",unpack=True,usecols=[0,z_index])
   hv = h*c/(lamb*1e-4)  # erg
   density = 1e-6*(4*pi/c)*lambdaI/(hv**2) /(erg_to_GeV*1e9) *(1+z[z_index])**3
   hv = hv *(erg_to_GeV*1e9)
   p=ax2.plot(hv,density*hv**2,"--")
   # CMB
   hv = np.logspace(-4,-1,1000)
   ax2.plot(hv,nCMB(hv,z[z_index])*hv**2,color=p[0].get_color(),linestyle='-',label="z="+str(z[z_index-1]))

ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.grid(b=True,which='major')
ax2.legend(loc="best",frameon=False,framealpha=0.5)
ax2.set_xlabel("energy [eV]")
ax2.set_ylabel("$n$ [photon.eV.cm$^{-3}$]")

plt.show()