Incluída convolução com F0. Fix: NaN no Vetor.

	Melhor aparência do gráfico

	Incluído programa pn_x.py que apresenta os coeficientes
           de Poisson para ions com espectro de colisão dado por
           decaimento exponencial.

gaussbesseltest.py

 #-*- 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 *
@@ -13,8 +15,28 @@ def CSRf(Z1, Z2, Ein, M1, M2, Theta):
         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 ####
@@ -32,17 +54,24 @@ Theta_s = PI - (Theta_in + Theta_out)*PI/180
 espessura=2
 
 # Cálculo:
-deltaE=25000
+deltaE=20000
 passoE=20
 passox=0.02
 e = arange(E0-deltaE, E0, passoE)
+totalbessel=e*0
+totalbessel_ls=e*0
+totalgaussiana=e*0
+
+
 ## 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
+elementos=[[72.,178.49,217.],[14.,28.,121.],[8.,16.,102]] # 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 
 ####################################
@@ -52,6 +81,7 @@ 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.))
@@ -73,15 +103,31 @@ for componente in elementos:
        
        # 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 # É multiplicada pelo passo para manter a escala do gráfico
-                                                                  # independente do mesmo.
+       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)))
+
+       ### 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:
@@ -91,10 +137,28 @@ for componente in elementos:
        sinalgaussiana=gaussianacamada*secchoque+sinalgaussiana
        #plt.plot(e,gaussianacamada)
 
-   plt.plot(e, sinalbessel)
-   plt.plot(e, sinalgaussiana)  
+#   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
+bes=plt.plot(e,totalbessel,lw=2,color="red")
+besls=plt.plot(e,totalbessel_ls,lw=2,color="blue")
+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,gaus],["Bessel","Bessel+lineshape","Gaussiana"],loc= "upper left", shadow=True)  
+plt.xlabel("Energia dos Ions (eV)")
+plt.ylabel("Rendimento (contagens)")
+#print totalbessel
 show() 
 ##################

pn_x.py

+#-*- coding: utf-8 -*-
+# Estatística de Poisson para feixes com espectro de colisão dado por decaimento exponencial
+# f*¹(E)=alpha*exp(-alpha*E) p/ E > 0
+#
+# Autor: Rafael Pezzi
+# Data: Set/2011
+# Licença: GPL 3.0
+
+def fatorial(numero):
+
+   resultado = 1
+
+   lista = range(1,n+1)
+
+   for x in lista:
+
+     resultado = x * resultado
+   return resultado
+
+from pylab import *
+ion()
+dedx=192. # eV/nm
+dw2dx=20740. # eV^2/nm
+nmax=45
+
+enes = arange(0,nmax,1)
+pn=zeros(nmax)
+profundidade=2
+width=0.8
+m = 2*dedx**2/dw2dx
+rodada=0.
+plotes=[]
+
+print "m =",m
+#print len(pn)
+#print len(enes)
+
+for profundidade in [.1,.25,.5,1,2.5,5,7.5]:
+
+ for n in enes:
+    pn[n]=((m*profundidade)**n)*exp(-m*profundidade)/fatorial(n)
+
+ figure()
+ subplot(111)
+ plt.title("Coeficientes de Poisson para m="+str(m)+"/nm")
+ plt.ylabel('$P_n$')
+ plt.xlabel('Numero de Colisoes')
+ if rodada==0: cor='red'
+ if rodada==1: cor='green'
+ if rodada==2: cor='orange'
+ if rodada==3: cor='yellow'
+ if rodada==4: cor='cyan'
+ if rodada==5: cor='black'
+ if rodada==6: cor='blue'
+
+ plote=bar(enes, pn,width, color=cor)
+ plotes+=plote
+ plt.ylim=(0.,1.)
+ draw()
+ rodada+=1
+ plt.legend( [plote[0]] ,[str(profundidade)+" nm"], title="Profundidade:")  
+
+show()