#!/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")

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}$ [$eV.cm^{-3}$]")

plt.show()