numerical.py 6.23 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
#-*- 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. 
# Referência: R. P. Pezzi, et al. Applied Physics Letters, 92 164102 (2008)

from pylab import *
from numpy import nan_to_num
from scipy.special import i1
from Ewindow import Edict

PI=math.pi

####################################################################################################################################################
# Seção de choque Rutherford
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

####################################################################################################################################################
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
# 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)
    
    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

####################################################################################################################################################

62
def spectro(p, modelo, ion, Emin, Emax, EPasso, passox, LScontrol, resol, espessura):
63 64

    e = arange(Emin, Emax, EPasso)
65 66
    if modelo == 'bessel':
        Y = e*0.
67
        '''try:
68 69 70 71
            YY
        except NameError: 
            YY, = plot(e,e*0., 'b')
        else:
72
            YY.set_xdata(e)'''
73 74 75

    elif modelo == 'gaussiana':
        X = e*0.
76
        '''try:
77 78
            XX
        except NameError: 
79
            #XX, = plot(e,e*0., 'r')
80
        else:
81
            #XX.set_xdata(e)'''
82

83
    dose = 12e34
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 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

    Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180
    mi = ion['mass']

    ## Calcula o espectro de cada elemento da amostra
    for componente in Edict:
   
        mt = Edict[componente]['mass']
        c = Edict[componente]['dist']
        alpha_lineshape = Edict[componente]['LineShape']

        k = ((sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2	
        dedxef=p['dedx']*(k/cos(p['Theta_in']*PI/180.) + 1/cos(p['Theta_out']*PI/180.))
        dw2dxef=p['dW2dx']*(k**2./cos(p['Theta_in']*PI/180.)+1./cos(p['Theta_out']*PI/180.))

        sinalbessel=e*0
        sinalgaussiana=e*0

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

            # Seção de choque - CSRf(Z1, Z2, Ein, M1, M2, Theta):
            CS = c[int(x/passox)]*k*CSRf(ion['Z'], Edict[componente]['Z'], k*p['E0']-dedxef*x, mi, mt, Theta_s)*passox 
            # É multiplicada pelo passo para manter a escala do gráfico independente do mesmo.

            # 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)
            if modelo == 'bessel':
                alpha=dedxef*(2./dw2dxef)
                m=alpha*dedxef
                lbd=m*x*alpha
                besselcamada = lbd*exp(-m*x-alpha*(k*p['E0']-e))*i1(2.*sqrt(lbd*(k*p['E0']-e)))/(sqrt(lbd*(k*p['E0']-e)))
                sinalbessel = besselcamada*CS+sinalbessel

            # Gaussiana:
            elif modelo == 'gaussiana':
                Em=k*p['E0']-x*dedxef
                sigma=x*dw2dxef
                gaussianacamada = exp((-(e-Em)**2)/(2.*sigma))/sqrt(2.*PI*sigma) 
                sinalgaussiana = gaussianacamada*CS+sinalgaussiana

        if modelo == 'bessel':
            besself = nan_to_num(sinalbessel) 
            if LScontrol == 1:
                besself = convoluiF0(besself, alpha_lineshape, EPasso)
            Y = besself + Y
136
            Y, gshiftY = convoluiGauss(Y,resol/2.35482,EPasso)
137 138 139
        elif modelo == 'gaussiana':
            if LScontrol == 1:     
                sinalgaussiana = convoluiF0(sinalgaussiana, alpha_lineshape, EPasso)           
140 141 142 143
            X = sinalgaussiana + X
            X, gshiftX = convoluiGauss(X,resol/2.35482,EPasso)

    if modelo == 'bessel':
144 145
        plt.plot(e-gshiftY,Y*dose, 'b')
        #YY.set_ydata(Y-*dose)
146 147
        ylim([0,max(Y*dose)*1.1])
    elif modelo == 'gaussiana':
148 149
        plt.plot(e-gshiftX,Y*dose, 'r')
        #XX.set_ydata((X)*dose)
150 151
        ylim([0,max(X*dose)*1.1])
    #plt.plot(e,Y*dose)
152 153 154 155
    suptitle('Signal received', fontsize=12)
    xlabel("Enegy (eV)")
    ylabel("Yield")

156 157 158
    amostra = np.loadtxt('H2.py')
    xxx= amostra[:,0]*1000
    ccc=amostra[:,1]
159

160
    plot (xxx, ccc, 'r.')
161 162 163 164 165

    show()

####################################################################################################################################################