calc.py 5.73 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
32
33
34
35

####################################################################################################################################################
# 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()

36
37
38
def set_lim(x, x2, y, shift):  
    ylim([0,y*1.1])          
    xlim([x*0.95,x2*1.05])
39
40
41
42

def variables(OFlabentries):
    EPasso = float(OFlabentries[16].get())
    passox = float(OFlabentries[17].get())
43
44
45
    Emin = float(OFlabentries[14].get())
    Emax = float(OFlabentries[15].get())
    e = arange(Emin, Emax, EPasso)
46
47
    resol = float(OFlabentries[13].get())
    dose = float(OFlabentries[24].get())
48
    return EPasso, passox, e, resol, dose, Emin, Emax
49

50
51
52
53
54
55
56
57
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)

58
59
60
####################################################################################################################################################
####################################################################################################################################################

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

63
    EPasso, passox, e, resol, dose, Emin, Emax = variables(OFlabentries)
64
    espectro.clear()
65
    espectro.eshift = float(emin)
66

67
68
69
    gshift=0.

    Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180.
70
71
72
73
74
    mi = ion['mass']

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

76
            mt = Edict[componente]['mass']
77
            c = Edict[componente]['dist'].data
78
79
80
81
            alpha_lineshape = Edict[componente]['LineShape']

            k = kinematic(mt, mi, Theta_s)

82
            dedxef=p['dedx']*(k/cos(p['Theta_in']*PI/180.) + 1./cos(p['Theta_out']*PI/180.))
83
84
            dw2dxef=p['dW2dx']*(k**2./cos(p['Theta_in']*PI/180.)+1./cos(p['Theta_out']*PI/180.))

85
86
            sinalbessel = e*0.
            sinalgaussiana = e*0.
87
88

            # Calcula a contribuição de cada camada
89
            for x in arange(0, Edict[componente]['profundidademax'], passox): 
90
 
91
92
93
94
95
96
97
98
99
100
101
102
                # 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':
103
                    alpha=2.*dedxef/dw2dxef
104
                    m=alpha*dedxef
105
106
107
108
                    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))) 
109
110
111
112
113

                # Gaussiana:
                elif modelo == 'gaussiana':
                    Em=k*p['E0']-x*dedxef
                    sigma=x*dw2dxef
114
115
116
117
118
119
                    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)
            
120
121
122
            if modelo == 'bessel':
                linecolor='r'

123
            elif modelo == 'gaussiana':         
124
125
126
127
128
129
                linecolor='b'

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

    espectro.multiply(dose)
130
131
    plote[indice].set_xdata(espectro.axisenergy())           
    plote[indice].set_ydata(espectro.spectro)
132
133
134

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

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

139
    labeler(plote,botoes,fig)
140
141
    fig.canvas.draw()
    fig.show()
142
143