numerical.py 13.4 KB
Newer Older
Matheus Müller's avatar
Matheus Müller committed
1
2
3
4
5
6
#-*- 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 *
7
from numpy import nan_to_num, zeros
Matheus Müller's avatar
Matheus Müller committed
8
9
from scipy.special import i1
from Ewindow import Edict
10
from layers import Ldict
Matheus Müller's avatar
Matheus Müller committed
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
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
# 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
    
51
    #eixo=arange(-6.*Sigma,6.*Sigma,passoE)
52
    eixo=arange(0,len(temp)*passoE,passoE)
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

    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


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

    areaantes=espectro.sum()
    temp=zeros(len(espectro),float)
    wg4_step=4.0*Sigma/passoE;
    wg2=Sigma**2
68
    
69
70
71
72
    #eixo=arange(-6.*Sigma,6.*Sigma,passoE)
    eixo=arange(0,len(temp)*passoE,passoE)
    gaussi = exp(-1.0*pow(passoE*(eixo-wg4_step),2.0)/(wg2*2.0))

73
74
    for i in range(0,len(eixo)):
        for j in range(0,i):
75
            temp[i]+=espectro[j]*gaussi[i-j]
76
77
78
            
    return temp*areaantes/temp.sum(), wg4_step*passoE

Matheus Müller's avatar
Matheus Müller committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
####################################################################################################################################################
####################################################################################################################################################
# Funcoes em comum

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

def set_lim(x, y, shift, dose):  
    ylim([0,max(y*dose)*1.1])          
    xlim([min(x-shift)*0.95,max(x-shift)*1.05])

93
####################################################################################################################################################
94
####################################################################################################################################################
95

96
def spectro(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, tempor):
Matheus Müller's avatar
Matheus Müller committed
97

98
99
100
101
    EPasso = float(OFlabentries[16].get())
    passox = float(OFlabentries[17].get())
    e = arange(float(OFlabentries[14].get()), float(OFlabentries[15].get()), EPasso)
    resol = float(OFlabentries[13].get())
102
    dose = float(OFlabentries[24].get())
103
    gshift=0
Matheus Müller's avatar
Matheus Müller committed
104
    X = e*0.
Matheus Müller's avatar
Matheus Müller committed
105
106
107
108
109
110
111
112
113
114
115

    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']

116
117
118
119
120
        k = ((sqrt (mt**2-mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2	

        #k = mi*cos(Theta_s) + sqrt (mt**2 - mi**2*sqrt(sin(Theta_s)))  
        #k = k*k/(mi+mt)/(mi+mt)

Matheus Müller's avatar
Matheus Müller committed
121
122
123
        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.))

124
125
        sinalbessel = e*0
        sinalgaussiana = e*0
Matheus Müller's avatar
Matheus Müller committed
126
127

        # Calcula a contribuição de cada camada
128
        for x in arange(passox, Edict[componente]['profundidademax'], passox): 
Matheus Müller's avatar
Matheus Müller committed
129
130
131
132
133
134

            # 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:  
135
            # Ainda falta incluir a delta de Dirac na origem
Matheus Müller's avatar
Matheus Müller committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
            # 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)
160
161
162
            X = besself + X
            linecolor='r'

Matheus Müller's avatar
Matheus Müller committed
163
164
165
        elif modelo == 'gaussiana':
            if LScontrol == 1:     
                sinalgaussiana = convoluiF0(sinalgaussiana, alpha_lineshape, EPasso)           
166
            X = sinalgaussiana + X
167
168
169
170
            linecolor='b'

    if CGausscontrol == 1:
        X, gshift = convoluiGauss(X,resol/2.35482,EPasso)    
171

172
173
174
    #tempor.set_xdata(e-gshift)           
    #tempor.set_ydata(X*dose)
    plot(e-gshift, X*dose)
Matheus Müller's avatar
Matheus Müller committed
175

Matheus Müller's avatar
Matheus Müller committed
176
177
178
    set_lim(e, X, gshift, dose)
    spect_file(e, X, gshift, dose)
    
179
    plt.show()
Matheus Müller's avatar
Matheus Müller committed
180
181

####################################################################################################################################################
182
####################################################################################################################################################
Matheus Müller's avatar
Matheus Müller committed
183

184
def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, tempor):
185
186
187
188
189

    EPasso = float(OFlabentries[16].get())
    passox = float(OFlabentries[17].get())
    e = arange(float(OFlabentries[14].get()), float(OFlabentries[15].get()), EPasso)
    resol = float(OFlabentries[13].get())
190
    dose = float(OFlabentries[24].get())
Matheus Müller's avatar
Matheus Müller committed
191
    X = e*0.
192
    gshift=0.
193
    curdepht=passox    #Avoid zero division
194
195
196
197
198
199

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

    for camada in Ldict:

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

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

205
206
207
            mt = camada.Elements[componente]['mass']
            c = float(camada.Spinboxes[componente].get())
            alpha_lineshape = camada.Elements[componente]['LineShape']
208
209
            dedx = float(camada.LandE[3].get())
            dW2dx = float(camada.LandE[4].get())
210

211
            k = ((sqrt (mt**2-mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2	
212
213
214
215
216
217
            dedxef=dedx*(k/cos(p['Theta_in']*PI/180.) + 1/cos(p['Theta_out']*PI/180.))
            dw2dxef=dW2dx*(k**2./cos(p['Theta_in']*PI/180.)+1./cos(p['Theta_out']*PI/180.))

            sinalbessel = e*0
            sinalgaussiana = e*0

218
219
            # Calcula a contribuição de cada passo
            for x in arange(curdepht, xfinal, passox): 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252

                # 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.

                # Bessel:  
                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)
                X = besself + X
                linecolor='r'

            elif modelo == 'gaussiana':
                if LScontrol == 1:     
                    sinalgaussiana = convoluiF0(sinalgaussiana, alpha_lineshape, EPasso)           
                X = sinalgaussiana + X
                linecolor='b'

253
254
        curdepht = xfinal

255
    if CGausscontrol == 1:
Matheus Müller's avatar
Matheus Müller committed
256
257
258
259
        X, gshift = convoluiGauss(X,resol/2.35482,EPasso)  
  
    tempor.set_xdata(e-gshift)           
    tempor.set_ydata(X*dose)
260

Matheus Müller's avatar
Matheus Müller committed
261
262
263
    set_lim(e, X, gshift, dose)
    spect_file(e, X, gshift, dose)
  
264
    plt.show()
265
266
267

####################################################################################################################################################
####################################################################################################################################################
268
269
270
271
272
273
274
275

def spectroRNA(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, tempor):

    EPasso = float(OFlabentries[16].get())
    passox = float(OFlabentries[17].get())
    e = arange(float(OFlabentries[14].get()), float(OFlabentries[15].get()), EPasso)
    resol = float(OFlabentries[13].get())
    dose = float(OFlabentries[24].get())
Matheus Müller's avatar
Matheus Müller committed
276
    X = e*0.
277
278
279
280
281
282
283
284
285
286
287
    gshift=0

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

    ## Calcula o espectro o elemento da amostra
   
    mt = Edict[0]['mass']
    c = Edict[0]['dist']
    alpha_lineshape = Edict[0]['LineShape']

288
    k = ((sqrt (mt**2-mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2	
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
    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, Edict[0]['profundidademax'], passox): 

        # Seção de choque - CSRf(Z1, Z2, Ein, M1, M2, Theta):
        CS = c[int(x/passox)]*k*CSRf(ion['Z'], Edict[0]['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:  
        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)
        X = besself + X
        linecolor='r'

    elif modelo == 'gaussiana':
        if LScontrol == 1:     
            sinalgaussiana = convoluiF0(sinalgaussiana, alpha_lineshape, EPasso)           
        X = sinalgaussiana + X
        linecolor='b'

    if CGausscontrol == 1:
        X, gshift = convoluiGauss(X,resol/2.35482,EPasso)    

    tempor.set_xdata(e-gshift)           
    tempor.set_ydata(X*dose)

Matheus Müller's avatar
Matheus Müller committed
336
337
    set_lim(e, X, gshift, dose)
    spect_file(e, X, gshift, dose)
338
     
339
    plt.show()
340
341
342

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