gaussbesseltest.py 7.95 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
#-*- 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
#########################################################


42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64

######## Função de convolução com Gaussiana #######
def convoluiGauss(espectro, Sigma,passoE):

    areaantes=espectro.sum()
    temp=zeros(len(espectro),float)
    wg4_step=4.0*Sigma/passoE;
    wg2=Sigma**2
    
    eixo=arange(-6*Sigma,6*Sigma,passoE)
    eixo=arange(0,len(temp)*passoE,passoE)
    print len(eixo)
    
    for i in range(0,len(eixo)):
        for j in range(0,i):
            temp[i]+=espectro[j]*exp(-1.0*pow(passoE*(i-j-wg4_step),2.0)/(wg2*2.0));
            
    return temp*areaantes/temp.sum(), wg4_step*passoE
#    return temp/temp.max(), wg4_step*passoE  # Incluido para teste com Gaussiana -> Altura = 1
#########################################################



65 66
#### Parâmetros iniciais ####
# Feixe:
67
E0=99000 #eV
68 69 70
mi=1.
Zi=1.
# Alvo:
71 72 73 74
dedx=230. # eV/nm
dw2dx=19000. # eV^2/nm
Theta_in=7.
Theta_out=58.67
75
Theta_s = PI - (Theta_in + Theta_out)*PI/180
76
espessura=19
77 78 79 80

# Cálculo:
deltaE=20000
passoE=20
81
passox=0.01
82 83 84 85
e = arange(E0-deltaE, E0, passoE)
totalbessel=e*0
totalbessel_ls=e*0
totalgaussiana=e*0
86
resolucao = 350 # Largura a meia algura da Gaussiana que representa a resolução experimental (FWHM) em eV
87 88 89 90 91 92 93

## 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 ##
94 95 96
elementos=[[57.,138.9,192.,"la.prof"],[38.,87.62,165.,"sr.prof"],[22.,47.867,142.,"sr.prof"]] # La, Sr,ti
#elementos=[[72.,178.49,217.,"hf.txt"],[14.,28.,121.,"si.txt"],[8.,16.,102,"o.txt"]] # Hf, Si e O
#elementos=[[72.,178.49,17.],[14.,28.,21.],[8.,16.,12]] # Hf, Si e O
97 98 99 100 101 102 103 104 105 106 107 108
#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]
109 110
   perfil = loadtxt(componente[3], unpack=True)
#   plot(perfil)
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136

   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):
137
       secchoque=CSRf(Zi, Zt,k*E0-dedxef*x,mi,mt, Theta_s)*passox*perfil[int(x/passox)]
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
                              # É 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
161
       sigma=(sqrt(x*dw2dxef)+resolucao/2.35482)**2
162 163 164 165 166 167 168 169 170 171 172 173 174
       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
175
   totalbessel_ls_gauss,gshift=convoluiGauss(totalbessel_ls,resolucao/2.35482,passoE)
176 177 178 179
   totalgaussiana=sinalgaussiana+totalgaussiana
#####################################

# Mostra o Gráfico
180 181 182
#bes=plt.plot(e,totalbessel,lw=2,color="red")
#besls=plt.plot(e,totalbessel_ls,lw=2,color="blue")
beslsgaus=plt.plot(e-gshift,totalbessel_ls_gauss,lw=2,color="orange")
183 184 185 186
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")
187 188 189
#plt.legend([bes,besls,beslsgaus,gaus],["Bessel","Bessel+lineshape","Bessel+ls+expres","Gaussiana"],loc= "upper left", shadow=True)  
plt.legend([beslsgaus,gaus],["Bessel+ls+expres","Gaussiana"],loc= "upper left", shadow=True)  
#plt.legend([bes,besls,gaus],["Bessel","Bessel+lineshape","Gaussiana"],loc= "upper left", shadow=True)  
190 191
plt.xlabel("Energia dos Ions (eV)")
plt.ylabel("Rendimento (contagens)")
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225




################################################
# Testando convolução com deltas de Dirac
#figure()
#xxy=e*0
#xxy[00]=1. # Correspondente a uma delta em 80000

#xxy[100]=1. # Correspondente a uma delta em 82000
#xxy[300]=1.
#xxy[400]=1.
#xxy[500]=1.

# Convoluindo as deltas com uma gaussiana de FWHM = 550 eV
#xxyz,g_shift=convoluiGauss(xxy,resolucao/2.35482,passoE)

# Convoluindo a convolução novamente -> FWHM Res = 777.81 eV = sqrt(550^2+550^2)
#xxyz2,g_shift2=convoluiGauss(xxyz,resolucao/2.35482,passoE)
#deltas = plot(e,xxy,lw=2)
#convo1 = plot(e-g_shift,xxyz+0.001,lw=2)
#convo2 = plot(e-g_shift-g_shift2,xxyz2+0.002,lw=2)
#convo2_analitica=plot(e,exp(-(e-82000)**2.0/(2.*(777.81/2.35482)**2))+0.002,lw=2)
#plt.legend([deltas,convo1,convo2,convo2_analitica],["Deltas","Conv Gauss 1x 550 eV","Conv Gauss 2x 550 eV","Conv Analitica (777.81)"],loc= "upper right", shadow=True)  

#print e-g_shift, xxyz
#print str(g_shift) + " <- g_shift"
#print len(e)
#print len(xxy)
#print len(xxyz)



226 227
show() 
##################