Analytic.py
1.78 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
from Constants import *
from numpy import sqrt, abs, sin, arcsin
from Read import ReadE_source, ReadD_source, ReadEGMF
# Value for analytic expression
Ee_default=1e3 #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**(1/2.) #GeV
def Analytic_theta(E_ic,fileId):
Esource=ReadE_source(fileId) #GeV
Dsource=ReadD_source(fileId) #Mpc
B=ReadEGMF(fileId) #Gauss
delta_to_theta=lambda_gg(Esource)/Dsource
RL0=RL(Esource/2,B)
Dic0=Dic(Esource/2)
Ee2= E_ic /Eic(1) #GeV^2
delta=Dic0/(2*RL0)*((Esource/2)**2 /Ee2 -1)
return abs(arcsin(delta_to_theta*sin(delta))*degre)
# Compton accumulation
def ECompton_threshold(Compton_threshold = 0.05):
return Compton_threshold/(4/3*2.7*Tcmbp) *me*1e-6 #GeV
# Compton scattering
def Dic(Ee=Ee_default): # Ee (GeV)
return 3*(me*1e-6)**2/(4*sigmaT*rhoCMB*1e-9*Ee) /Mpc #Mpc
def lambdaIC():
return 1/(nCMB*sigmaT*Mpc) #Mpc
def Eic(Ee=Ee_default):
return 4*Ecmb*Ee**2/(3*me**2)*1e3 #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): # Egamma (GeV)
return 0.8 /(Egamma*1e-3) *1e3#Mpc (from Durrer and Neronov 2013)