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

Classes spectro e profile adicionadas, programa adaptado ao uso de spectro

parent 2fa6f27f
......@@ -2,14 +2,15 @@
import Tkinter as tk
from Ewindow import *
from numerical import *
from calc import *
from spec import *
import tkFileDialog
fig = plt.figure()
ax = fig.add_subplot(111)
fig.canvas.set_window_title("Signal")
xx = arange(0,1, 0.1)
temp, = plot(xx,xx*0)
espectro = spectro(100, 1.)
temp, = plot(espectro.spectro,espectro.spectro*0)
xlabel("Enegy (eV)")
ylabel("Yield")
doseold = 12e34
......@@ -61,12 +62,12 @@ def init_calc():
if int(distvar.get()) == 1:
for i in range(len(Edict)):
Edict[i]['dist'] = np.loadtxt(Edict[i]['symbol']+".prof")
spectro(param, str(modelvar.get()), ionb, OFlabentries, int(LSvar.get()), int(CGvar.get()), temp)
calc(param, str(modelvar.get()), ionb, OFlabentries, int(LSvar.get()), int(CGvar.get()), temp, 'draw', espectro)
elif int(distvar.get()) == 0:
spectrolayer(param, str(modelvar.get()), ionb, OFlabentries, int(LSvar.get()), int(CGvar.get()), temp)
calc(param, str(modelvar.get()), ionb, OFlabentries, int(LSvar.get()), int(CGvar.get()), temp, 'layer', espectro)
elif int(modevar.get()) == 1:
spectroRNA(param, str(modelvar.get()), ionb, OFlabentries, int(LSvar.get()), int(CGvar.get()), temp)
pass #spectroRNA(param, str(modelvar.get()), ionb, OFlabentries, int(LSvar.get()), int(CGvar.get()), temp)
##############################################################################
# Main window
......
#-*- 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
from scipy.special import i1
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:
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
####################################################################################################################################################
# Fator cinematico
def kinematic(mt, mi, esp):
return ((sqrt (mt**2-mi**2*(sin(esp)**2)) + mi*cos(esp) ) / (mi+mt) )**2
#k = mi*cos(Theta_s) + sqrt (mt**2 - mi**2*sqrt(sin(Theta_s)))
#k = k*k/(mi+mt)/(mi+mt)
####################################################################################################################################################
# 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])
def variables(OFlabentries):
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())
return EPasso, passox, e, resol, dose
####################################################################################################################################################
####################################################################################################################################################
def calc(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, tempor, metodo, espectro):
EPasso, passox, e, resol, dose = variables(OFlabentries)
espectro.clear()
gshift=0
Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180
mi = ion['mass']
if metodo == 'layer': pass
elif metodo == 'draw':
for componente in Edict:
mt = Edict[componente]['mass']
c = Edict[componente]['dist']
alpha_lineshape = Edict[componente]['LineShape']
k = kinematic(mt, mi, Theta_s)
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[componente]['profundidademax'], 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
# 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':
espectro.soma( list(nan_to_num(sinalbessel)) )
if LScontrol == 1:
espectro.ConvolveF0(alpha_lineshape)
linecolor='r'
elif modelo == 'gaussiana':
espectro.soma( list(nan_to_num(sinalgaussiana)) )
if LScontrol == 1:
espectro.ConvoluiF0(alpha_lineshape)
linecolor='b'
if CGausscontrol == 1:
espectro.ConvolveGauss(resol/2.35482)
#tempor.set_xdata(espectro.axisenergy())
#tempor.set_ydata(espectro.spectro*dose)
espectro.multiply(dose)
plot(espectro.axisenergy(), espectro.spectro)
#set_lim(espectro.axisenergy(), espectro.spectro, espectro.eshift, dose)
spect_file(espectro.axisenergy(), espectro.spectro, espectro.eshift, dose)
plt.show()
......@@ -113,7 +113,11 @@ def spectro(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, tempor):
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
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)
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.))
......@@ -204,7 +208,7 @@ def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, tempor)
dedx = float(camada.LandE[3].get())
dW2dx = float(camada.LandE[4].get())
k = ((sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2
k = ((sqrt (mt**2-mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2
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.))
......@@ -281,7 +285,7 @@ def spectroRNA(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, tempor):
c = Edict[0]['dist']
alpha_lineshape = Edict[0]['LineShape']
k = ((sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2
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.))
......
from pylab import *
import scipy
from scipy.weave import *
from scipy.special import *
from matplotlib import *
PI = math.pi
class profile:
def __init__(self, N, x):
self.shift = 0.
self.xmin = 0.
self.xmax = N*x
self.stepsize = x
self.data = zeros(N)
self.size = N
self.contam = 0.
def axisdepth(self):
return arange(self.contam, self.size*self.stepsize, self.stepsize)
def clear(self):
self.data = zeros(self.size)
def resize(self, size):
N = self.size
if M > N:
for j in range(M-N):
self.data.insert(N+j,0)
elif N > M:
for j in range(N-M):
self.data.pop()
self.size = M
self.xmax = M*self.stepsize
def setcontam(self, newcontam):
delta = newcontam-self.contam
dx = math.ceil(delta/self.stepsize)
if (delta < 0):
for i in arange(math.ceil(newcontam*self.stepsize),self.xmax):
self.data[i] = self.data[i-dx]
resize(self.size+delta)
elif (delta > 0):
self.resize(self.size + delta)
datatemp = profile(self.xmax,self.stepsize)
datatemp.data = self.data
for i in arange(math.floor(newcontam/self.stepsize),self.xmax):
self.data[i] = datatemp.data[i-dx]
contam = newcontam
for i in range(math.ceil(newcontam/self.stepsize)):
data[i]=0
def drawline(self, xa, xb, ca, cb):
contam = self.contam
stepsize = self.stepsize
if (xb >= self.xmax):
xb = self.xmax
if (xa < 0):
xa = 0.
if (xa >= self.xmax):
xa = self.xmax
if (xb < 0):
xb = 0
if (xa < xb):
for i in arange(math.floor(xa/stepsize), math.floor(xb/stepsize)):
if (i > math.floor(contam/stepsize)):
a = (cb - ca) / (xb - xa)
b = -a*xa + ca
conc = a*(1.*i*stepsize) + b
if (conc < 0):
conc = 0.
self.data[i] = conc;
if (xb < xa):
for i in arange(math.floor(xb/stepsize), math.floor(xa/stepsize)):
if (i > math.floor(contam/stepsize)):
a = (ca - cb) / (xa - xb)
b = -a*xb + cb
conc = a*(1.*i*stepsize) + b
if(conc < 0):
conc = 0
self.data[i] = conc
def adderfc(self, qtd, sigma):
xsigma = self.xmax
ratio = self.contam/self.stepsize
if (4.*sigma/self.stepsize < self.xmax):
xsigma = 4.*sigma/self.stepsize
for i in arange(math.ceil(ratio)+1, xsigma+ratio):
j = i - math.floor(ratio)
self.data[i] = self.data[i] + qtd*scipy.special.erfc(j*self.stepsize/sqrt(2.)/sigma);
def remerfc(self, qtd, sigma):
xsigma = self.xmax
ratio = self.contam/self.stepsize
if (4.*sigma/self.stepsize < xmax):
xsigma = 4.*sigma/self.stepsize
for i in arange(math.ceil(ratio)+1, xsigma+ratio):
j = i - math.floor(ratio)
self.data[i] = self.data[i] - qtd*scipy.special.erfc(j*self.stepsize/sqrt(2.)/sigma)
def erfc(self, x):
passo = .005
valor = 0.
for i in range(x/passo + 1):
valor += exp(-(i*passo)**2)
erf = 1. - 2.*valor*passo/sqrt(PI)
if erf<0:
erf = 0.
return erf
def quant(self):
soma = 0.
for i in range(self.size):
soma = soma + self.data[i]
soma = soma*self.stepsize
return soma
def Kn(self, n, m, tltcorr): #Nao compila
tltcorr = 1.
stepsize = self.stepsize
data = tuple(self.data)
xmax = self.xmax
datasize = int(self.size)
code = """
long double fatorial(int a)
{
int aa;
long double valor=1.0;
for(aa=1;aa<=a;aa++)
valor=valor*aa;
return valor;
}
double datab[datasize];
for (int k=0;k<datasize;k++)
datab[k] = data[k];
long double kk=0.0;
if(n==0){
for(int x=0;x<xmax;x++)
kk = kk + 1.0*exp(-1.0*M*x*stepsize*tltcorr)*datab[x];
kk = kk*pow(stepsize,1.0);}
else{
for(int x=0;x<xmax;x++)
kk = kk + 1.0*powl(1.0*x*tltcorr/(1.0/stepsize),n)*exp(-1.0*M*x*tltcorr/(1.0/stepsize))*datab[x];
kk = kk*1.0*(1.0*powl(1.0*M,n)/(1.0*fatorial(n)*1.0/stepsize));}
return_val = kk/stepsize*1000.0;}
"""
return inline_tools.inline(code,['n','m','tltcorr','stepsize','data','xmax', 'datasize'],type_converters=converters.blitz,compiler = 'gcc')
def smooth(self, qts):
for n in range(qts):
for i in arange(math.ceil(self.contam/self.stepsize)+4,self.xmax-6):
self.data[i] = (self.data[i-3]+self.data[i-2]+self.data[i-1]+self.data[i]+self.data[i+1]+self.data[i+2]+self.data[i+3])/7.
prof = profile(1000, 1)
prof.adderfc(50, 500)
plot(prof.axisdepth(),prof.data)
prof.drawline(450,600,10,200)
plot(prof.axisdepth(),prof.data)
print prof.quant()
prof.smooth(20)
plot(prof.axisdepth(),prof.data)
print prof.quant()
show()
import scipy
from matplotlib import *
from pylab import *
from scipy.weave import *
from scipy.signal import *
PI=math.pi
class spectro:
def __init__(self, N, x):
self.eshift = 0.
self.mshift = 0
self.Emin = 0.
self.Emax = 0.
self.stepsize = x
self.spectro = ones(N)
self.size = N
def norma(self):
return sum(self.spectro)/self.stepsize
def axisenergy(self):
return arange(self.eshift,(self.size/self.stepsize)+self.eshift,self.stepsize)
def resize(self,M):
N = self.size
self.spectro = list(self.spectro)
if M > N:
for j in range(M-N):
self.spectro.insert(N+j+1,0)
elif N > M:
for j in range(N-M):
self.spectro.pop()
self.size = M
def clear(self):
self.spectro = zeros(self.size)
self.eshift = 0.
def multiply(self, factor):
for j in range(self.size):
self.spectro[j] = self.spectro[j]*factor
def soma(self, lista):
if len(lista) > self.size:
self.resize(len(lista))
for j in range(len(lista)):
self.spectro[j] = self.spectro[j] + lista[j]
def setnorma(self, norman):
normaold = self.norma()
for n in range(self.size):
self.spectro[n] = self.spectro[n]*nan_to_num(norman/normaold)
def normalize(self):
self.spectro = self.spectro/self.norma()
def value(self, x):
return self.spectro[self.eshift + x/self.stepsize]
def venergy(self, energy):
i = int(math.floor((self.eshift - energy + self.mshift)/self.stepsize))
if i < 0 or i > self.size:
return 0
else:
return (self.spectro[i+1]-self.spectro[i])/(self.stepsize*(energy+self.eshift)+self.spectro[i]-i*self.spectro[i+1]+self.spectro[i])
def setprimary(self, N, step, Emin, Emax):
imin = math.ceil(Emin/step)
imax = math.floor(Emax/step)
self.clear();
self.resize(Npontos);
for i in arange(imin,imax,1):
self.Data[i]=1.*(1./(i+0.5)-1./(i-0.5))*step**(-2.)
if((imin-0.5)*step<Emin):
self.Data[imin]-=1.*(1./Emin-1./(step*(imin-0.5)))/step
if((imin-0.5)*step>Emin):
self.Data[imin-1]=1.*(1./(step*(imin-0.5))-1./Emin)/step
if(Emax<(imax+0.5)*step):
self.Data[imax]-=1.*(1./(step*(imax+0.5))-1./Emax)/step
if(Emax>(imax+0.5)*step):
self.Data[imax+1]=1.*(1./Emax-1./(step*(imax+0.5)))/step
self.normalize()
def firstmoment(self):
moment=0;
for i in arange(0,self.size,1):
moment = moment + (1.*i*self.stepsize)*self.spectro[i]
return moment/self.norma()
def secondmoment(self):
moment=0;
for i in arange(0,self.size,1):
moment = moment + (1.*i*self.stepsize)**2*self.spectro[i]
return moment/self.norma()
def thirdmoment(self):
moment=0;
for i in arange(0,self.size,1):
moment = moment + (1.*i*self.stepsize)**3*self.spectro[i]
return moment/self.norma()
def ConvolveGauss(self, wg):
normaantes = self.norma()
stepsize = self.stepsize
self.resize(self.size + int(math.ceil(8.*wg/stepsize)))
datasize = int(self.size)
data = tuple(self.spectro)
code = """
py::tuple temp(datasize);
double datadois[datasize];
for (int k=0;k<datasize;k++)
datadois[k] = data[k];
wg/=2.3548;
double wg2=wg*wg;
double wg4_step=4.0*wg/stepsize;
for(int i=0;i<datasize;i++){
temp[i] = 0.0;
for(int j=0;j<=i;j++)
temp[i] = double(temp[i]) + datadois[j]*exp(-1.0*pow(stepsize*(i-j-wg4_step),2.0)/wg2/2.0);
}
return_val = temp;
"""
self.spectro = list(inline_tools.inline(code,['wg','datasize','stepsize','data'],type_converters=converters.blitz,compiler = 'gcc'))
self.setnorma(normaantes)
self.eshift =+4.*wg
def ConvolveLorentz(self, wl):
normaantes = self.norma()
stepsize = self.stepsize
self.resize(self.size + int(math.ceil(20.*wl/stepsize)))
datasize = int(self.size)
data = tuple(self.spectro)
code = """
py::tuple temp(datasize);
double datadois[datasize];
for (int k=0;k<datasize;k++)
datadois[k] = data[k];
double wl2=wl*wl, wl10_step=10.0*wl/stepsize;
for(int i=0;i<datasize;i++){
temp[i] = 0.0;
for(int j=0;j<=i;j++)
temp[i] = double(temp[i]) + datadois[j]/float(pow(stepsize*(i-j-wl10_step),2.0)+wl2);
}
return_val = temp;
"""
self.spectro = list(inline_tools.inline(code,['wl','datasize','stepsize','data'],type_converters=converters.blitz,compiler = 'gcc'))
self.setnorma(normaantes)
self.eshift =+ 10.*wl
def ConvolveF0(self, alpha):
normaantes = self.norma()
datasize = int(self.size)
stepsize = self.stepsize
data = tuple(self.spectro)
code = """
py::tuple temp(datasize);
double datadois[datasize];
for (int k=0;k<datasize;k++)
datadois[k] = data[k];
for(int i=0;i<datasize;i++){
temp[i] = 0.0;
for(int j=0;j<=i;j++)
temp[i] = double(temp[i]) + 1.0*datadois[j]*alpha*stepsize*exp(-1.0*alpha*(1.0*stepsize*(i-j)));
}
return_val = temp;
"""
self.spectro = list(inline_tools.inline(code,['alpha','datasize','stepsize','data'],type_converters=converters.blitz,compiler = 'gcc'))
self.setnorma(normaantes)
def Convolve(self, spectwo):
normaantes = self.norma()
stepsize = self.stepsize
if spectwo.size <= self.size:
datasize = int(self.size)
spectwo.resize(datasize)
else:
datasize = int(spectwo.size)
self.resize(datasize)
data = tuple(self.spectro)
datadois = tuple(spectwo.spectro)
code = """
py::tuple temp(datasize);
double datab[datasize];
double datadoisb[datasize];
for (int k=0;k<datasize;k++)
datab[k] = data[k];
for (int l=0;l<datasize;l++)
datadoisb[l] = datadois[l];
for (int i=0;i<datasize;i++){
temp[i] = 0.0;
for(int j=0;j<=i;j++)
temp[i] = double(temp[i]) + datab[j]*datadoisb[i-j]*stepsize;}
return_val = temp;
"""
self.spectro = list(inline_tools.inline(code,['datasize','stepsize','data','datadois'],type_converters=converters.blitz,compiler = 'gcc'))
self.eshift = self.eshift + spectwo.eshift
self.setnorma(normaantes)
def SetGauss(self, step, size, wg):
self.clear()
self.resize(size)
wg2 = wg**2.
wg4_step = 4.*wg/step
for i in range(size):
self.spectro[i] = exp(-2*(step*(i-wg4_step)**2)/wg2)
self.spectro = list(self.spectro)
self.eshift = -4*wg
def SetLorentz(self, step, size, wl):
self.clear()
self.resize(size)
wl2 = wl**2
wl10_step = 10*wl/step
for i in range(size):
self.spectro[i] = 1/(4*(step*(i-wl10_step)**2)+wl2)
self.spectro = list(self.spectro)
self.eshift = -10*wl
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