spectro.py 5.31 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
from pylab import *

PI=math.pi

## Primeiro defino algumas funcoes que sao usadas pela classe spectro
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 modbesselp(e, param,x):
    # Funcao que retorna a perda de energia de acordo com a funcao 
    #   modificada de bessel do primeiro tipo de ordem 1
    # **** NAO FUNCIONA!! ****
    alpha=2.*param['dedx']/param['dW2dx']
    m = alpha * param['dedx']
    lbd=m*x*alpha
    return (lbd*exp(-m*x-alpha*e)/(sqrt(lbd*e)))*i1(2.*sqrt(lbd*e))


def profile(pares, passo):
    # Funcao que retorna os vetores espessura e perfil de concentracao
    #  a partir dos pares de pontos fornecidos.
    # pares deve ser uma matriz de duas colunas: [espessura, concentracao] 
    #  Modos: degraus - distribuicao do tipo escada. Ex. [10 , 1]
    #         representa uma amostra com 10 angstroms de espessura
    #         com concentracao constante    
    depth = sum(pares[:,0]) # A espessura total eh a soma da primeira coluna
    x = arange(0,depth, passo) # Monta o eixo X
    c = [] # Inicia o vetor concentracao vazio
    for i in range(len(pares)): # popula o vetor concentracao
       temp = ones(pares[i][0]/passo)*pares[i][1]
       c.extend(temp.tolist())
    return x,c



## 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.

    def addbessel(self,param,peso,x):
        # 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, param, x) 

    
    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
        #  Ainda falta colocar os nomes nos eixos e outras descricoes
        # 
        plot(self.e,self.Y)
        show()

    def setparam (self, parametros):
        #Define os parametros fixos do spectro
        # parametros eh um dicionario definido com algo como
        # param = dict(dedx=287. , dW2dx=35000 , E0= 100000, Theta_in=0, Theta_out=70, fwhm0=300 )
        self.E0 = parametros['E0']
        self.dedx = parametros['dedx']
        self.dW2dx = parametros['dW2dx']
        self.Theta_in = parametros['Theta_in']
        self.Theta_out = parametros['Theta_out']
        self.fwhm0 = parametros['fwhm0']
        self.sigma0 = pow(self.fwhm0/2.35,2.)
       # self. = parametros['']

    def addelements(self,arquivos, param,k,modelo):
        # Adiciona ao spectro os elementos definidos nos arquivos com 
        #  os parametros dados 
        #  *** Eh preciso incluir o fator cinematico apropriadamente ***
        #  Argumentos para addelements sao:
        #       Lista de arquivos do tipo ['file1.dat','file2.dat']
        #       Parametros
        #       Modelo de perda de energia: 'gaussiana' ou 'bessel'

        # ** Falta incluir informacoes de massa do elemento no arquivo
        # ** quem sabe na primeira linha incluir numero atomico e massa
        # ** e ao carregar o arquivo com np.loadtxt(), retirar estes valores
        # ** dos vetores x e c

        for arquivo in arquivos:
            amostra = np.loadtxt(arquivo)
            x=amostra[:,0]*0.01
            c=amostra[:,1]*1.
            self.setparam(param)
            for i in x:
               sigma=i*self.dW2dx*(k*k/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.))
               if modelo == 'gaussiana':
                   self.addgauss(Em, sigma, c[i*100]*k)

               if modelo == 'bessel':
                   self.addbessel(param, c[i*100]*k,x)
                   
            k=k-0.21


    def addelements2(self,arquivos, param,k,modelo,passo):
        # Adiciona ao spectro os elementos definidos nos arquivos com 
        #  os parametros dados 
        #  *** Eh preciso incluir o fator cinematico apropriadamente ***
        #  Argumentos para addelements sao:
        #       Lista de arquivos do tipo ['file1.dat','file2.dat']
        #       Parametros
        #       Modelo de perda de energia: 'gaussiana' ou 'bessel'

        # ** Falta incluir informacoes de massa do elemento no arquivo
        # ** quem sabe na primeira linha incluir numero atomico e massa
        # ** e ao carregar o arquivo com np.loadtxt(), retirar estes valores
        # ** dos vetores x e c

        for arquivo in arquivos:
            amostra = np.loadtxt(arquivo)
            x,c = profile(amostra,passo)
            self.setparam(param)
            for i in x:
               sigma=i*self.dW2dx*(k*k/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.))
               if modelo == 'gaussiana':
                   self.addgauss(Em, sigma, c[int(i/passo)]*k)

               if modelo == 'bessel':
                   self.addbessel(param, c[i*100]*k,x)
                   
            k=k-0.21