Package pysls :: Module sp
[frames] | no frames]

Source Code for Module pysls.sp

  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   
16 -class spectra:
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 #an array of wavelengths (in angstroms) 21 self.fluxes = fluxes #flux (with continuum) 22 self.fluxes_cs = fluxes_cs # flux (without continuum) 23 self.fluxes_csc = np.array([]) #flux (without continuum, calculated) 24 self.noise = noise 25 self.zsdss = zsdss # redshift calculated by SDSS 26 self.secz = 0 27 self.mask = np.array([]) 28 self.plotcolour = "black" 29 self.coeff0 = coeff0 # Center wavelength (log10) of first pixel 30 self.coeff1 = coeff1 # Log10 dispersion per pixel 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
36 - def initmask(self):
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
41 - def copy(self):
42 """Returns a deep copy of itself""" 43 return pythoncopy.deepcopy(self)
44
45 - def pixtoangst(self, index):
46 """converts a pixel to a wavelenght in angstroems, using a formula from SDSS site""" 47 return 10**(self.coeff0 + self.coeff1 * index)
48
49 - def angsttopix(self, wavelength):
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
53 - def pixtoz(self, index):
54 """converts a pixel to a redshift""" 55 return 10**(self.coeff1*index)-1
56
57 - def ztopix(self, z):
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) #individual pixel number array 64 self.angst =np.array( map(self.pixtoangst,i) )
65
66 - def setz(self, z):
67 """sets the wavelengths corrsponding to the redshift""" 68 self.angst =(1 + z) / (1 + self.zsdss) * self.angst
69
70 - def calc_csc(self):
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
100 - def mask_angst(self, tupleslist, z = 0):
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
110 - def mask_height(self, height, fluxchoice):
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
156 -class crosscorrspec:
157 """This is a crosscorrelation spectra : it contains the origin spectra, a template which we use for the crosscorrelation and the result"""
158 - def __init__(self, spectra = spectra(), template = tem.template(), corr = np.array([])):
159 self.spectra = spectra 160 self.template = template 161 self.corr = corr
162
163 - def crosscorr(self, fluxchoice):
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
190 -def spfactory(directory):
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