Calculo por metodo de camadas operacional

numerical.py

@@ -7,6 +7,7 @@ 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
 
@@ -172,4 +173,100 @@ def spectro(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol):
 ####################################################################################################################################################
 
 def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol):
-    return 0
+
+    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())
+    X = e*0.
+    dose = 12e34
+    gshift=0.
+    antizerodiv=0.
+
+    Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180
+    mi = ion['mass']
+
+    for camada in Ldict:
+
+        ## Calcula o espectro de cada elemento da amostra
+        for componente in range(len(camada.Elements)):
+                          
+            mt = camada.Elements[componente]['mass']
+            c = float(camada.Spinboxes[componente].get())
+            alpha_lineshape = camada.Elements[componente]['LineShape']
+            dedx = float(camada.LandE[4].get())
+            dW2dx = float(camada.LandE[5].get())
+
+            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.))
+
+            sinalbessel = e*0
+            sinalgaussiana = e*0
+
+            xinit = float(camada.LandE[6].get())
+            xfinal = float(camada.LandE[7].get())
+            if  xinit == 0:
+                antizerodiv = passox
+            else:
+                antizerodiv = 0
+
+            # Calcula a contribuição de cada camada
+            for x in arange(xinit+antizerodiv, xfinal, passox): 
+
+                # 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:  
+                # 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':
+                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)    
+    plt.plot(e-gshift,X*dose, linecolor)
+    ylim([0,max(X*dose)*1.1])
+    suptitle('Signal received', fontsize=12)
+    xlabel("Enegy (eV)")
+    ylabel("Yield")
+
+    data2 = open('temp.sim', 'w')
+    for i in range(len(X)):
+        print>>data2, e[i]-gshift,X[i]*dose
+    data2.close()
+     
+    show()
+
+####################################################################################################################################################
+####################################################################################################################################################