1 """
2 Spectra
3 """
4
5
6 import numpy as np
7 import pyfits as pf
8 import math
9 import glob
10 import lib
11 import copy as pythoncopy
12 import tem
13 import scipy.signal as spsig
14 import sys
15
17 """a spectra is a flux assiociated to wavelengths, the other parameters can be readfrom the SDSS data """
18
19 - def __init__(self, fluxes_cs = np.array([]), fluxes=np.array([]), noise =np.array([]), zsdss = 0, coeff0 = 3.5793, coeff1 = 0.0001, ra = 0, dec = 0, angst =np.array([])):
20 self.angst = angst
21 self.fluxes = fluxes
22 self.fluxes_cs = fluxes_cs
23 self.fluxes_csc = np.array([])
24 self.noise = noise
25 self.zsdss = zsdss
26 self.secz = 0
27 self.mask = np.array([])
28 self.plotcolour = "black"
29 self.coeff0 = coeff0
30 self.coeff1 = coeff1
31 self.ra = ra
32 self.dec = dec
33 self.name = 'not from fit file'
34 self.allfluxes = {'fluxes' : self.fluxes, 'fluxes_cs' : self.fluxes_cs, 'fluxes_csc' : self.fluxes_csc}
35
37 """initialize the mask to an array containing only 'False' values, this fonction should be called only once the spectra is non-empty"""
38 self.mask = np.ones(len(self.angst))
39 self.mask = np.array(self.mask != 1)
40
42 """Returns a deep copy of itself"""
43 return pythoncopy.deepcopy(self)
44
46 """converts a pixel to a wavelenght in angstroems, using a formula from SDSS site"""
47 return 10**(self.coeff0 + self.coeff1 * index)
48
50 """converts a wavelength in angstroems to a pixel (it is the inverse function 'pixtoangst')"""
51 return math.floor((math.log(wavelength,10) - self.coeff0)/self.coeff1)
52
54 """converts a pixel to a redshift"""
55 return 10**(self.coeff1*index)-1
56
58 """converts a redshift to a pixel (it is the inverse function 'pixtoz')"""
59 return int(math.log((z + 1), 10)/self.coeff1)
60
61 - def setangst(self, coeff0 = 3.5793, coeff1 = 0.0001):
62 """defines the wavelength of the spectra"""
63 i =np.arange(1, len(self.fluxes_cs)+1)
64 self.angst =np.array( map(self.pixtoangst,i) )
65
67 """sets the wavelengths corrsponding to the redshift"""
68 self.angst =(1 + z) / (1 + self.zsdss) * self.angst
69
71 """substracts the continuum of the spectra, using a median filter"""
72 kernel = 31
73 smooth = spsig.medfilt(self.fluxes, kernel)
74 self.fluxes_csc = self.fluxes - smooth
75 self.allfluxes['fluxes_csc'] = self.fluxes_csc
76
77
78 - def read(self, directory, filename):
79 """creates a spectra from a SDSS file (i.e. reads the redshift, the flux, the noise...etc."""
80 hdulist = pf.open(filename)
81 name = filename[len(directory):len(filename)]
82 fluxes = np.array(hdulist[0].data[0])
83 fluxes_cs =np.array( hdulist[0].data[1])
84 noise =np.array( hdulist[0].data[2])
85 coeff0 = hdulist[0].header['COEFF0']
86 coeff1 = hdulist[0].header['COEFF1']
87 ra = hdulist[0].header['RA']
88 dec = hdulist[0].header['DEC']
89 zsdss = hdulist[0].header['Z']
90 hdulist.close()
91
92 sp = spectra(fluxes_cs, fluxes, noise, zsdss, coeff0, coeff1, ra, dec)
93 sp.name = name
94 sp.setangst(coeff0, coeff1)
95 sp.calc_csc()
96
97 return sp
98
99
101 """for a spectra, we give a list of tuples (each one containing a wavelength and a width) and we hide the pixels at the given wavelength"""
102 if tupleslist != [] and tupleslist != None:
103 for tuple in tupleslist:
104 n = self.angsttopix((1 + z) * tuple[0])
105 i1 = self.angsttopix((1 + z) * tuple[0])-self.angsttopix((1 + z) * (tuple[0]-tuple[1]/2.0))
106 i2 = self.angsttopix((1 + z) * (tuple[0]+tuple[1]/2.0)) - self.angsttopix((1 + z) * tuple[0])
107 self.mask[n - i1:n + i2 + 1] = True
108
109
111 """hides the peaks higher than a given value [not used]"""
112 a = self.allfluxes[fluxchoice] < height
113 self.mask[a] = True
114
115 - def mask_borders(self, lambda1 = 5007 + 20, lambda2 = 6562.81 - 160):
116 """eliminates the flux outside the two given wavelengths"""
117 if lambda2 < lambda1:
118 x = lambda1
119 lambda1 = lambda2
120 lambda2 = x
121
122 lambda1 = lambda1 * (self.zsdss + 1)
123 lambda2 = lambda2 * (self.zsdss + 1)
124
125 n1 = self.angsttopix(lambda1)
126 n2 = self.angsttopix(lambda2)
127 size = len(self.fluxes_cs)
128
129 if n1 > 0:
130 self.mask[0:n1] = True
131 if n2 < size:
132 self.mask[n2:] = True
133
134
135 - def specsim(self, z1, z2, model, noiseintensity, secintensity = 1):
136 """simulates a spectra containing one main and one secondary spectra : you give one spectra as a model, the noise intensity, the two redshift and the ratio between the two spectra flux [not used]"""
137 self.zsdss = model.zsdss
138 self.secz = z2
139 self.fluxes_cs = model.fluxes_cs
140 self.noise = model.noise
141 self.setangst()
142 self.setz(z1)
143 self.zsdss = z1
144
145 size = len(self.fluxes_cs)
146 shifting = self.ztopix(z2)-self.ztopix(z1)
147 b = np.zeros((1,size))[0]
148 b[shifting:size-1] = self.fluxes_cs[0:size-shifting-1]
149 b *= secintensity
150 self.fluxes_cs = self.fluxes_cs + b
151 noise = noiseintensity*(np.random.rand(size)-0.5)
152 self.fluxes_cs = self.fluxes_cs + noise
153 self.calc_csc()
154
155
157 """This is a crosscorrelation spectra : it contains the origin spectra, a template which we use for the crosscorrelation and the result"""
162
164 """makes the crosscorrelation between the spectra and the templates : it starts when both arrays are 'face to face', then it shifts pixel by pixel until the last one """
165 corr = []
166
167 model = np.array(map(self.template.flux_gauss, self.spectra.angst))
168 a = self.spectra.allfluxes[fluxchoice].copy()
169 a[self.spectra.mask] = 0
170 index1 = len(a)
171 index2 = len(model)
172 """
173 #this part was made to avoid useless calculations, but with standards parameters, it didn't make it quicker, so I abandonned it. Maybe it can improved to be effective.
174 i = 0
175 mod = model[0:index2-1-i]
176 while np.any(mod > 0.01) :
177 corr = np.append(corr, np.sum(self.spectra.fluxes_cs[i:index1-1]*mod))
178 i += 1
179 mod = model[0:index2-1-i]
180 #print list(mod)
181 #sys.exit()
182 self.corr = corr
183 """
184
185 for i in range(index1):
186 corr.append(np.sum(a[i:index1-1]*model[0:index2-1-i]))
187 self.corr = np.array(corr)
188
189
191 """factory function for a list of spectra : if you give the name of a directory, it will return a list of spectra obtained from every fits files in the directory [not used]"""
192
193 spectralist = []
194
195 filelist = glob.glob(directory + '*.fit')
196 for file in filelist:
197 sp = spectra()
198 sp = sp.read(directory, file)
199
200 spectralist.append(sp)
201
202 return spectralist
203
204
205 -def cspfactory(spectralist, template, fluxchoice):
206 """factory function for a list of crosscorrelate spectra [not used]"""
207 crosssplist = []
208 for sp in spectralist:
209 c = crosscorrspec(sp, template)
210 c.crosscorr(fluxchoice)
211 crosssplist.append(c)
212 return crosssplist
213