spec.py 10.4 KB
Newer Older
1 2 3 4
######################################################################################
#####  Classe que define o vetor que representa um espectro de perda de energia  #####
######################################################################################

5
import scipy 
6 7 8 9
from pylab import *              ### Contem funcoes matematicas                    ###
from scipy.weave import *        ### Biblioteca para implementar codigo em C       ###

PI = math.pi
10 11 12

class spectro:
    def __init__(self, N, x):
13 14
        self.eshift = 0.                   ### Deslocamento em energia             ###
        self.mshift = 0                    ### MEIS shift                          ###
15 16
        self.Emin = 0.
        self.Emax = 0.
17 18 19
        self.stepsize = x                  ### Passo em energia                    ###
        self.spectro = zeros(N)            ### Representa as contagens do espectro ###
        self.size = N                      ### Tamanho do Vetor                    ###
20

21
    def norma(self): ### Retorna a area total do espectro ###
22
        return sum(self.spectro)*self.stepsize
23

24 25 26
    def axisenergy(self): ### Retorna o vetor com as energias analogas a self.spectro ###
        ae = arange(self.eshift,((self.size*self.stepsize)+self.eshift-1),self.stepsize)
        return ae
27

28
    def resize(self,M): ### Altera o tamanho do vetor para M ###
29 30 31 32 33 34 35
        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):
36
                del self.spectro[-1]
37 38
        self.size = M

39
    def clear(self): ### Limpa (zera) o espectro ### 
40 41 42
        self.spectro = zeros(self.size)
        self.eshift = 0.

43 44 45 46
    def reverse(self):
        temp = list()
        for j in range(self.size):
            temp.insert(j, self.spectro[int(self.size-j-1)])
47
        self.spectro = temp
48
    
49
    def multiply(self, factor): ### Multiplica o espectro por um fator ###
50 51 52
        for j in range(self.size):
            self.spectro[j] = self.spectro[j]*factor

53
    def soma(self, lista): ### Soma um vetor ou lista ao espectro ###
54 55 56 57 58
        if len(lista) > self.size:
            self.resize(len(lista))
        for j in range(len(lista)):
            self.spectro[j] = self.spectro[j] + lista[j] 

59
    def setnorma(self, norman): ### Ajusta o espectro para conter a area dada ### 
60 61 62 63
        normaold = self.norma()
        for n in range(self.size):
            self.spectro[n] = self.spectro[n]*nan_to_num(norman/normaold)

64
    def normalize(self): ### Ajusta o espectro para conter a area = 1 ###
65 66
        self.spectro = self.spectro/self.norma()

67
    def value(self, x): ### Retorna o valor ###
68 69
        return self.spectro[self.eshift + x/self.stepsize]

70
    def venergy(self, energy): ### Retorna o valor em energia ###
71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
        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()

95
    def firstmoment(self): ### Retorna o primeiro momento da distribuicao ###
96 97 98 99 100 101 102 103
        moment=0.;
        i=0;
        #for i in arange(0,self.size,1):
        #    moment = moment + (1.*i*self.stepsize)*self.spectro[i]
        while moment < self.norma()/2.:
            moment = moment +self.spectro[i]*self.stepsize
            i= i +1
        return   i*self.stepsize
104

105
    def secondmoment(self): ### Retorna o segundo momento da distribuicao ###
106 107 108
        moment=0;
        for i in arange(0,self.size,1):
             moment = moment + (1.*i*self.stepsize)**2*self.spectro[i]
109 110 111 112 113 114
        #a = list()
        #for k in range(self.size):
        #    a.insert(k, [self.axisenergy()[k],self.spectro[k]])
        #return  np.var(a)#
        return moment/self.norma()
        #print self.axisenergy()
115

116
    def thirdmoment(self): ### Retorna o terceiro momento da distribuicao ###
117 118 119 120 121
        moment=0;
        for i in arange(0,self.size,1):
            moment = moment + (1.*i*self.stepsize)**3*self.spectro[i]
        return  moment/self.norma()

122
    def ConvolveGauss(self, wg): ### Convolui o espectro com uma curva normal (gaussiana) ###
123 124 125 126
        normaantes = self.norma()
        stepsize = self.stepsize
        self.resize(self.size + int(math.ceil(8.*wg/stepsize)))
        datasize = int(self.size)
127
        data = tuple(self.spectro) ### Utiliza inline tools para rodar em C e obter um ganho em performance ###
128 129 130
        code = """
            py::tuple temp(datasize);
            double datadois[datasize];
131 132 133 134
            double expo[datasize];
            double wg2=double(wg)*double(wg);
            double wg4_step=4.0*double(wg)/stepsize;
            double temptemp[datasize];
135
            for (int k=0;k<datasize;k++)
136 137 138
                datadois[k] = data[datasize-k-1];
            for (int k=0;k<datasize;k++)
                expo[k] = exp(-1.0*pow(stepsize*(k-wg4_step),2.0)/(wg2*2.0));
139
            for(int i=0;i<datasize;i++){
140
                temptemp[i] = 0.0;
141
                for(int j=0;j<=i;j++)
142
                    temptemp[i] = double(temptemp[i]) + datadois[j]*expo[i-j];
143
                }
144 145
            for (int k=0;k<datasize;k++) 
                temp[k] = temptemp[datasize-k-1];
146 147 148 149
            return_val = temp;
            """
        self.spectro = list(inline_tools.inline(code,['wg','datasize','stepsize','data'],type_converters=converters.blitz,compiler = 'gcc'))
        self.setnorma(normaantes)
150
        self.eshift = float(self.eshift) + 4.*wg
151

152
    def ConvolveLorentz(self, wl): ### Convolui o espectro com uma curva lorentziana ###
153 154 155 156
        normaantes = self.norma()
        stepsize = self.stepsize
        self.resize(self.size + int(math.ceil(20.*wl/stepsize)))
        datasize = int(self.size)
157
        data = tuple(self.spectro) ### Utiliza inline tools para rodar em C e obter um ganho em performance ###
158 159 160 161 162
        code = """
            py::tuple temp(datasize);
            double datadois[datasize];
            for (int k=0;k<datasize;k++)
                datadois[k] = data[k];
163
            double wl2=double(wl)*double(wl), wl10_step=10.0*double(wl)/stepsize;
164 165 166 167 168 169 170 171 172
            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)
173
        self.eshift = float(self.eshift) + 10.*wl
174

175
    def ConvolveF0(self, alpha): ### Convolui o espectro com a forma de linha ###
176 177
        normaantes = self.norma()
        stepsize = self.stepsize
178 179 180
        datasize = int(self.size)
        data = tuple(self.spectro) ### Utiliza inline tools para rodar em C e obter um ganho em performance ###
        code = """ 
181 182
            py::tuple temp(datasize);
            double datadois[datasize];
183 184 185 186
            double expo[datasize];
            double temptemp[datasize];
            for (int k=0;k<datasize;k++) 
                datadois[k] = data[datasize-k-1];
187
            for (int k=0;k<datasize;k++)
188
                expo[k] = exp(-1.0*alpha*k*stepsize);
189
            for(int i=0;i<datasize;i++){
190
                temptemp[i] = 0.0;
191
                for(int j=0;j<=i;j++)
192
                    temptemp[i] = double(temptemp[i]) + 1.0*datadois[j]*alpha*stepsize*expo[i-j];
193
                }
194 195
            for (int k=0;k<datasize;k++) 
                temp[k] = temptemp[datasize-k-1];
196 197 198 199 200
            return_val = temp;
            """
        self.spectro = list(inline_tools.inline(code,['alpha','datasize','stepsize','data'],type_converters=converters.blitz,compiler = 'gcc'))
        self.setnorma(normaantes)

201
    def Convolve(self, spectwo): ### Convolui dois espectros ###
202 203 204 205 206 207 208 209 210
        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)
211
        datadois = tuple(spectwo.spectro) ### Utiliza inline tools para rodar em C e obter um ganho em performance ###
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
        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')) 
227
        self.eshift = float(self.eshift) + float(spectwo.eshift)
228 229
        self.setnorma(normaantes)

230
    def SetGauss(self, step, size, wg): ### Define o espectro como uma curva normal (gaussiana) ###
231 232 233 234 235 236 237
        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)
238
        self.eshift = float(self.eshift) - 4*wg
239

240
    def SetLorentz(self, step, size, wl): ### Define o espectro como uma curva lorentziana ###
241 242 243 244 245 246 247
        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)
248
        self.eshift = float(self.eshift) - 10*wl
249