calc.py 8.97 KB
Newer Older
1 2 3 4 5 6 7
#-*- 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, zeros
8
from scipy.special import i1, jn
9 10 11 12 13 14 15 16
from Ewindow import Edict
from layers import Ldict
PI=math.pi

####################################################################################################################################################
# Seção de choque Rutherford
def CSRf(Z1, Z2, Ein, M1, M2, Theta):
    if M1 < M2:
17
        return ((Z1*Z2*4.8e-20*1.e8)/(4.*Ein))**2. * sin(Theta)**-4.
18
    else:
19
        return 0.
20 21 22 23

####################################################################################################################################################
# Fator cinematico
def kinematic(mt, mi, esp):
24 25
    MM = mi/mt
    return ((sqrt (1. - MM**2.*sin(esp)**2.) + MM*cos(esp) ) / (1. + MM) )**2.
26

27 28 29 30 31
####################################################################################################################################################
# Coeficientes de Poisson
def kpoisson(n, m, x):
    return (m*x)**n*exp(-m*x)/math.factorial(n)  

32 33 34 35 36 37
####################################################################################################################################################
# Funcoes em comum

def spect_file(x, y, shift, dose):
    data2 = open('temp.sim', 'w')
    for i in range(len(y)):
38
        print>>data2, x[i]/1000.,y[i]
39 40
    data2.close()

41 42 43
def set_lim(x, x2, y, shift):  
    ylim([0,y*1.1])          
    xlim([x*0.95,x2*1.05])
44 45 46 47

def variables(OFlabentries):
    EPasso = float(OFlabentries[16].get())
    passox = float(OFlabentries[17].get())
48 49 50
    Emin = float(OFlabentries[14].get())
    Emax = float(OFlabentries[15].get())
    e = arange(Emin, Emax, EPasso)
51 52
    resol = float(OFlabentries[13].get())
    dose = float(OFlabentries[24].get())
53
    return EPasso, passox, e, resol, dose, Emin, Emax
54

55 56 57 58 59 60 61 62
def labeler(plotes, botoes, fig):
    labels = list()
    grafs = list()
    for l in range(len(botoes)):
        labels.insert(l, botoes[l]['text'])
        grafs.insert(l, plotes[l])
    fig.legend(grafs, labels)

63 64 65
####################################################################################################################################################
####################################################################################################################################################

66
def calc(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo, espectro, emin, plote, fig, firsttime, botoes, indice):
67

68
    EPasso, passox, e, resol, dose, Emin, Emax = variables(OFlabentries)
69
    espectro.clear()
70
    espectro.eshift = float(emin)
71

72 73 74
    gshift=0.

    Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180.
75 76 77 78 79
    mi = ion['mass']

    if metodo == 'layer': pass
    elif metodo == 'draw':
        for componente in Edict:
80

81
            mt = Edict[componente]['mass']
82
            c = Edict[componente]['dist'].data
83 84 85 86
            alpha_lineshape = Edict[componente]['LineShape']

            k = kinematic(mt, mi, Theta_s)

87
            dedxef=p['dedx']*(k/cos(p['Theta_in']*PI/180.) + 1./cos(p['Theta_out']*PI/180.))
88 89
            dw2dxef=p['dW2dx']*(k**2./cos(p['Theta_in']*PI/180.)+1./cos(p['Theta_out']*PI/180.))

90 91
            sinalbessel = e*0.
            sinalgaussiana = e*0.
92 93

            # Calcula a contribuição de cada camada
94
            for x in arange(0, Edict[componente]['profundidademax'], passox): 
95
 
96 97 98 99 100 101 102 103 104 105 106 107
                # 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:  
                # 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':
108
                    alpha=2.*dedxef/dw2dxef
109
                    m=alpha*dedxef
110 111 112 113
                    lbd=m*x*alpha
                    espectro.spectro[int(k*p['E0']/EPasso-1)] = CS*exp(-m*x)+x*m*alpha*exp(-m*x)*i1(0)/sqrt(0.00000000001)*passox*passox
                    besselcamada = lbd*exp(-m*x-alpha*(k*p['E0']-e))*i1(2.*sqrt(lbd*(k*p['E0']-e)))/(sqrt(lbd*(k*p['E0']-e)+.0000001))
                    espectro.soma( list(nan_to_num(besselcamada*CS))) 
114 115 116 117 118

                # Gaussiana:
                elif modelo == 'gaussiana':
                    Em=k*p['E0']-x*dedxef
                    sigma=x*dw2dxef
119 120 121 122 123 124
                    gaussianacamada = exp((-(e-Em)**2.)/(2.*sigma))/sqrt(2.*PI*sigma) 
                    espectro.soma( list(nan_to_num(gaussianacamada*CS)) )
                 
            if LScontrol == 1:
                espectro.ConvolveF0(1./alpha_lineshape)
            
125 126 127
            if modelo == 'bessel':
                linecolor='r'

128
            elif modelo == 'gaussiana':         
129 130 131 132 133 134
                linecolor='b'

        if CGausscontrol == 1:
            espectro.ConvolveGauss(resol/2.35482)    

    espectro.multiply(dose)
135 136
    plote[indice].set_xdata(espectro.axisenergy())           
    plote[indice].set_ydata(espectro.spectro)
137 138 139

    spect_file(espectro.axisenergy(), espectro.spectro, espectro.eshift, dose)

140 141 142 143
    ylim([0,max(espectro.spectro)*1.1])          
    if firsttime == True:
        xlim([Emin,Emax])

144
    labeler(plote,botoes,fig)
145 146
    fig.canvas.draw()
    fig.show()
147

148 149 150
####################################################################################################################################################
####################################################################################################################################################

151
def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo, espectro, emin, plote, fig, firsttime, botoes, indice, aux, aux2):
152 153

    EPasso, passox, e, resol, dose, Emin, Emax = variables(OFlabentries)
154
    passox = 0.005
155 156 157
    gshift=0.
    curdepht=passox    #Avoid zero division
    espectro.clear()
158 159 160
    aux.clear()
    aux2.clear()
    aux2.spectro[0] = 1./EPasso
161

162
    Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180.
163 164 165 166 167 168 169 170 171 172 173 174 175 176
    mi = ion['mass']

    for camada in Ldict:

        xfinal = curdepht + float(camada.LandE[5].get())

        ## Calcula o espectro de cada elemento da amostra
        for componente in range(len(camada.Elements)):

            mt = camada.Elements[componente]['mass']
            c = float(camada.Spinboxes[componente].get())
            alpha_lineshape = camada.Elements[componente]['LineShape']
            dedx = float(camada.LandE[3].get())
            dW2dx = float(camada.LandE[4].get())
177

178 179
            k = ((sqrt (mt**2.-mi**2.*(sin(Theta_s)**2.)) + mi*cos(Theta_s) ) / (mi+mt) )**2.	
            dedxef=dedx*(k/cos(p['Theta_in']*PI/180.) + 1./cos(p['Theta_out']*PI/180.))
180 181 182 183 184 185 186 187 188
            dw2dxef=dW2dx*(k**2./cos(p['Theta_in']*PI/180.)+1./cos(p['Theta_out']*PI/180.))
            
            # Calcula a contribuição de cada passo
            for x in arange(curdepht, xfinal, passox): 

                # Seção de choque - CSRf(Z1, Z2, Ein, M1, M2, Theta):
                CS = c*k*CSRf(ion['Z'], camada.Elements[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.

189 190
                alpha=dedxef*(2./dw2dxef)
                m=alpha*dedxef
191

192 193 194 195
                aux.clear()
                for y in range(20):
                    aux.soma(kpoisson((y+1.),m,passox)*alpha*exp(-alpha*aux.axisenergy()))
                aux.spectro[0] = kpoisson(0.,m,passox)*1./EPasso
196

197 198 199 200 201 202
                aux2.Convolve(aux)
                aux2.multiply(CS)
                espectro.soma(aux2.spectro)
                aux2.multiply(1./CS)
         
            linecolor='b'
203 204 205

        curdepht = xfinal

206
    espectro.reverse()
207 208 209 210 211
    if LScontrol == 1:
        espectro.ConvolveF0(1./alpha_lineshape)
    if CGausscontrol == 1:
        espectro.ConvolveGauss(resol/2.35482) 

212 213
    espectro.eshift = espectro.eshift - (1-k)*p['E0'] - EPasso
    espectro.multiply(dose*0.93)
214 215 216 217 218 219 220 221 222 223 224 225
    plote[indice].set_xdata(espectro.axisenergy())           
    plote[indice].set_ydata(espectro.spectro)

    spect_file(espectro.axisenergy(), espectro.spectro, espectro.eshift, dose)

    ylim([0,max(espectro.spectro)*1.1])          
    #if firsttime == True:
     #   xlim([Emin,Emax])

    labeler(plote,botoes,fig)
    fig.canvas.draw()
    fig.show()