#-*- 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 ######################################################### ######## 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 ######################################################### #### Parâmetros iniciais #### # Feixe: E0=99000 #eV mi=1. Zi=1. # Alvo: dedx=230. # eV/nm dw2dx=19000. # eV^2/nm Theta_in=7. Theta_out=58.67 Theta_s = PI - (Theta_in + Theta_out)*PI/180 espessura=19 # Cálculo: deltaE=20000 passoE=20 passox=0.01 e = arange(E0-deltaE, E0, passoE) totalbessel=e*0 totalbessel_ls=e*0 totalgaussiana=e*0 resolucao = 350 # Largura a meia algura da Gaussiana que representa a resolução experimental (FWHM) em eV ## 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=[[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 #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] perfil = loadtxt(componente[3], unpack=True) # plot(perfil) 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): secchoque=CSRf(Zi, Zt,k*E0-dedxef*x,mi,mt, Theta_s)*passox*perfil[int(x/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))) ### 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 sigma=(sqrt(x*dw2dxef)+resolucao/2.35482)**2 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 totalbessel_ls_gauss,gshift=convoluiGauss(totalbessel_ls,resolucao/2.35482,passoE) totalgaussiana=sinalgaussiana+totalgaussiana ##################################### # Mostra o Gráfico #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") 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") #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) plt.xlabel("Energia dos Ions (eV)") plt.ylabel("Rendimento (contagens)") ################################################ # 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) show() ##################