EBL-spectrums.py
1.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#!/bin/python
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
from Constantes import *
fig = plt.figure()
ax = fig.add_subplot(111)
z=0
# ==== 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,"--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,"--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,"--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,"--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,"--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,"--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}$ [$cm^{-3} eV^{-1}$]")
plt.show()