gaussbesseltest.py 1.69 KB
Newer Older
1 2
#-*- coding: utf-8 -*-
# Teste para a perda de energia baseada na função modificada de Bessel do primeiro tipo de ordem 1.
3
# Apresenta os graficos para a perda de energia para algumas espessuras na geometria do experimento
4
# comparadas com as correspondentes perdas de energias Gaussianas nas mesmas espessuras.
5
# Referência: R. P. Pezzi, et al. Applied Physics Letters, 92 164102 (2008)
6 7 8

from pylab import *
from scipy.special import i1
9
PI=math.pi
10

11 12
E0=100000 #eV
deltaE=25000
13
passo=30
14

15 16
dedx=192. # eV/nm
dw2dx=20740. # eV^2/nm
17 18 19 20
Theta_in=0.
Theta_out=70.
Theta_s = PI - (Theta_in + Theta_out)*PI/180
mt = 178.49
21 22
mi = 1.
k = ((sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2	
23 24 25 26


dedxef=dedx*(k/cos(Theta_in*PI/180.) + 1/cos(Theta_out*PI/180.))
dw2dxef=dw2dx*(k**2./cos(Theta_in*PI/180.)+1./cos(Theta_out*PI/180.))
27
e = arange(E0-deltaE, E0, passo)
28

29
for x in arange(0,30,1): 
30 31 32 33 34 35 36

    # Bessel:  
    # Ainda falta incluir a delta de Dirac na origem e a convolução com a resolução experimental
    # Existe um valor máximo em energia no qual a função de Bessel pode ser computada em função
    # da precisão numérica. Uma possível solução pode ser trabalhar em outras unidades que não
    # resultem em argumentos tão grandes para a função exponencial.

37
    # Perceba que os cálculos de Bessel correspondentes às camadas mais fundas estão truncados.
38

39 40
    alpha=dedxef*(2./dw2dxef)
    m=alpha*dedxef
41
    lbd=m*x*alpha
42 43
    bessel = lbd*exp(-m*x-alpha*(k*E0-e))*i1(2.*sqrt(lbd*(k*E0-e)))/(sqrt(lbd*(k*E0-e)))
    plt.plot(e, bessel)  
44 45

    # Gaussiana:
46
    Em=k*E0-x*dedxef
47
    sigma=x*dw2dxef
48
    gaussian = exp((-(e-Em)**2)/(2.*sigma))/sqrt(2.*PI*sigma) 
49
    plt.plot(e,gaussian)
50 51

show()