numerical.py 13.3 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
116
117
118
119

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

120
121
        sinalbessel = e*0
        sinalgaussiana = e*0
Matheus Müller's avatar
Matheus Müller committed
122
123

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

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

Matheus Müller's avatar
Matheus Müller committed
159
160
161
        elif modelo == 'gaussiana':
            if LScontrol == 1:     
                sinalgaussiana = convoluiF0(sinalgaussiana, alpha_lineshape, EPasso)           
162
            X = sinalgaussiana + X
163
164
165
166
            linecolor='b'

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

    tempor.set_xdata(e-gshift)           
    tempor.set_ydata(X*dose)
Matheus Müller's avatar
Matheus Müller committed
170

Matheus Müller's avatar
Matheus Müller committed
171
172
173
    set_lim(e, X, gshift, dose)
    spect_file(e, X, gshift, dose)
    
Matheus Müller's avatar
Matheus Müller committed
174
175
176
    show()

####################################################################################################################################################
177
####################################################################################################################################################
Matheus Müller's avatar
Matheus Müller committed
178

179
def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, tempor):
180
181
182
183
184

    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())
185
    dose = float(OFlabentries[24].get())
Matheus Müller's avatar
Matheus Müller committed
186
    X = e*0.
187
    gshift=0.
188
    curdepht=passox    #Avoid zero division
189
190
191
192
193
194

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

    for camada in Ldict:

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

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

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

            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.))
            dw2dxef=dW2dx*(k**2./cos(p['Theta_in']*PI/180.)+1./cos(p['Theta_out']*PI/180.))

            sinalbessel = e*0
            sinalgaussiana = e*0

213
214
            # Calcula a contribuição de cada passo
            for x in arange(curdepht, xfinal, passox): 
215
216
217
218
219
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

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

248
249
        curdepht = xfinal

250
    if CGausscontrol == 1:
Matheus Müller's avatar
Matheus Müller committed
251
252
253
254
        X, gshift = convoluiGauss(X,resol/2.35482,EPasso)  
  
    tempor.set_xdata(e-gshift)           
    tempor.set_ydata(X*dose)
255

Matheus Müller's avatar
Matheus Müller committed
256
257
258
    set_lim(e, X, gshift, dose)
    spect_file(e, X, gshift, dose)
  
259
260
261
262
    show()

####################################################################################################################################################
####################################################################################################################################################
263
264
265
266
267
268
269
270

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
271
    X = e*0.
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
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
    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']

    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, 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
331
332
    set_lim(e, X, gshift, dose)
    spect_file(e, X, gshift, dose)
333
334
335
336
337
     
    show()

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