-
Notifications
You must be signed in to change notification settings - Fork 0
/
SolNuFluxoptimized.py
128 lines (96 loc) · 3.94 KB
/
SolNuFluxoptimized.py
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
from __future__ import division
import numpy as np
import scipy.integrate as integrate
from scipy.integrate import solve_ivp
from scipy.interpolate import UnivariateSpline
import matplotlib.pyplot as plt
Es_npts = 10
log10_Es_min = 2.0
log10_Es_max = 10.0
Es=np.power(10., np.linspace(log10_Es_min, log10_Es_max, Es_npts))#GeV
Es=np.append(Es, [5e3, 4.5e4, 5e5, 5e7 ])
Es=np.sort(Es)
c=299792458*100 # cm/s
def L0(E, z, k=1, gamma=2, E_max=1.0e7):
return k*np.power(E*(1+z),-gamma)*np.exp(-E*(1+z)/E_max*(1+z))
def W(z, a=3.4 , b=-0.3 , c1=-3.5 , B=5000 , C=9 , eta=-10):
return ((1+z)**(a*eta)+((1+z)/B)**(b*eta)+((1+z)/C)**(c1*eta))**(1/eta)
def L(z,E):
return W(z)*L0(E, z)
def H(z, H0=70*1e5/3.085677581e24, OM=0.3089, OL=0.6911):
return H0*np.sqrt(OM*(1.+z)**3. + OL)
def sigma(E, g, M, m): # cm^2
return (g**4/(4*np.pi))*(2*E*m)/((2*E*m-M**2)**2+((M**4*g**4)/(16*np.pi**2)))* 0.389376e-27
def dsigma(Ep, E, g, M, m):
if Ep < E:
return 0
else:
return sigma(Ep, g, M, m)*3/Ep*((E/Ep)**2+(1-(E/Ep))**2)
def A(z, E, n, g, M, m): # cm^-3*GeV^-1*s^-1
return c*56*(1+z)**3*sigma(E*(1+z), g, M, m)*n
def fun(z, E, n, g, M, m):
return (-A(z, E, n, g, M, m)+(1+z)*L(z,E))/(-(1+z)*H(z))
def First_guess(E, g, M, m):
sol=np.zeros(E.size)
for i in np.arange(E.size):
sol_return=solve_ivp(fun=lambda z, n: fun(z, E[i], n, g, M, m),t_span=[4.0, 0.0], y0=[0.0], t_eval=[0.0],
rtol=1e-8, atol=1e-8, method='LSODA') #method='RK45'
sol[i] = sol_return.y[-1,-1]
return sol
def Re(E, z, g, M, m, interp):
def integrand(x, E, z, g, M, m):
return dsigma(x, E*(1+z), g, M, m) *interp(x)
def integrator(E, z, g, M, m):
return integrate.quad(lambda x: integrand(x, E, z,g, M, m), E*(1+z), 1e9, epsabs=1e-2, epsrel=1e-2)[0]
return (1+z)*c*56*(1+z)**3*integrator(E, z, g, M ,m)
def fun_full(z, E, n, g, M, m, interp):
return (-A(z, E, n, g, M, m)+(1+z)*L(z,E)+Re(E, z, g, M, m, interp))/(-(1+z)*H(z))
def final_sol(E, g, M, m):
F1=First_guess(E, g, M, m)
interpF1_linear = UnivariateSpline(E, F1, k=1, ext=0)
sol_array=np.zeros(E.size)
count = 0
while 1:
for i in np.arange(E.size):
sol_return=solve_ivp(fun=lambda z, n: fun_full(z, E[i], n, g, M, m, interpF1_linear),t_span=[4.0, 0.0], y0=[0.0], t_eval=[0.0],
rtol=1e-8, atol=1e-8, method='LSODA') #method='RK45'
sol_array[i]=sol_return.y[-1,-1]
error = sum(np.subtract(sol_array,F1)**2)
print(error)
F1 = np.copy(sol_array)
interpF1_linear = UnivariateSpline(E, F1, k=1, ext=0)
count = count + 1
if error < 1 or count > 10 :
break
return sol_array
Final = final_sol(Es, 0.03, 0.01, 1e-10)
#log10_E_test_min = 2.0
#log10_E_test_max = 10.0
#E_test_npts = 10
#E_test=np.power(10., np.linspace(log10_E_test_min, log10_E_test_max, E_test_npts))
#E_test=np.append(E_test, [5e3, 4.5e4, 5e5, 5e7 ])
#E_test=np.sort(E_test)
#
#flux = E_test*E_test * final_sol(E_test, 0.03, 0.01, 1e-10)
#renorm_flux = max(flux)
#flux = flux/renorm_flux
#
#plt.rcParams['xtick.labelsize']=26
#plt.rcParams['ytick.labelsize']=26
#plt.rcParams['legend.fontsize']=18
#plt.rcParams['legend.borderpad']=0.4
#plt.rcParams['axes.labelpad']=10
#plt.rcParams['ps.fonttype']=42
#plt.rcParams['pdf.fonttype']=42
#
#fig = plt.figure(figsize=[9,9])
#ax = fig.add_subplot(1,1,1)
#
#ax.plot(E_test, flux, label="Flux")
#ax.set_xlabel(r'Neutrino energy $E$ [GeV]', fontsize=25)
#ax.set_ylabel(r'Neutrino flux [$10^{-8}$ GeV cm$^{-2}$ s$^{-1}$ sr$^{-1}$]', fontsize=25)
#plt.xscale('log')
#ax.set_xlim((10.**3., 10.**8.))
#ax.set_ylim((0.0, 1.0))
#plt.legend(loc='upper right')
#plt.savefig('First_sol_flux.png', bbox_inches='tight', dpi=300)