EBL-spectrums.py
3.05 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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#!/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()