Commit 6658267c authored by Matheus Müller's avatar Matheus Müller

Forma de linha e melhorias

parent 9c99ee57
from numpy import *
SLOT = dict(e0 = dict(name = "New", symbol = "Hf", mass = 178.00, Z = 72, dist = [1./3.]*200 + [0.]*1800, control = 1),
e1 = dict(name = "New", symbol = "O", mass = 16.00, Z = 8, dist = [2./3.]*200 + [0.]*1800, control = 1),
e2 = dict(name = "New", symbol = "Si", mass = 28.00, Z = 14, dist = [0.]*200 + [1.]*1800, control = 1),
e3 = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
e4 = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
e5 = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
e6 = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
e7 = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
e8 = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
e9 = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
ELEk = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
ELEl = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
ELEm = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
ELEn = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0),
ELEo = dict(name = "New", symbol = "-", mass = 0.00, Z = 0, dist = zeros((1000), float), control = 0))
SLOT = dict(e0 = dict(name = "New", symbol = "Hf", mass = 178.00, Z = 72., dist = [0.]*60 + [1./3.]*400 + [0.]*3600, control = 1),
e1 = dict(name = "New", symbol = "O", mass = 16.00, Z = 8., dist = [0.]*60 + [2./3.]*400 + [0.]*3600, control = 1),
e2 = dict(name = "New", symbol = "Si", mass = 28.00, Z = 14., dist = [0.]*460 + [1.]*3600, control = 1))
import Tkinter as tk
from PTE import *
from numpy import ones
Edict = list()
Edict = dict()
Ebuttons = list()
i=0
......@@ -9,20 +10,20 @@ points = []
spline = 0
tag1 = "theline"
#######################################
##############################################################################
# Element callback
def elem_callback(edict , ebutton, frm):
#plot_dist(edict, frm)
PTable(edict , ebutton)
#######################################
##############################################################################
# Element addition
def ewin_build(window):
def ewin_build(window, xmax, xstep):
def create():
global i, Edict, Ebuttons
Edict.insert(i, dict(name = "New", symbol = "Hf", mass = 178.00, Z = 72, dist = [1./3.]*200 + [0.]*1800, control = 0))
Edict[i] = dict(name = "New", symbol = "Hf", mass = 178.00, Z = 72, dist = ones(int(xmax/xstep)), control = 0, LineShape = 200.)
Ebuttons.insert(i, tk.Button(Eframe,text=i,width=1,height=1, command=lambda i=i : elem_callback(Edict[i] , Ebuttons[i], Dframe)))
Ebuttons[i].pack(side='left')
i=i+1
......@@ -31,7 +32,7 @@ def ewin_build(window):
but.pack()
but.place(anchor='n', y=5, x=44)
#######################################
##############################################################################
# Element Frames
Label1 = tk.LabelFrame(window, text = 'Elements', relief='raised', bd=2)
......@@ -47,7 +48,7 @@ def ewin_build(window):
Dframe = tk.Frame(Label2)
Dframe.pack()
#######################################
##############################################################################
# Canvas callback
def point(event):
Distcanvas.create_oval(event.x, event.y, event.x+1, event.y+1, fill="black")
......@@ -73,7 +74,7 @@ def ewin_build(window):
print event.x, event.y
return spline
#######################################
##############################################################################
# Drawing Canvas
Distcanvas = tk.Canvas(Dframe, bg="white", width=600, height= 300)
......@@ -85,9 +86,8 @@ def ewin_build(window):
Distcanvas.bind("<Button-3>", graph)
Distcanvas.bind("<Button-2>", toggle)
#######################################
##############################################################################
# Init
create()
print Distcanvas.grid_size()
#######################################
##############################################################################
1.5200e+02 4.1800e+02
1.7400e+02 4.2900e+02
1.9000e+02 4.4600e+02
1.9900e+02 4.5200e+02
2.0600e+02 4.5200e+02
2.1200e+02 4.4500e+02
2.1800e+02 4.3300e+02
2.2400e+02 4.2000e+02
2.3200e+02 4.1100e+02
2.4800e+02 4.0700e+02
2.7500e+02 4.0200e+02
2.9100e+02 4.0200e+02
3.0900e+02 4.0400e+02
3.2400e+02 4.0400e+02
3.3400e+02 4.0100e+02
3.4300e+02 3.9300e+02
3.4900e+02 3.8700e+02
3.5400e+02 3.7900e+02
3.6000e+02 3.6600e+02
3.6600e+02 3.5300e+02
3.7000e+02 3.4100e+02
3.7600e+02 3.2200e+02
3.8000e+02 3.0900e+02
3.8600e+02 2.9300e+02
3.9200e+02 2.8400e+02
3.9900e+02 2.7400e+02
4.1100e+02 2.7100e+02
4.2500e+02 2.7100e+02
4.4000e+02 2.7700e+02
4.5300e+02 2.8100e+02
4.6500e+02 2.7900e+02
4.7500e+02 2.7100e+02
4.8300e+02 2.5900e+02
4.8700e+02 2.4700e+02
4.9100e+02 2.3300e+02
4.9500e+02 2.1600e+02
4.9900e+02 1.9900e+02
5.0400e+02 1.7800e+02
5.0900e+02 1.5700e+02
5.1300e+02 1.4500e+02
5.1500e+02 1.3700e+02
5.1900e+02 1.3300e+02
5.2200e+02 1.3100e+02
5.2600e+02 1.3100e+02
5.3000e+02 1.3800e+02
5.3400e+02 1.4800e+02
5.3800e+02 1.6500e+02
5.4000e+02 1.8100e+02
5.4400e+02 2.0300e+02
5.4600e+02 2.2000e+02
5.5000e+02 2.5000e+02
5.5600e+02 2.9400e+02
5.5600e+02 3.1500e+02
5.6100e+02 3.3600e+02
5.6500e+02 3.4800e+02
5.6700e+02 3.5000e+02
5.6900e+02 3.4800e+02
5.7200e+02 3.4000e+02
5.7600e+02 3.2200e+02
5.7800e+02 3.0100e+02
5.8100e+02 2.8100e+02
5.8300e+02 2.5700e+02
5.8700e+02 2.3300e+02
5.8800e+02 2.1200e+02
5.9100e+02 1.9000e+02
5.9400e+02 1.6600e+02
5.9700e+02 1.3700e+02
6.0100e+02 1.1100e+02
6.0500e+02 9.1000e+01
6.0900e+02 7.9000e+01
6.1400e+02 7.2000e+01
6.2000e+02 7.1000e+01
6.2800e+02 7.0000e+01
......@@ -2,16 +2,15 @@
import Tkinter as tk
from Ewindow import *
from numpy import *
from spectro import *
from numerical import *
import sys
#######################################
##############################################################################
# Initial parameters
param = dict(dedx=192. ,Theta_out=70., dW2dx=20740., E0= 100000., Theta_in=0., FWHM0=180.)
param = dict(dedx=192. ,Theta_out=70., dW2dx=20740., E0= 100000., Theta_in=0.) #, FWHM0=180.)
ionb = dict(Z=1, mass=1.0079)
#######################################
##############################################################################
# Calculation callback
# Defines parameters
......@@ -22,23 +21,21 @@ def init_calc():
param['Theta_out'] = float(EntryT_out.get())
param['Theta_in'] = float(EntryT_in.get())
param['E0'] = float(EntryE0.get())
param['FWHM0'] = float(EntryFWHM0.get())
#param['FWHM0'] = float(EntryFWHM0.get())
ionb['mass'] = float(Entryionmass.get())
ionb['Z'] = float(EntryionZ.get())
data=spectro_espalhamento(float(EntryEmin.get()),float(EntryEmax.get()),20)
data.sweepelements(param, str(modelvar.get()), ionb)
data.plot()
spectro(param, str(modelvar.get()), ionb, float(EntryEmin.get()),float(EntryEmax.get()), float(EntryEstep.get()), float(EntryDstep.get()), int(LSvar.get()) )
#######################################
##############################################################################
# Main window
root = tk.Tk()
root.minsize(385,260)
root.minsize(400,290)
root.title('Open Flatus')
root.geometry('385x260+200+400')
root.geometry('400x290+200+400')
#######################################
##############################################################################
# Elements window
ewin = tk.Tk()
......@@ -46,7 +43,7 @@ ewin.minsize(600,400)
ewin.title('Elements')
ewin.geometry('600x400+600+400')
#######################################
##############################################################################
# Frames
MasterFrame1 = tk.Frame(root)
MasterFrame1.pack(side='left')
......@@ -63,56 +60,64 @@ Pframel.pack(side='left')
Pframee = tk.Frame(Label2)
Pframee.pack(side='right')
Iframel = tk.Frame(Label3, width=8)
Iframel = tk.Frame(Label3, width=9)
Iframel.pack(side='left')
Iframee = tk.Frame(Label3, width=12)
Iframee = tk.Frame(Label3, width=11)
Iframee.pack(side='right')
#######################################
##############################################################################
# Parameters
# Labels
Labeldedx = tk.Label(Pframel, width=8, text='dƐ/dx')
LabeldW2dx = tk.Label(Pframel, width=8, text='dω²/dx')
LabelT_out = tk.Label(Pframel, width=8, text='θ out')
LabelT_in = tk.Label(Pframel, width=8, text='θ in')
LabelE0 = tk.Label(Pframel, width=8, text='E0')
LabelFWHM0 = tk.Label(Pframel, width=8, text='FWHM0')
LabelEmin = tk.Label(Pframel, width=8, text='Emin')
LabelEmax = tk.Label(Pframel, width=8, text='Emax')
Labeldedx = tk.Label(Pframel, width=9, text='dƐ/dx')
LabeldW2dx = tk.Label(Pframel, width=9, text='dω²/dx')
LabelT_out = tk.Label(Pframel, width=9, text='θ out')
LabelT_in = tk.Label(Pframel, width=9, text='θ in')
#LabelFWHM0 = tk.Label(Pframel, width=9, text='FWHM0')
LabelEmin = tk.Label(Pframel, width=9, text='Emin')
LabelEmax = tk.Label(Pframel, width=9, text='Emax')
LabelEstep = tk.Label(Pframel, width=9, text='Estep')
LabelDstep = tk.Label(Pframel, width=9, text='Depht step')
Labeldedx.pack()
LabeldW2dx.pack()
LabelT_out.pack()
LabelT_in.pack()
LabelE0.pack()
LabelFWHM0.pack()
#LabelFWHM0.pack()
LabelEmin.pack()
LabelEmax.pack()
LabelEstep.pack()
LabelDstep.pack()
# Entries
Entrydedx = tk.Entry(Pframee, width=12)
EntrydW2dx = tk.Entry(Pframee, width=12)
EntryT_out = tk.Entry(Pframee, width=12)
EntryT_in = tk.Entry(Pframee, width=12)
EntryE0 = tk.Entry(Pframee, width=12)
EntryFWHM0 = tk.Entry(Pframee, width=12)
EntryEmin = tk.Entry(Pframee, width=12)
EntryEmax = tk.Entry(Pframee, width=12)
Entrydedx = tk.Entry(Pframee, width=11)
EntrydW2dx = tk.Entry(Pframee, width=11)
EntryT_out = tk.Entry(Pframee, width=11)
EntryT_in = tk.Entry(Pframee, width=11)
#EntryFWHM0 = tk.Entry(Pframee, width=11)
EntryEmin = tk.Entry(Pframee, width=11)
EntryEmax = tk.Entry(Pframee, width=11)
EntryEstep = tk.Entry(Pframee, width=11)
EntryDstep = tk.Entry(Pframee, width=11)
Entrydedx.pack()
EntrydW2dx.pack()
EntryT_out.pack()
EntryT_in.pack()
EntryE0.pack()
EntryFWHM0.pack()
#EntryFWHM0.pack()
EntryEmin.pack()
EntryEmax.pack()
EntryEstep.pack()
EntryDstep.pack()
#Ion
Labelionmass = tk.Label(Iframel, width=8, text='Mass')
LabelionZ = tk.Label(Iframel, width=8, text='Z')
LabelE0 = tk.Label(Iframel, width=9, text='Energy')
Labelionmass = tk.Label(Iframel, width=9, text='Mass')
LabelionZ = tk.Label(Iframel, width=9, text='Z')
LabelE0.pack()
Labelionmass.pack()
LabelionZ.pack()
Entryionmass = tk.Entry(Iframee, width=12)
EntryionZ = tk.Entry(Iframee, width=12)
EntryE0 = tk.Entry(Iframee, width=11)
Entryionmass = tk.Entry(Iframee, width=11)
EntryionZ = tk.Entry(Iframee, width=11)
EntryE0.pack()
Entryionmass.pack()
EntryionZ.pack()
......@@ -122,13 +127,15 @@ EntrydW2dx.insert(0,param['dW2dx'])
EntryT_out.insert(0,param['Theta_out'])
EntryT_in.insert(0,param['Theta_in'])
EntryE0.insert(0,param['E0'])
EntryFWHM0.insert(0,param['FWHM0'])
#EntryFWHM0.insert(0,param['FWHM0'])
EntryEmin.insert(0,70000)
EntryEmax.insert(0,100000)
EntryEstep.insert(0,20)
EntryDstep.insert(0,0.01)
Entryionmass.insert(0, ionb['mass'])
EntryionZ.insert(0, ionb['Z'])
#######################################
##############################################################################
# Energy loss models
modelvar = tk.StringVar()
......@@ -137,7 +144,14 @@ R1.pack()
R2 = tk.Radiobutton(Label4, text="Bessel", variable=modelvar, value='bessel')
R2.pack()
#######################################
##############################################################################
# Line shape checkbox
LSvar = tk.IntVar()
C1 = tk.Checkbutton(Label4, text="Use Line Shape", variable=LSvar)
C1.pack()
##############################################################################
CalcButton = tk.Button(Label4, text = 'Calculate', command = init_calc, bd=4, width=20, height=4)
CalcButton.pack()
......@@ -145,18 +159,16 @@ CalcButton.pack()
ExitButton = tk.Button(Label4, text = 'Exit', command = sys.exit, bd=4, width=20, height=2)
ExitButton.pack(side='bottom')
#######################################
##############################################################################
# Elements callback
ewin_build(ewin)
ewin_build(ewin, 15, float(EntryDstep.get()))
#######################################
##############################################################################
data=spectro_espalhamento(float(EntryEmin.get()),float(EntryEmax.get()),20)
data.sweepelements(param, 'bessel', ionb)
data.plot()
spectro(param, 'gaussiana', ionb, float(EntryEmin.get()), float(EntryEmax.get()), float(EntryEstep.get()), float(EntryDstep.get()), int(LSvar.get()) )
root.mainloop()
ewin.mainloop()
#######################################
##############################################################################
#-*- 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
####################################################################################################################################################
def spectro(p, modelo, ion, Emin, Emax, EPasso, passox, LScontrol):
e = arange(Emin, Emax, EPasso)
Y = e*0.
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
elif modelo == 'gaussiana':
if LScontrol == 1:
sinalgaussiana = convoluiF0(sinalgaussiana, alpha_lineshape, EPasso)
Y = sinalgaussiana + Y
plt.plot(e,Y)
suptitle('Signal received', fontsize=12)
xlabel("Enegy (eV)")
ylabel("Yield")
amostra = np.loadtxt('H.py')
xxx= (amostra[:,0]*15000/650)+85000
ccc=amostra[:,1]*15e-37 -1e-34
plot (xxx, ccc)
show()
####################################################################################################################################################
......@@ -15,31 +15,30 @@ def gaussian(e, Em, sigma):
# Funcao que retorna uma gaussiana centrada em Em e largura sigma
return exp((-(e-Em)**2)/(2.*sigma))/sqrt(2.*PI*sigma)
def modbessel(e, alpha, m, x, kE0):
#def modbessel(e, alpha, m, x, k, E0):
# Funcao que retorna a perda de energia de acordo com a funcao
# modificada de bessel do primeiro tipo de ordem 1
# N FUNCIONA!
lbd = m*alpha*x
#print lbd*exp((-m*x-alpha*(e-kE0))/1000)*i1(2*sqrt(lbd*(e-kE0)/1000))/(sqrt(lbd*(e-kE0)/1000 +.0000001))
return lbd*exp((-m*x-alpha*(e-kE0))/1000)*i1(2*sqrt(lbd*(e-kE0)/1000))/(sqrt(lbd*(e-kE0)/1000 +.0000001))
## Definindo classe referente a objetos de espectro de espalhamento
class spectro_espalhamento:
def __init__(self, Emin, Emax, EPasso):
# Inicia o spectro no limite dado, com o passo EPasso
self.e = arange(Emin, Emax, EPasso)
self.Y= e*0.
self.Y = e*0.
def addbessel(self, peso, alpha, m, x, kE0):
def addbessel(self, peso, alpha, m, x, k, E0):
# Adiciona uma funcao modificada de bessel do primeiro tipo:
# energia media Em, com variancia sigma com o peso dado
#self.Y=self.Y +peso*modbesselp(self.e, Em, sigma, dedx, dw2dx, x)
self.Y=self.Y +peso*modbessel(self.e, alpha, m, x, kE0)
lbd=m*x*alpha
bessel = lbd*exp(-m*x-alpha*(k*E0-e))*i1(2.*sqrt(lbd*(k*E0-e)))/(sqrt(lbd*(k*E0-e)))
self.Y = self.Y + peso*bessel #modbessel(self.e, alpha, m, x, k, E0)
print len(self.Y), len(self.e)
def addgauss(self,Em,sigma,peso):
# Adiciona uma gaussiana ao spectro:
# centrada em Em, com variancia sigma com o peso dado
self.Y=self.Y+peso*gaussian(self.e, Em, sigma)
self.Y = self.Y+peso*gaussian(self.e, Em, sigma)
def plot(self):
# Plota o spectro
......@@ -64,9 +63,9 @@ class spectro_espalhamento:
if SLOT[ATOM]['control']==0:
pass
else:
passo = 0.01
xstep = 0.01
c = SLOT[ATOM]['dist']
x = arange(0,20,passo)
x = arange(0.1,20,xstep)
self.setparam(param)
Theta_s = PI - (self.Theta_in + self.Theta_out)*PI/180
......@@ -74,17 +73,21 @@ class spectro_espalhamento:
mi = ion['mass']
k = ( (sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2
for i in x:
for i in [0.1]: #x:
if modelo == 'gaussiana':
sigma=i*self.dW2dx*(k**2./cos(self.Theta_in*PI/180.)+1./cos(self.Theta_out*PI/180.))+self.sigma0
Em = k*self.E0 - i*self.dedx*(k/cos(self.Theta_in*PI/180.) + 1/cos(self.Theta_out*PI/180.))
CS = CSRf(ion['Z'], SLOT[ATOM]['Z'], Em, ion['mass'], SLOT[ATOM]['mass'], Theta_s)
if modelo == 'gaussiana':
self.addgauss(Em, sigma, c[int(i/passo)]*k*CS)
self.addgauss(Em, sigma, c[int(i/xstep)]*k*CS)
elif modelo == 'bessel':
print Em
alpha = (2.* param['dedx']/10. * Em ) / ((param['dW2dx']/10.) *(sigma/i))
m = alpha * param['dedx']/10. * Em
self.addbessel(c[int(i/passo)]*k*CS, alpha, m, i,k*self.E0)
dedxef = self.dedx*(k/cos(self.Theta_in*PI/180.) + 1/cos(self.Theta_out*PI/180.))
dw2dxef = self.dW2dx*(k**2./cos(self.Theta_in*PI/180.)+1./cos(self.Theta_out*PI/180.))
alpha = dedxef*(2./dw2dxef)
m = alpha*dedxef
Em = k*self.E0-dedxef*x
CS = CSRf(ion['Z'], SLOT[ATOM]['Z'], Em , ion['mass'], SLOT[ATOM]['mass'], Theta_s)
self.addbessel(c[int(i/xstep)]*k*CS*xstep, alpha, m, i, k, self.E0)
from Ewindow import *
#-*- 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)
root = tk.Tk()
root.maxsize(800,600)
from pylab import *
from mpmath import *
ewin_build(root)
PI=math.pi
root.mainloop()
#########################################################
#### Parâmetros iniciais ####
# Feixe:
E0=100000 #eV
mi=1.
Zi=1.
# Alvo:
dedx=192. # eV/nm
dw2dx=20740. # eV^2/nm
Theta_in=0.
Theta_out=70.
Theta_s = PI - (Theta_in + Theta_out)*PI/180
espessura=2
# Cálculo:
deltaE=20000
passoE=20
passox=0.02
e = arange(E0-deltaE, E0, passoE)
besselcamada= arange(E0-deltaE, E0, passoE)*0
## 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=[[72.,178.49,217.],[14.,28.,121.],[8.,16.,102]] # 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
####################################
mt=178
Zt=72
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.))
alpha=dedxef*(2./dw2dxef)
m=alpha*dedxef
lbd = limit(lambda x: m*x*alpha, 0)
#bespecial = besseli(1, 2.*sqrt(lbd*(k*E0-e)) )
#i1(2.*sqrt(lbd*(k*E0-e)))
print lbd
for i in e:
besselcamada[int(i)] = lbd*exp(-lbd/alpha-alpha*(k*E0-int(i)))*besseli(1, 2.*sqrt(lbd*(k*E0-int(i))) ) /(sqrt(lbd*(k*E0-int(i))))
### Removendo Not A Number do vetor ###
ondeestaoNaN = isnan(besselcamada);
besselcamada[ondeestaoNaN] = 0
plot(e, besselcamada)
#####################################
#print totalbessel
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