gaussbesseltest.py 5.51 KB
Newer Older
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
#-*- coding: utf-8 -*-
# Comparação do espectros de espalhamento baseados na função modificada de Bessel do primeiro tipo de ordem 1 
# com o modelo Gaussiano. 
# Incluída convolução com forma de linha originária da colisão frontal
#
# Referência: R. P. Pezzi, et al. Applied Physics Letters, 92 164102 (2008)

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

######## Definindo a Seção de Choque #######
def CSRf(Z1, Z2, Ein, M1, M2, Theta):
    if M1 < M2:
        return ((Z1*Z2*4.8e-20*1e8)/(4*Ein))**2 * sin(Theta)**-4 * (sqrt(1-((M1/M2)*sin(Theta))**2)+cos(Theta))**2 / sqrt(1-((M1/M2)*sin(Theta))**2) 
    else:
        return 0
############################################

######## Função de convolução com forma de linha #######
def convoluiF0(espectro, alpha,passoE):

    areaantes=espectro.sum()

    temp=espectro[::-1]
    eixo=arange(0,12.*alpha,passoE)
    lineshape=alpha*exp(-eixo/alpha)

#    ls_plt=plt.plot(lineshape/lineshape.sum(),lw=2)
#    sig_plt=plt.plot(temp/temp.sum(),lw=2)

    temp = convolve(temp/temp.sum(), lineshape/lineshape.sum())
    temp.resize(len(espectro))

#    result_plt=plt.plot(temp/temp.sum(),'o')
#    legend([ls_plt,sig_plt,result_plt],["Lineshape","Sinal","Resultado"])

    return temp[::-1]*areaantes
#########################################################


#### Parâmetros iniciais ####
# Feixe:
44
E0=100000 #eV
45 46
mi=1.
Zi=1.
47

48
# Alvo:
49 50 51 52
dedx=192. # eV/nm
dw2dx=20740. # eV^2/nm
Theta_in=0.
Theta_out=70.
53
Theta_s = PI - (Theta_in + Theta_out)*PI/180
54
espessura=2
55 56 57 58

# Cálculo:
deltaE=20000
passoE=20
59
passox=0.02
60 61 62 63
e = arange(E0-deltaE, E0, passoE)
totalbessel=e*0
totalbessel_ls=e*0
totalgaussiana=e*0
64

65 66 67 68 69 70 71

## TODO: - Otimizar os cálculos de cada componente: evitar cálculos (de zeros) em energias desnecessários
##       - Definir o intervalo de plotagem do espectro independente do deltaE
#############################


## Define a composição da amostra ##
72
elementos=[[72.,178.49,217.],[14.,28.,121.],[8.,16.,102]] # Hf, Si e O
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
#elementos=[[72.,178.49,217.]]#,[14.,28.,121.],[8.,16.,102]] # Hf, Si e O

# Esta lista deve incluir a distribuição
# em profundidade do elemento 
####################################

## Calcula o espectro de cada elemento da amostra
for componente in elementos:
   
   mt=componente[1];
   Zt=componente[0];
   alpha_lineshape=componente[2]

   k = ((sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2	
   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.))

   sinalbessel=e*0
   sinalgaussiana=e*0

   # Calcula a contribuição de cada camada
   for x in arange(passox,espessura,passox): 

       # 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.

       # Também temos um problema com a primeira camada, que tem profundidade zero.
       
       # Perceba que os cálculos de Bessel correspondentes às camadas mais fundas estão truncados
       #  para uma perda de energia maior que 25 keV. (Observado com Theta_in=70, para x > 30)
       # Este truncamento resulta de valores NaN (Not a Number) que aparecem no vetor. Foi incluida uma
       #  operação para remover os NaN do vetor de Bessel
       # 

       #       CSRf(Z1, Z2, Ein, M1, M2, Theta):
111
       secchoque=CSRf(Zi, Zt,k*E0-dedxef*x,mi,mt, Theta_s)*passox
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
                              # É multiplicada pelo passo para manter a escala do gráfico
                              # independente do mesmo.

       alpha=dedxef*(2./dw2dxef)
       m=alpha*dedxef
       lbd=m*x*alpha
       besselcamada = lbd*exp(-m*x-alpha*(k*E0-e))*i1(2.*sqrt(lbd*(k*E0-e)))/(sqrt(lbd*(k*E0-e)))

       ### Removendo Not A Number do vetor ###
       ondeestaoNaN = isnan(besselcamada);
       besselcamada[ondeestaoNaN] = 0

       ### Incluir o sinal correspondente à ausência de colisões (Delta de Dirac em k*E0)
       # Pode ser feito adicionando uma gaussiana de largura muito pequeno em k*E0 para circundar o 
       #  fato de kE0 poder não cair exatamente em uma posição discreta do vetor

       sinalbessel=besselcamada*secchoque+sinalbessel

       ## Plotar o sinal da camada
       #plt.plot(e,besselcamada)

       # Gaussiana:
       Em=k*E0-x*dedxef
135
       sigma=x*dw2dxef
136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
       gaussianacamada = exp((-(e-Em)**2)/(2.*sigma))/sqrt(2.*PI*sigma) 
       sinalgaussiana=gaussianacamada*secchoque+sinalgaussiana
       #plt.plot(e,gaussianacamada)

#   plt.plot(e, sinalbessel)
#   plt.plot(e, sinalgaussiana)  

   ### Convoluir com forma de linha da colisão frontal:
   sinalbessel_ls=convoluiF0(sinalbessel,alpha_lineshape,passoE)


   totalbessel=sinalbessel+totalbessel
   totalbessel_ls=sinalbessel_ls+totalbessel_ls
   totalgaussiana=sinalgaussiana+totalgaussiana
#####################################

# Mostra o Gráfico
153 154
bes=plt.plot(e,totalbessel,lw=2,color="red")
besls=plt.plot(e,totalbessel_ls,lw=2,color="blue")
155 156 157 158
gaus=plt.plot(e,totalgaussiana,lw=2,color="green")

#plt.plot(e,totalbessel*10,lw=2, color="blue")
#plt.plot(e,totalgaussiana*10,lw=2, color="green")
159
plt.legend([bes,besls,gaus],["Bessel","Bessel+lineshape","Gaussiana"],loc= "upper left", shadow=True)  
160 161
plt.xlabel("Energia dos Ions (eV)")
plt.ylabel("Rendimento (contagens)")
162
#print totalbessel
163 164
show() 
##################