spectro.py 3.57 KB
 Matheus Muller committed Aug 15, 2011 1 ``````from pylab import * `````` Matheus Müller committed Aug 22, 2011 2 3 ``````from scipy.special import i1 from Edict import * `````` Matheus Muller committed Aug 15, 2011 4 5 6 `````` PI=math.pi `````` Matheus Müller committed Sep 09, 2011 7 8 9 10 11 12 ``````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 `````` Matheus Müller committed Aug 22, 2011 13 ``````## funcoes usadas pela classe spectro `````` Matheus Muller committed Aug 15, 2011 14 15 16 17 ``````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) `````` Matheus Müller committed Sep 09, 2011 18 ``````def modbessel(e, alpha, m, x, kE0): `````` Matheus Muller committed Aug 15, 2011 19 20 `````` # Funcao que retorna a perda de energia de acordo com a funcao # modificada de bessel do primeiro tipo de ordem 1 `````` Matheus Müller committed Sep 09, 2011 21 22 23 24 `````` # 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)) `````` Matheus Muller committed Aug 15, 2011 25 26 27 28 29 30 31 32 `````` ## 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. `````` Matheus Müller committed Sep 09, 2011 33 `````` def addbessel(self, peso, alpha, m, x, kE0): `````` Matheus Muller committed Aug 15, 2011 34 35 `````` # Adiciona uma funcao modificada de bessel do primeiro tipo: # energia media Em, com variancia sigma com o peso dado `````` Matheus Müller committed Sep 09, 2011 36 37 `````` #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) `````` Matheus Muller committed Aug 15, 2011 38 39 40 41 42 43 44 45 `````` 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) def plot(self): # Plota o spectro `````` Matheus Müller committed Sep 09, 2011 46 `````` plt.plot(self.e,self.Y) `````` Matheus Müller committed Aug 22, 2011 47 48 49 `````` suptitle('Signal received', fontsize=12) xlabel("Enegy (eV)") ylabel("Yield") `````` Matheus Muller committed Aug 15, 2011 50 51 52 `````` show() def setparam (self, parametros): `````` Matheus Müller committed Aug 22, 2011 53 `````` # Define os parametros fixos do spectro `````` Matheus Muller committed Aug 15, 2011 54 55 56 57 58 `````` self.E0 = parametros['E0'] self.dedx = parametros['dedx'] self.dW2dx = parametros['dW2dx'] self.Theta_in = parametros['Theta_in'] self.Theta_out = parametros['Theta_out'] `````` Matheus Müller committed Aug 22, 2011 59 60 61 62 63 64 65 66 67 68 `````` self.FWHM0 = parametros['FWHM0'] self.sigma0 = pow(self.FWHM0/2.35,2.) def sweepelements(self, param, modelo, ion): for ATOM in SLOT: if SLOT[ATOM]['control']==0: pass else: passo = 0.01 c = SLOT[ATOM]['dist'] `````` Matheus Müller committed Sep 09, 2011 69 `````` x = arange(0,20,passo) `````` Matheus Müller committed Aug 22, 2011 70 71 72 73 74 `````` self.setparam(param) Theta_s = PI - (self.Theta_in + self.Theta_out)*PI/180 mt = SLOT[ATOM]['mass'] mi = ion['mass'] `````` Matheus Müller committed Sep 09, 2011 75 `````` k = ( (sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2 `````` Matheus Müller committed Aug 22, 2011 76 `````` `````` Matheus Müller committed Sep 09, 2011 77 78 79 80 81 `````` for i in [0.1]: #x: 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) `````` Matheus Müller committed Aug 22, 2011 82 `````` if modelo == 'gaussiana': `````` Matheus Müller committed Sep 09, 2011 83 84 85 86 87 88 89 `````` self.addgauss(Em, sigma, c[int(i/passo)]*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) `````` Matheus Muller committed Aug 15, 2011 90