numerical.py 6.14 KB
Newer Older
Matheus Müller's avatar
Matheus Müller committed
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
62
# 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

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

def spectro(p, modelo, ion, Emin, Emax, EPasso, passox, LScontrol, resol):
Matheus Müller's avatar
Matheus Müller committed
63
64

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

    elif modelo == 'gaussiana':
        X = e*0.
        try:
            XX
        except NameError: 
            XX, = plot(e,e*0., 'r')
        else:
            XX.set_xdata(e)

Matheus Müller's avatar
Matheus Müller committed
83
    dose = 2e34
Matheus Müller's avatar
Matheus Müller committed
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
136

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

    ## 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
137
            Y, gshiftY = convoluiGauss(Y,resol/2.35482,EPasso)
Matheus Müller's avatar
Matheus Müller committed
138
139
140
        elif modelo == 'gaussiana':
            if LScontrol == 1:     
                sinalgaussiana = convoluiF0(sinalgaussiana, alpha_lineshape, EPasso)           
141
142
143
144
145
146
147
148
149
150
            X = sinalgaussiana + X
            X, gshiftX = convoluiGauss(X,resol/2.35482,EPasso)

    if modelo == 'bessel':
        YY.set_ydata(Y*dose)
        ylim([0,max(Y*dose)*1.1])
    elif modelo == 'gaussiana':
        XX.set_ydata(X*dose)
        ylim([0,max(X*dose)*1.1])
    #plt.plot(e,Y*dose)
Matheus Müller's avatar
Matheus Müller committed
151
152
153
154
    suptitle('Signal received', fontsize=12)
    xlabel("Enegy (eV)")
    ylabel("Yield")

155
156
157
    #amostra = np.loadtxt('H2.py')
    #xxx= amostra[:,0]*1000
    #ccc=amostra[:,1]
Matheus Müller's avatar
Matheus Müller committed
158

159
    #plot (xxx, ccc, 'r.')
Matheus Müller's avatar
Matheus Müller committed
160
161
162
163
164

    show()

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