numerical.py 13.4 KB
Newer Older
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
8 9
from scipy.special import i1
from Ewindow import Edict
10
from layers import Ldict
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

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):
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
104
    X = e*0.
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)

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
126 127

        # Calcula a contribuição de cada camada
128
        for x in arange(passox, Edict[componente]['profundidademax'], passox): 
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
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'

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

176 177 178
    set_lim(e, X, gshift, dose)
    spect_file(e, X, gshift, dose)
    
179
    plt.show()
180 181

####################################################################################################################################################
182
####################################################################################################################################################
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())
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:
256 257 258 259
        X, gshift = convoluiGauss(X,resol/2.35482,EPasso)  
  
    tempor.set_xdata(e-gshift)           
    tempor.set_ydata(X*dose)
260

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

336 337
    set_lim(e, X, gshift, dose)
    spect_file(e, X, gshift, dose)
338
     
339
    plt.show()
340 341 342

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