analytic.py
4.29 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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
from constants import *
from numpy import sqrt, abs, cos, sin, arcsin, searchsorted, loadtxt, log
from scipy.integrate import quad
from read import ReadProfile,resultdir
# Value for analytic expression
Ee_default=1e3 #GeV
Egamma_default=1 #GeV
B_default=1e-17 #Gauss
lambda_B_default=0.3 #Mpc
def Analytic_flux(E_ic):
return me*1e3/4*sqrt(3e9/Ecmb)*1e-9*E_ic**(-3/2.) #GeV
def Analytic_theta_vs_energy(E_ic,fileId):
B,Esource,Dsource=ReadProfile(fileId,[0,1,2]) #Gauss #GeV #Mpc
Esource=1e5
#delta_to_theta=lambda_gg(Esource)/Dsource
delta_to_theta=1.9440974051313842/Dsource
RL0=RL(Esource/2,B)
Dic0=Dic(Esource/2)
delta=Dic0/(2*RL0)*((Esource/(2*Ee(E_ic)))**2 -1)
return abs(arcsin(delta_to_theta*sin(delta)))*degre
def Analytic_delay_vs_theta(theta,fileId):
Esource,distSource=ReadProfile(fileId,[1,2]) #GeV #Mpc
lgg = lambda_gg(Esource)[0]
E0_source = Esource*1e-3 #TeV
delta_ic = arcsin(distSource/lgg*sin(theta))
Dic0=Dic(E0_source/2)
c_delta_t = (lgg*(1-cos(delta_ic)) - distSource*(1-cos(theta)))*log(lgg/theta)
return c_delta_t *Mpc/c # sec.
def Analytic_delay_vs_energy(Egamma, fileId):
B,Esource,Dsource=ReadProfile(fileId,[0,1,2]) #Gauss #GeV #Mpc
RL0=RL(Esource/2,B)
Dic0=Dic(Esource/2)
E_e = Ee(Egamma)
lgg = lambda_gg(Esource)[0]
#lgg=1.9440974051313842
delta_ic = Dic0/(2*RL0)*((Esource/2/E_e)**2 -1)
theta = arcsin(lgg/Dsource*sin(delta_ic))
c_delta_t = (lgg*(1-cos(delta_ic)) - Dsource*(1-cos(theta)))*log(lgg/theta)
return c_delta_t *Mpc/c # sec.
# Compton accumulation
def ECompton_threshold(Compton_threshold = 0.005,z=0):
return Compton_threshold/(4/3*Ecmb*(1+z)/me*1e-3) *me*1e-6 #GeV
# Compton scattering
def Dic(Ee=Ee_default,z=0): # Ee (GeV)
return 3*(me*1e3)**2/(4*sigmaT*Ee*1e9*rhoCMB*(1+z)**4) /Mpc #Mpc
def lambdaIC(z=0):
return 1/(sigmaT*Mpc*nCMB*(1+z)**3) #Mpc
def Eic(Ee=Ee_default,z=0):
return 4*Ecmb*(1+z) *Ee**2/(3*me**2)*1e3 #GeV
def Ee(Egamma=Egamma_default,z=0):
return me*sqrt((3*Egamma*1e-3 )/(4*Ecmb*(1+z))) #GeV
def tIC():
return lambdaIC()/(c*yr/Mpc) #yr
# Larmor radius
def RL(Ee=Ee_default,B=B_default):
return (Ee/erg_to_GeV)/(e*B) /Mpc #Mpc
# Magnetic deflection
def delta(Ee=Ee_default,B=B_default):
return lambdaIC()/RL(Ee,B)*degre
def Delta(Ee=Ee_default,B=B_default,lambda_B=lambda_B_default):
Delta1=Dic(Ee)/RL(Ee,B)*degre #degre if Dic >> lambda_B
Delta2=sqrt(Dic(Ee)*lambda_B)/RL(Ee,B)*degre #degre if Dic << lambda_B
return Delta1, Delta2
# Threshold energy when Dic = RL
def Ethreshold_ic(Ee=Ee_default,B=B_default):
return Eic(Ee)*Dic(Ee)/RL(Ee,B) # GeV
# Pair production
def Ethreshold_gg():
return (me)**2/Eebl *1e-3 #GeV
def lambda_gg(Egamma=1e3,z=0): # Egamma (GeV)
# Bilinear interpolation
z_tab = [0,0.5,1,2,3]
E_tab = loadtxt("lambda_e.dat",unpack=True,usecols=[0])*me*1e-6
i2 = searchsorted(z_tab,z)
if z<=0:
fy=1
i1=i2
elif z>3:
fy=1
i2-=1
i1=i2
else:
i1 = i2-1
fy=(z-z_tab[i1])/(z_tab[i2]-z_tab[i1])
j2 = searchsorted(E_tab,Egamma)
j1 = j2-1
fx=(Egamma-E_tab[j1])/(E_tab[j2]-E_tab[j1])
lambda_e = loadtxt("lambda_e.dat",unpack=True,usecols=[i1+1,i2+1,i1+6,i2+6])
lambda11 = lambda_e[0,j1]*((1-fx)*(1-fy))
lambda12 = lambda_e[1,j1]*((1-fx)*fy)
lambda21 = lambda_e[0,j2]*(fx*(1-fy))
lambda22 = lambda_e[1,j2]*(fx*fy)
lambda_proper = lambda11 + lambda12 + lambda21 + lambda22
lambda11 = lambda_e[2,j1]*((1-fx)*(1-fy))
lambda12 = lambda_e[3,j1]*((1-fx)*fy)
lambda21 = lambda_e[2,j2]*(fx*(1-fy))
lambda22 = lambda_e[3,j2]*(fx*fy)
lambda_comobile = lambda11 + lambda12 + lambda21 + lambda22
return lambda_proper, lambda_comobile, 800.e3/Egamma #Mpc (from Durrer and Neronov 2013)
def comobileTime(z):
return -1/(H0*(1+z)*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
def distPhoton(z):
return -c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
def distLepton(z,E):
beta = sqrt(1-m**2*c**4/E**2)
return -beta*c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
def properIntegrand(z):
return -c/(H0*(1+z)*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
def comobileIntegrand(z):
return -c/(H0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
def distance(z=0):
return quad(properIntegrand,z,0)[0]/Mpc, quad(comobileIntegrand,z,0)[0]/Mpc