ich habe (mal wieder) eine Frage und hoffe auf eure Hilfe, da ich wirklich nicht mehr weiter weiß.
Ich bin gerade dabei eine Funktion zu integrieren, die an beiden Grenzen eine Singularität beinhalten.
Da es einfacher ist, wenn ich mein Code poste hier erstmal meine Idee
Code: Alles auswählen
from scipy import integrate
from pylab import *
#Alle Energien sind in TeV. Laengeneinheit in cm und Zeit in s. Teilchendichte in cm^3
c=3e10 #cm/s
d=1.2e25 #Abstand zur Quelle
R=1e15 #Radius des Plasmoids
q=1e-4 #injection
t=1e2*q*R/(28*math.pi*d**2) #Vorfaktor
a=np.pi**2.0/700.0
E0=1e8
n=1 #Teilchendichte
A=1.0 #Normierung des Beispielspektrums
mb=1e-27 #Umrechnung von mb in cm^2
p=1.0 #proton spectral index
Eth=1.22e-3 #Threshold for pion production
epsilon=1e-3
def sigma(Ep):
L=np.log(Ep/1.0)
if Ep<=Eth:
return 0.0
return (34.3+1.88*L+0.25*L*L)*(1-(Eth/Ep)**4)**2
def J(Ep):
return A*Ep**(-2.1)*np.exp(-Ep/E0)
def F_gamma(x,Ep):
L=np.log(Ep/1.0)
B=1.3+0.14*L+0.011*L*L
b=1.0/(1.79+0.11*L+0.008*L*L)
k=1.0/(0.801+0.049*L+0.014*L*L)
xb=(x+epsilon)**b
lx=np.log(x+epsilon)
return B*lx/x*((1-xb)/(1+k*xb*(1-xb)))*((1-xb)/(1+k*xb*(1-xb)))*\
((1-xb)/(1+k*xb*(1-xb)))*((1-xb)/(1+k*xb*(1-xb)))*\
(1/lx-4*b*xb/(1-xb)-4*k*b*xb*(1-2*xb)/(1+k*xb*(1-xb)))
def delta(x,Emin):
return
def F1_gamma(Ep,Eg):
return n*c*sigma(Ep)*F_gamma(Ep,Eg/Ep)*J(Ep)/Ep*mb
def flux2(Eg):
f2=0.0
e0 = np.logspace(math.log10(Eg),12, 131)
for n in range(0, len(e0)-2):
DE = e0[n+1] - e0[n]
EC = e0[n] + 0.5*DE
if J(EC) == 0:
break
f2 += J(EC) * DE *n*c*sigma(EC)*F_gamma(Eg/EC,EC)/EC*mb
return f2
energy=np.logspace(-1,12,131)
f=np.zeros(len(energy))
for n,nrg in enumerate(energy):
f[n]=flux2(nrg)
solution=np.column_stack((energy,f*energy**2.0))
np.savetxt("kelner1.csv",solution,delimiter=";")
print solution
Was ich bisher gemacht habe:
1) ich habe mittels integrate.quad versucht die Funktion zu integrieren. Allerdings musste ich realisieren, dass quad mit singularitäten nicht gut klar kommt.
2) Ich habe die Obersumme berechnet. Deshalb sind die Abweichungen..
Ich würde gerne wissen, ob ihr vllt einen Tipp hat, wie man das Problem am besten lösen könnte.
Für solche Tipps wäre ich sehr dankbar.
Die besten Grüße