Commit c71439ae authored by Matheus Müller's avatar Matheus Müller

Renomeado para profiler, agora abre distribuicao salva no arquivo caso exista

parent d84799e0
#-*- 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.
# Incluída convolução com forma de linha originária da colisão frontal
#
# Referência: R. P. Pezzi, et al. Applied Physics Letters, 92 164102 (2008)
from pylab import *
from scipy.special import i1
PI=math.pi
######## Definindo a Seção de Choque #######
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
#########################################################
######## 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)
print len(eixo)
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
# return temp/temp.max(), wg4_step*passoE # Incluido para teste com Gaussiana -> Altura = 1
#########################################################
#### Parâmetros iniciais ####
# Feixe:
E0=99000 #eV
mi=1.
Zi=1.
# Alvo:
dedx=215. # eV/nm
dw2dx=20000. # eV^2/nm
Theta_in=7.
Theta_out=58.67
Theta_s = PI - (Theta_in + Theta_out)*PI/180
espessura=19
# Cálculo:
deltaE=20000
passoE=20
passox=0.01
e = arange(E0-deltaE, E0, passoE)
totalbessel=e*0
totalbessel_ls=e*0
totalgaussiana=e*0
resolucao = 250 # Largura a meia algura da Gaussiana que representa a resolução experimental (FWHM) em eV
## TODO: - Otimizar os cálculos de cada componente: evitar cálculos (de zeros) em energias desnecessários
## - Definir o intervalo de plotagem do espectro independente do deltaE
#############################
## Define a composição da amostra ##
elementos=[[57.,138.9,192.,"la.prof"],[38.,87.62,165.,"sr.prof"],[22.,47.867,142.,"sr.prof"]] # La, Sr,ti
#elementos=[[72.,178.49,217.,"hf.txt"],[14.,28.,121.,"si.txt"],[8.,16.,102,"o.txt"]] # Hf, Si e O
#elementos=[[72.,178.49,17.],[14.,28.,21.],[8.,16.,12]] # Hf, Si e O
#elementos=[[72.,178.49,217.]]#,[14.,28.,121.],[8.,16.,102]] # Hf, Si e O
# Esta lista deve incluir a distribuição
# em profundidade do elemento
####################################
## Calcula o espectro de cada elemento da amostra
for componente in elementos:
mt=componente[1];
Zt=componente[0];
alpha_lineshape=componente[2]
perfil = loadtxt(componente[3], unpack=True)
# plot(perfil)
k = ((sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2
dedxef=dedx*(k/cos(Theta_in*PI/180.) + 1/cos(Theta_out*PI/180.))
dw2dxef=dw2dx*(k**2./cos(Theta_in*PI/180.)+1./cos(Theta_out*PI/180.))
sinalbessel=e*0
sinalgaussiana=e*0
# Calcula a contribuição de cada camada
for x in arange(passox,espessura,passox):
# 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)
# Este truncamento resulta de valores NaN (Not a Number) que aparecem no vetor. Foi incluida uma
# operação para remover os NaN do vetor de Bessel
#
# CSRf(Z1, Z2, Ein, M1, M2, Theta):
secchoque=CSRf(Zi, Zt,k*E0-dedxef*x,mi,mt, Theta_s)*passox*perfil[int(x/passox)]
# É multiplicada pelo passo para manter a escala do gráfico
# independente do mesmo.
alpha=dedxef*(2./dw2dxef)
m=alpha*dedxef
lbd=m*x*alpha
besselcamada = lbd*exp(-m*x-alpha*(k*E0-e))*i1(2.*sqrt(lbd*(k*E0-e)))/(sqrt(lbd*(k*E0-e)))
### Removendo Not A Number do vetor ###
ondeestaoNaN = isnan(besselcamada);
besselcamada[ondeestaoNaN] = 0
### Incluir o sinal correspondente à ausência de colisões (Delta de Dirac em k*E0)
# Pode ser feito adicionando uma gaussiana de largura muito pequeno em k*E0 para circundar o
# fato de kE0 poder não cair exatamente em uma posição discreta do vetor
sinalbessel=besselcamada*secchoque+sinalbessel
## Plotar o sinal da camada
#plt.plot(e,besselcamada)
# Gaussiana:
Em=k*E0-x*dedxef
sigma=(sqrt(x*dw2dxef)+resolucao/2.35482)**2
gaussianacamada = exp((-(e-Em)**2)/(2.*sigma))/sqrt(2.*PI*sigma)
sinalgaussiana=gaussianacamada*secchoque+sinalgaussiana
#plt.plot(e,gaussianacamada)
# plt.plot(e, sinalbessel)
# plt.plot(e, sinalgaussiana)
### Convoluir com forma de linha da colisão frontal:
sinalbessel_ls=convoluiF0(sinalbessel,alpha_lineshape,passoE)
totalbessel=sinalbessel+totalbessel
totalbessel_ls=sinalbessel_ls+totalbessel_ls
totalbessel_ls_gauss,gshift=convoluiGauss(totalbessel_ls,resolucao/2.35482,passoE)
totalgaussiana=sinalgaussiana+totalgaussiana
#####################################
# Mostra o Gráfico
#bes=plt.plot(e,totalbessel,lw=2,color="red")
#besls=plt.plot(e,totalbessel_ls,lw=2,color="blue")
beslsgaus=plt.plot(e-gshift,totalbessel_ls_gauss,lw=2,color="orange")
gaus=plt.plot(e,totalgaussiana,lw=2,color="green")
#plt.plot(e,totalbessel*10,lw=2, color="blue")
#plt.plot(e,totalgaussiana*10,lw=2, color="green")
#plt.legend([bes,besls,beslsgaus,gaus],["Bessel","Bessel+lineshape","Bessel+ls+expres","Gaussiana"],loc= "upper left", shadow=True)
plt.legend([beslsgaus,gaus],["Bessel+ls+expres","Gaussiana"],loc= "upper left", shadow=True)
#plt.legend([bes,besls,gaus],["Bessel","Bessel+lineshape","Gaussiana"],loc= "upper left", shadow=True)
plt.xlabel("Energia dos Ions (eV)")
plt.ylabel("Rendimento (contagens)")
################################################
# Testando convolução com deltas de Dirac
#figure()
#xxy=e*0
#xxy[00]=1. # Correspondente a uma delta em 80000
#xxy[100]=1. # Correspondente a uma delta em 82000
#xxy[300]=1.
#xxy[400]=1.
#xxy[500]=1.
# Convoluindo as deltas com uma gaussiana de FWHM = 550 eV
#xxyz,g_shift=convoluiGauss(xxy,resolucao/2.35482,passoE)
# Convoluindo a convolução novamente -> FWHM Res = 777.81 eV = sqrt(550^2+550^2)
#xxyz2,g_shift2=convoluiGauss(xxyz,resolucao/2.35482,passoE)
#deltas = plot(e,xxy,lw=2)
#convo1 = plot(e-g_shift,xxyz+0.001,lw=2)
#convo2 = plot(e-g_shift-g_shift2,xxyz2+0.002,lw=2)
#convo2_analitica=plot(e,exp(-(e-82000)**2.0/(2.*(777.81/2.35482)**2))+0.002,lw=2)
#plt.legend([deltas,convo1,convo2,convo2_analitica],["Deltas","Conv Gauss 1x 550 eV","Conv Gauss 2x 550 eV","Conv Analitica (777.81)"],loc= "upper right", shadow=True)
#print e-g_shift, xxyz
#print str(g_shift) + " <- g_shift"
#print len(e)
#print len(xxy)
#print len(xxyz)
show()
##################
......@@ -6,17 +6,23 @@ from scipy.interpolate import interp1d
from scipy import interpolate
from pylab import *
from concentracao import *
import sys
listay = list()
listax = list()
drs = []
fig = plt.figure()
ax = fig.add_subplot(111)
filename = "data.txt" #raw_input()
filename2 = sys.argv[1]+".prof" #raw_input()
filename = sys.argv[1]+".dat" #raw_input()
print "nome dos arquivos "+filename+" e "+ filename2
fig.canvas.set_window_title(filename)
dxstep = 0.01
profundidademax = 20.
xx = arange(0,20, dxstep) # we'll create an x-axis from 0 to 2 pi
temp, = plot(xx,xx*0) # this is our initial plot, and does nothing
#modo = 'step'
modo = 'reta'
#######################################################################################################
# Add square
......@@ -28,12 +34,12 @@ def add_square(x, y):
dr = DraggableRectangle(rec)
dr.connect()
drs.append(dr)
rec.set_y(y-0.01)
rec.set_y(y)#-0.01)
listay.append(dr.rect.xy[1])
listax.append(dr.rect.xy[0])
dr.updateplot()
ylim([0,1])
xlim([0,20])
xlim([0,profundidademax])
def add_immovable_square(x, y):
global listay, drs, ax
......@@ -41,11 +47,11 @@ def add_immovable_square(x, y):
for rec in rect:
dr = DraggableRectangle(rec)
drs.append(dr)
rec.set_y(y-0.01)
rec.set_y(y)#-0.01)
listay.append(dr.rect.xy[1])
listax.append(dr.rect.xy[0])
ylim([0,1])
xlim([0,20])
xlim([0,profundidademax])
def remove_square():
global listay, drs
......@@ -53,6 +59,10 @@ def remove_square():
drs.pop()
print listay, drs
def clear_figure():
global fig
fig.clf()
#######################################################################################################
class voidfunctions:
......@@ -65,22 +75,30 @@ class voidfunctions:
'key_release_event', self.on_a)
self.cidr = self.rect.figure.canvas.mpl_connect(
'key_release_event', self.on_r)
self.cidc = self.rect.figure.canvas.mpl_connect(
'key_release_event', self.on_c)
def on_a(self, event):
if event.key=='a':
if not event.inaxes: return
add_square(float(event.xdata), float(event.ydata))
# NAO FUNCIONA (DELETA TUDO)
# NAO FUNCIONA
def on_r(self, event):
if event.key=='r':
if not event.inaxes: return
remove_square()
def on_c(self, event):
if event.key=='c':
if not event.inaxes: return
clear_figure()
def disconnect(self):
'disconnect all the stored connection ids'
self.rect.figure.canvas.mpl_disconnect(self.cida)
self.rect.figure.canvas.mpl_disconnect(self.cidr)
self.rect.figure.canvas.mpl_disconnect(self.cidc)
#######################################################################################################
......@@ -92,7 +110,7 @@ class DraggableRectangle:
self.background = None
def updateplot(self):
global listax, listay, temp, drs, dxstep
global listax, listay, temp, drs, dxstep, modo
xytemp = list(drs)
for i in range(len(drs)):
xytemp[i] = drs[i].rect.xy
......@@ -107,12 +125,17 @@ class DraggableRectangle:
data.close()
xydata = np.loadtxt(filename)
modo = 'step'
modo2 = 'reta'
xy = montaperfil( xydata, dxstep, modo)
data2 = open(filename2 , 'w')
for i in range(len(xy)):
print>>data2, xy[i]
data2.close()
temp.set_ydata(montaperfil( xydata, dxstep, modo)) # update the plot data
draw() # redraw the canvas
ylim([0,1])
xlim([0,20])
xlim([0,profundidademax])
def connect(self):
'connect to all the events we need'
......@@ -199,6 +222,14 @@ for void in voidrect:
add_immovable_square(0,0)
add_immovable_square(20,0)
try:
f = np.loadtxt(filename)
except IOError:
pass
else:
for initial in range( len(f) ):
add_square(f[initial][0],f[initial][1])
#######################################################################################################
plt.show()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment