Commit 9f28f796 authored by Rafael Peretti Pezzi's avatar Rafael Peretti Pezzi

Comparação dos espectros Bessel e Gaussiano, incluido Seção de Choque

parent 324d41fb
#-*- coding: utf-8 -*- #-*- coding: utf-8 -*-
# Teste para a perda de energia baseada na função modificada de Bessel do primeiro tipo de ordem 1. # Comparação do espectros de espalhamento baseados na função modificada de Bessel do primeiro tipo de ordem 1
# Apresenta os graficos para a perda de energia para algumas espessuras na geometria do experimento # com o modelo Gaussiano.
# comparadas com as correspondentes perdas de energias Gaussianas nas mesmas espessuras.
# Referência: R. P. Pezzi, et al. Applied Physics Letters, 92 164102 (2008) # Referência: R. P. Pezzi, et al. Applied Physics Letters, 92 164102 (2008)
from pylab import * from pylab import *
from scipy.special import i1 from scipy.special import i1
PI=math.pi 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
###########################################
#### Parâmetros iniciais ####
# Feixe:
E0=100000 #eV E0=100000 #eV
deltaE=25000 mi=1.
passo=30 Zi=1.
# Alvo:
dedx=192. # eV/nm dedx=192. # eV/nm
dw2dx=20740. # eV^2/nm dw2dx=20740. # eV^2/nm
Theta_in=0. Theta_in=0.
Theta_out=70. Theta_out=70.
Theta_s = PI - (Theta_in + Theta_out)*PI/180 Theta_s = PI - (Theta_in + Theta_out)*PI/180
mt = 178.49 espessura=2
mi = 1.
k = ((sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2 # Cálculo:
deltaE=25000
passoE=20
passox=0.02
e = arange(E0-deltaE, E0, passoE)
## 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 ##
elementos=[[72.,178.49],[14.,28.],[8.,16.]] # 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];
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
dedxef=dedx*(k/cos(Theta_in*PI/180.) + 1/cos(Theta_out*PI/180.)) # Calcula a contribuição de cada camada
dw2dxef=dw2dx*(k**2./cos(Theta_in*PI/180.)+1./cos(Theta_out*PI/180.)) for x in arange(passox,espessura,passox):
e = arange(E0-deltaE, E0, passo)
for x in arange(0,30,1): # 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.
# Bessel: # Também temos um problema com a primeira camada, que tem profundidade zero.
# 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 # Perceba que os cálculos de Bessel correspondentes às camadas mais fundas estão truncados
# da precisão numérica. Uma possível solução pode ser trabalhar em outras unidades que não # para uma perda de energia maior que 25 keV. (Observado com Theta_in=70, para x > 30)
# resultem em argumentos tão grandes para a função exponencial.
# Perceba que os cálculos de Bessel correspondentes às camadas mais fundas estão truncados. # CSRf(Z1, Z2, Ein, M1, M2, Theta):
secchoque=CSRf(Zi, Zt,k*E0-dedxef*x,mi,mt, Theta_s)*passox # É 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)))
sinalbessel=besselcamada*secchoque+sinalbessel
#plt.plot(e,besselcamada)
alpha=dedxef*(2./dw2dxef) # Gaussiana:
m=alpha*dedxef Em=k*E0-x*dedxef
lbd=m*x*alpha sigma=x*dw2dxef
bessel = lbd*exp(-m*x-alpha*(k*E0-e))*i1(2.*sqrt(lbd*(k*E0-e)))/(sqrt(lbd*(k*E0-e))) gaussianacamada = exp((-(e-Em)**2)/(2.*sigma))/sqrt(2.*PI*sigma)
plt.plot(e, bessel) sinalgaussiana=gaussianacamada*secchoque+sinalgaussiana
#plt.plot(e,gaussianacamada)
# Gaussiana: plt.plot(e, sinalbessel)
Em=k*E0-x*dedxef plt.plot(e, sinalgaussiana)
sigma=x*dw2dxef #####################################
gaussian = exp((-(e-Em)**2)/(2.*sigma))/sqrt(2.*PI*sigma)
plt.plot(e,gaussian)
# Mostra o Gráfico
show() show()
##################
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment