Create s1d outputs

README.md

@@ -7,3 +7,8 @@ Python packages:
  - NumPy v1.12
  - SciPy v0.19
  - watchdog v0.8.2
+
+Output:
+ - intermediate: *_str.fits, 2D images with orders straightened and bad pixel removed
+ - final output: *_e2ds.fits. The spectra are flat-fielded, extracted, calibrated, and the cosmic rays are removed. The echelle orders are not merged. The FITS files contain three images 50 orders x 2048 pixels: the primary image contains the fluxes, the second contains the wavelength calibration and the third contains the signal-to-noise ratio.
+ - final output (optional): *_s1d.fits. The spectra are flat-fielded, extracted, calibrated, and the cosmic rays are removed. The echelle orders are merged. These are monodimensional spectra with constant step in wavelength. 

drslib/config.py

@@ -188,6 +188,8 @@ CONFIG['A_POS'] = 8. # guess pixel position for A
 CONFIG['B_POS'] = 23. # guess pixel position for B
 CONFIG['C_POS'] = 15. # guess pixel position for C
 CONFIG['HWTM'] = 6. # hwtm limit for contamination between A and B
+CONFIG['S1D'] = True # flag for creation of s1d output (True/False)
+CONFIG['S1D_STEP'] = 0.002 # constant step to use in s1d spectra
 
 #
 # Wavelength calibration
@@ -205,22 +207,35 @@ CONFIG['L0_GUESS'] = {32:2398.32, 33:2325.63, 34:2257.24, 35:2192.68, 36:2131.75
 
 CONFIG['DO_CALIB'] = {'dark':True, 'flat':True, 'une':True, 'fp':True, 'only_calib':False}
 #CONFIG['DO_CALIB'] = {'dark':False, 'flat':False, 'une':False, 'fp':False, 'only_calib':False}
+
+CONFIG['USE_FLAT'] = {'global':False, 'order':True, 'nor':False} # global: flat normalized by global mean value - order: flat normalized order by order by mean value - nor: flat normalized row by row by continuum fitting, removing the blaze function - No True value: no flat division
+
 CONFIG['FLAT_EXPT'] = 100.0
 CONFIG['WLCAL_EXPT'] = 100.0
 CONFIG['FP_EXPT'] = 60.0
 CONFIG['DARKLIST'] = {10:[],30:[],60:[],100:[]} # exptime dei dark
 CONFIG['NDARK'] = 4 # quanti tipi di dark (esempio: 10,30,60,100 sec = 4) meglio sovrastimare
+
+#
+# CCD info
+#
+
 CONFIG['RON'] = 5.0 # electrons
 CONFIG['GAIN'] = 2.2 # e/ADU
 CONFIG['XCCD'] = 2048
 CONFIG['YCCD'] = 2048
+CONFIG['N_ORD'] = 50 # number of echelle orders
+CONFIG['W_ORD'] = 41 # width of straighten order in pixel
+
+#
+# Quality checks
+#
+
 CONFIG['DARK_MEAN'] = {10:[5.0,5.0,10.0,8.0],30:[15.0,12.0,25.0,18.0],60:[25.0,20.0,40.0,32.0],100:[30.0,25.0,65.0,45.0]} # Dark quality check: mean signal below this number (4 quadrant for 4 exposure time)
-#CONFIG['FLATCHECK'] = [985,1030,575,590] # region x1:x2,y1:y2 for flats' quality check
 CONFIG['FLATSIGNAL'] = 300 # lowest flat signal expected to pass quality check
 CONFIG['SCIENCECHECK'] = [500,540,140,160] # region x1:x2,y1:y2 for images' quality check
 #CONFIG['NODSIGNAL'] = 2*CONFIG['DARK_MEAN']['quadrant'][0] # lowest image signal expected to pass quality check
 CONFIG['NODSIGNAL'] = 0.0 # lowest image signal expected to pass quality check
-CONFIG['USE_FLAT'] = {'global':False, 'order':True, 'nor':False} # global: flat normalized by global mean value - order: flat normalized order by order by mean value - nor: flat normalized row by row by continuum fitting, removing the blaze function - No True value: no flat division
 
 
 

drslib/flatframes.py

@@ -202,10 +202,10 @@ class GBFlats():
         xvalues = np.arange(CONFIG['XCCD'])
         np.seterr(divide='ignore')
 
-        for x in xrange(50):
+        for x in xrange(CONFIG['N_ORD']):
             # order limits
-            start = x*41
-            end = min(start+41,CONFIG['YCCD'])
+            start = x*CONFIG['W_ORD']
+            end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])
 
             # evaluate average background value for each order
             #background = []

drslib/nodding.py

@@ -295,14 +295,15 @@ class GBNodding():
 
         # prepare the structure for the calibrated results
         heacal = OrderedDict()
-        #stdSpectrum = np.zeros((50,CONFIG['YCCD']))
-        optSpectrum = np.zeros((50,CONFIG['YCCD']))
-        fsnr = np.zeros((50,CONFIG['YCCD']))
-        snr = np.zeros((50,CONFIG['YCCD']))
+        #stdSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
+        optSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
+        fsnr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
+        snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
 
-        for x in xrange(50):
-            start = x*41
-            end = min(start+41,CONFIG['YCCD'])
+        all_cosmics = 0
+        for x in xrange(CONFIG['N_ORD']):
+            start = x*CONFIG['W_ORD']
+            end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])
 
             # select only the rows wit the signal using the appropriate mask
 
@@ -323,7 +324,8 @@ class GBNodding():
             ordermasked = np.ma.MaskedArray(order,mask=omask)
             goodorder = np.ma.compress_rows(ordermasked)
 
-            optSpectrum[x],varOptFlux,profile,x1,x2 = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
+            optSpectrum[x],varOptFlux,profile,x1,x2,cosmics = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
+            all_cosmics = all_cosmics + cosmics
 
             #snr.append(round(max( (np.mean(optSpectrum[x][1000:1050])/np.std(optSpectrum[x][1000:1050])) ,0),2))
             #optSpectrum[x] = optSpectrum[x][::-1]
@@ -377,6 +379,11 @@ class GBNodding():
         optSpectrum = np.asarray(optSpectrum, dtype='float32')
 
         self.messages.append('The nodding %s was extracted.' % str(slit),)
+        if all_cosmics == 1050:
+            self.messages.append('%s cosmics were removed (maximum iteration reached).' % str(all_cosmics),)
+        else:
+            self.messages.append('%s cosmics were removed.' % str(all_cosmics),)
+
 
         if slit_pos == CONFIG['A_POS']:
             keyA = ''.join((CONFIG['SPEC_USED'][0],'1','A'))
@@ -384,10 +391,11 @@ class GBNodding():
             nomebase = os.path.splitext(aname)[0]
             #qui = aname.rindex('.')
             if 'grp' in fitsfile:
-                #calname = '_'.join((aname[0:qui],'Agrp.fits'))
-                calname = '_'.join((nomebase,'Agrp.fits'))
+                calname = '_'.join((nomebase,'Agrp_e2ds.fits'))
+                calname1d = '_'.join((nomebase,'Agrp_s1d.fits'))
             else:
-                calname = '_'.join((nomebase,'A.fits'))
+                calname = '_'.join((nomebase,'A_e2ds.fits'))
+                calname1d = '_'.join((nomebase,'A_s1d.fits'))
 
             #calname = str(os.path.basename(fitsfile)).replace('_AB','_A')
         elif slit_pos == CONFIG['B_POS']:
@@ -396,9 +404,11 @@ class GBNodding():
             nomebase = os.path.splitext(bname)[0]
             #qui = bname.rindex('.')
             if 'grp' in fitsfile:
-                calname = '_'.join((nomebase,'Bgrp.fits'))
+                calname = '_'.join((nomebase,'Bgrp_e2ds.fits'))
+                calname1d = '_'.join((nomebase,'Bgrp_s1d.fits'))
             else:
-                calname = '_'.join((nomebase,'B.fits'))
+                calname = '_'.join((nomebase,'B_e2ds.fits'))
+                calname1d = '_'.join((nomebase,'B_s1d.fits'))
             #calname = str(os.path.basename(fitsfile)).replace('_AB','_B')
         else: 
             print 'Wrong slit position!'
@@ -419,7 +429,7 @@ class GBNodding():
         self.messages.append('Nodding %s: SNR[K band, order=35, wl=2200 nm] = %s' % (str(slit),str(round(np.mean(snr[3]),2))),)
 
 
-        for o in xrange(50):
+        for o in xrange(CONFIG['N_ORD']):
             key_snr = ''.join((CONFIG['SNR'][0],str(o+32)))
             heaspe[key_snr] = (round(np.mean(snr[o]),2),CONFIG['SNR'][1])
 
@@ -431,8 +441,8 @@ class GBNodding():
         heaspe[CONFIG['BERV'][0]] = (barycorr,CONFIG['BERV'][1])
         heaspe[CONFIG['HJD'][0]] = (hjd,CONFIG['HJD'][1])
 
-        waves = np.zeros((50,CONFIG['YCCD']))
-        for o in xrange(50):
+        waves = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
+        for o in xrange(CONFIG['N_ORD']):
             waves[o] = varie.wcalib(heaspe,o)
 
         spefits = fits.PrimaryHDU(optSpectrum,header=heaspe)
@@ -445,6 +455,28 @@ class GBNodding():
         calname = os.path.join(CONFIG['RED_DIR'],calname)
         results.writeto(calname,clobber=True)
 
+        #t1 = time.time()
+
+        #print 's1d'
+        #print calname1d
+
+        if CONFIG['S1D']:
+            #s1d = varie.create_s1d(optSpectrum,snr,heaspe)
+            s1d = varie.create_s1d(optSpectrum,heaspe)
+            heaspe[CONFIG['KEYS']['FILENAME']] = calname1d
+            heaspe['CRPIX1'] = (1.,'Reference pixel')
+            heaspe['CRVAL1'] = (s1d.keys()[0],'Coordinates at reference pixel')
+            heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')
+
+            heaspe['CTYPE1'] = ('Nanometers','Units of coordinates')
+            heaspe['BUNIT'] = ('Relative Flux','Units of data values')
+
+            s1dfits = fits.PrimaryHDU(np.asarray(s1d.values()),header=heaspe)
+            calname1d = os.path.join(CONFIG['RED_DIR'],calname1d)
+            s1dfits.writeto(calname1d,clobber=True)
+
+        #t2 = time.time()
+        #print 's1d spectrum: %s s' %  str(t2-t1)
 
         if hea_ima[CONFIG['KEYS']['EXTMODE']] == CONFIG['EXTPAIR']:
             return calname, fsnr, straight
@@ -457,8 +489,12 @@ class GBNodding():
 
     def combine(self,acalib,bcalib,fsnr):
 
-        abnome = acalib.replace('_A.fits','_AB.fits')
-        abnome = abnome.replace('_Agrp.fits','_ABgrp.fits')
+        abnome = acalib.replace('_A_e2ds.fits','_AB_e2ds.fits')
+        abnome1d = acalib.replace('_A_e2ds.fits','_AB_s1d.fits')
+        abnome = abnome.replace('_Agrp_e2ds.fits','_ABgrp_e2ds.fits')
+        abnome1d = abnome1d.replace('_Agrp_e2ds.fits','_ABgrp_s1d.fits')
+        #print abnome
+        #print abnome1d
 
         acal = fits.open(acalib)
         afluxes = acal[0].data
@@ -471,11 +507,11 @@ class GBNodding():
         bfluxes = bcal[0].data
         bwaves = bcal[0].header
 
-        abcalib = np.zeros((50,CONFIG['YCCD']))
+        abcalib = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
 
-        snr = np.zeros((50,CONFIG['YCCD']))
+        snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
 
-        for o in xrange(50):
+        for o in xrange(CONFIG['N_ORD']):
             bshift = varie.rebin(awaves,bfluxes[o],bwaves,o)
             abcalib[o] = (afluxes[o]+bshift)/2.0
             #snr.append(max(np.mean(abcalib[o][1000:1050])/np.std(abcalib[o][1000:1050]),0))
@@ -502,7 +538,7 @@ class GBNodding():
         snr[snr==-np.inf] = 0
         snr = np.nan_to_num(snr)
 
-        for o in xrange(50):
+        for o in xrange(CONFIG['N_ORD']):
             key_snr = ''.join((CONFIG['SNR'][0],str(o+32)))
             abhea[key_snr] = (round(np.mean(snr[o]),2),CONFIG['SNR'][1])
 
@@ -531,8 +567,8 @@ class GBNodding():
             heaspe[CONFIG['AIRMASS'][0]] = (heaspe[CONFIG['KEYS']['AM']],CONFIG['AIRMASS'][1])
 
 
-        waves = np.zeros((50,CONFIG['YCCD']))
-        for o in xrange(50):
+        waves = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
+        for o in xrange(CONFIG['N_ORD']):
             waves[o] = varie.wcalib(heaspe,o)
 
         spefits = fits.PrimaryHDU(abcalib,header=heaspe)
@@ -545,11 +581,28 @@ class GBNodding():
 
         results.writeto(abnome,clobber=True)
 
+        if CONFIG['S1D']:
+
+            #s1d = varie.create_s1d(abcalib,snr,heaspe)
+            s1d = varie.create_s1d(abcalib,heaspe)
+            heaspe['CRPIX1'] = (1.,'Reference pixel')
+            heaspe['CRVAL1'] = (s1d.keys()[0],'Coordinates at reference pixel')
+            heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')
+
+            heaspe['CTYPE1'] = ('Nanometers','Units of coordinates')
+            heaspe['BUNIT'] = ('Relative Flux','Units of data values')
+
+            s1dfits = fits.PrimaryHDU(np.asarray(s1d.values()),header=heaspe)
+
+            s1dfits.writeto(abnome1d,clobber=True)
+
+
         self.messages.append('Nodding merged AB: SNR[Y band, order=73, wl=1050 nm] = %s' % (str(round(np.mean(snr[41]),2))),)
         self.messages.append('Nodding merged AB: SNR[J band, order=61, wl=1250 nm] = %s' % (str(round(np.mean(snr[29]),2))),)
         self.messages.append('Nodding merged AB: SNR[H band, order=46, wl=1650 nm] = %s' % (str(round(np.mean(snr[14]),2))),)
         self.messages.append('Nodding merged AB: SNR[K band, order=35, wl=2200 nm] = %s' % (str(round(np.mean(snr[3]),2))),)
 
+
         #calibrated = np.vstack((np.concatenate(np.flipud(awaves)),np.concatenate(np.flipud(abcalib))))
         #print calibrated
 
@@ -626,6 +679,7 @@ class GBNodding():
         nodAfits = fits.HDUList([hdu])
         nodAfits.writeto(Atmpnome,clobber=True)
 
+
         for n in self.group['noddings']:
             os.remove(n)
 

drslib/stare.py

@@ -435,15 +435,16 @@ class GBStare():
 
         # prepare the structure for the calibrated results
         heacal = OrderedDict() # header for the calibration table
-        #stdSpectrum = np.zeros((50,CONFIG['YCCD']))
-        optSpectrum = np.zeros((50,CONFIG['YCCD']))
-        fsnr = np.zeros((50,CONFIG['YCCD']))
-        snr = np.zeros((50,CONFIG['YCCD']))
+        #stdSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
+        optSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
+        fsnr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
+        snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
 
-        for x in xrange(50):
+        all_cosmics = 0
+        for x in xrange(CONFIG['N_ORD']):
 
-            start = x*41
-            end = min(start+41,CONFIG['YCCD'])
+            start = x*CONFIG['W_ORD']
+            end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])
 
             # select only the rows wit the signal using the appropriate mask
 
@@ -463,7 +464,8 @@ class GBStare():
             goodorder = np.ma.compress_rows(ordermasked)
 
             # call optimal extraction
-            optSpectrum[x],varOptFlux,profile,x1,x2 = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
+            optSpectrum[x],varOptFlux,profile,x1,x2,cosmics = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
+            all_cosmics = all_cosmics + cosmics
 
             olamp = mlamp.data[start:end]
             orderlamp = np.ma.MaskedArray(olamp,mask=omask)
@@ -512,9 +514,14 @@ class GBStare():
 
 
         self.messages.append('The spectrum %s was extracted.' % str(os.path.basename(straight)),)
+        if all_cosmics == 1050:
+            self.messages.append('%s cosmics were removed (maximum iteration reached).' % str(all_cosmics),)
+        else:
+            self.messages.append('%s cosmics were removed.' % str(all_cosmics),)
 
-        calname = os.path.join(CONFIG['RED_DIR'],str(os.path.basename(fitsfile)))
-
+        redname = os.path.join(CONFIG['RED_DIR'],str(os.path.basename(fitsfile)))
+        calname = redname.replace('.fits','_e2ds.fits')
+        calname1d = redname.replace('.fits','_s1d.fits')
 
         heaspe = fits.Header(imstr.header)
 
@@ -527,7 +534,7 @@ class GBStare():
         self.messages.append('Stare image: SNR[H band, order=46, wl=1650 nm] = %s' % (str(round(np.mean(snr[14]),2))),)
         self.messages.append('Stare image: SNR[K band, order=35, wl=2200 nm] = %s' % (str(round(np.mean(snr[3]),2))),)
 
-        for o in xrange(50):
+        for o in xrange(CONFIG['N_ORD']):
             key_snr = ''.join((CONFIG['SNR'][0],str(o+32)))
             heaspe[key_snr] = (round(np.mean(snr[o]),2),CONFIG['SNR'][1])
 
@@ -539,8 +546,8 @@ class GBStare():
         heaspe[CONFIG['BERV'][0]] = (barycorr,CONFIG['BERV'][1])
         heaspe[CONFIG['HJD'][0]] = (hjd,CONFIG['HJD'][1])
 
-        waves = np.zeros((50,CONFIG['YCCD']))
-        for o in xrange(50):
+        waves = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
+        for o in xrange(CONFIG['N_ORD']):
             waves[o] = varie.wcalib(heaspe,o)
 
         spefits = fits.PrimaryHDU(optSpectrum,header=heaspe)
@@ -554,6 +561,20 @@ class GBStare():
         calname = os.path.join(CONFIG['RED_DIR'],calname)
         results.writeto(calname,clobber=True)
 
+        if CONFIG['S1D']:
+            #s1d = varie.create_s1d(optSpectrum,snr,heaspe)
+            s1d = varie.create_s1d(optSpectrum,heaspe)
+            heaspe['CRPIX1'] = (1.,'Reference pixel')
+            heaspe['CRVAL1'] = (s1d.keys()[0],'Coordinates at reference pixel')
+            heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')
+
+            heaspe['CTYPE1'] = ('Nanometers','Units of coordinates')
+            heaspe['BUNIT'] = ('Relative Flux','Units of data values')
+
+            s1dfits = fits.PrimaryHDU(np.asarray(s1d.values()),header=heaspe)
+            calname1d = os.path.join(CONFIG['RED_DIR'],calname1d)
+            s1dfits.writeto(calname1d,clobber=True)
+
 
         if imstr.header[CONFIG['KEYS']['EXTMODE']] == CONFIG['EXTPAIR']:
             return calname, straight

drslib/varie.py

 """
-Last modified: 2017-03-07
+Last modified: 2017-06-07
 
 - badpix: bad pixels removals
 - stdcombine: weights for flat and dark combiner
@@ -7,8 +7,13 @@ Last modified: 2017-03-07
 - extract: extraction with pre-defined profiles
 - UNe_linelist: read the lists of UNe lines
 - UNe_calibrate: calibration with UNe lamps
-- rebin: rebin the B nodding on the A wavelengths prior to combine them
+- wcalib: apply wavelength calibration
+- rebin_linear: linearly rebin the B nodding on the A wavelengths prior to combine them
+- rebin2deg: parabolically rebin the B nodding on the A wavelengths prior to combine them
+- rebin: call either rebin_linear or rebin2deg (to change manually)
+- check_keyraw/check_keywords: check for keyword existence 
 - berv: computation of barycentric velocity correction (to be updated)
+- create_s1d: create s1d output
 """
 
 from drslib.config import CONFIG
@@ -89,10 +94,10 @@ def buildMaskC(fdata):
 
     maskC = np.ones((CONFIG['YCCD'],CONFIG['XCCD']), dtype='int')
 
-    for x in xrange(50):
+    for x in xrange(CONFIG['N_ORD']):
         # order limits
-        start = x*41
-        end = min(start+41,CONFIG['YCCD'])
+        start = x*CONFIG['W_ORD']
+        end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])
 
         # evaluate average background value for each order
         #background = []
@@ -346,11 +351,17 @@ def optExtract(data,gain,ron,slit_pos,ordine):
             #model[rworst,cworst] = 0
             variance = var(model)
 
+            stop += 1
+            if stop > 20:
+                cosmic = False
+
         else:
             cosmic = False               
-        stop += 1
-        if stop > 20:
-            cosmic = False
+        #stop += 1
+        #if stop > 20:
+        #    cosmic = False
+
+    #print stop
 
     # clean the extracted spectra from NaN and infinite values
 
@@ -374,7 +385,7 @@ def optExtract(data,gain,ron,slit_pos,ordine):
     #warnings.simplefilter('default',RuntimeWarning)
     warnings.simplefilter("default", optimize.OptimizeWarning)
 
-    return OptFlux, varOptFlux, profile, x1, x2
+    return OptFlux, varOptFlux, profile, x1, x2, stop
 
 
 #--------------------- Extraction with a pre-determined profile -------------------
@@ -402,7 +413,7 @@ def extract(data,optflux,x1,x2,profile,gain,ron):
 
 
     #plt.plot(OptFlux)
-    #plt.axis([0,2048,-20,2000])
+    #plt.axis([0,CONFIG['YCCD'],-20,2000])
     #plt.show()
 
     return OptFlux
@@ -727,19 +738,14 @@ def UNe_calibrate(lamp,order,select_lines,all_lines):
     return calib_failed, coeffs, messages
 
 
-def rebin(heawave,flux_old,heawave_old,o):
-    """
-    Rebin spectrum from wave_old to wave, interpolating linearly
-    """
+#--------------------- Apply wavelength calibration -------------------
 
+def wcalib(heawave,o):
     def lambdafit(x,lambda0,xc0):
-        return lambda0 + k1*(len(pixrange)-x-xc0) + k2*(len(pixrange)-x-xc0)**2 + k3*(len(pixrange)-x-xc0)**3
+        return lambda0 + k1*(x-xc0) + k2*(x-xc0)**2 + k3*(x-xc0)**3
 
     keyk1 = ''.join((CONFIG['WLCOEFFS']['k1'][0],str(o+32)))
-    try:
-        k1 = float(heawave[keyk1])
-    except:
-        print heawave[keyk1]
+    k1 = float(heawave[keyk1])
     keyk2 = ''.join((CONFIG['WLCOEFFS']['k2'][0],str(o+32)))
     k2 = float(heawave[keyk2])
     keyk3 = ''.join((CONFIG['WLCOEFFS']['k3'][0],str(o+32)))
@@ -749,26 +755,26 @@ def rebin(heawave,flux_old,heawave_old,o):
     keyxc = ''.join((CONFIG['WLCOEFFS']['xc'][0],str(o+32)))
     xc = float(heawave[keyxc])
 
-
-    pixrange = np.arange(len(flux_old))+1
+    #pixrange = -np.arange(CONFIG['YCCD'])+CONFIG['YCCD']
+    pixrange = np.arange(CONFIG['YCCD'])+1
     wave = lambdafit(pixrange,l0,xc)
 
-    keyk1 = ''.join((CONFIG['WLCOEFFS']['k1'][0],str(o+32)))
-    k1 = float(heawave_old[keyk1])
-    keyk2 = ''.join((CONFIG['WLCOEFFS']['k2'][0],str(o+32)))
-    k2 = float(heawave_old[keyk2])
-    keyk3 = ''.join((CONFIG['WLCOEFFS']['k3'][0],str(o+32)))
-    k3 = float(heawave_old[keyk3])
-    keyl0 = ''.join((CONFIG['WLCOEFFS']['l0'][0],str(o+32)))
-    l0 = float(heawave_old[keyl0])
-    keyxc = ''.join((CONFIG['WLCOEFFS']['xc'][0],str(o+32)))
-    xc = float(heawave_old[keyxc])
+    return wave
 
-    wave_old = lambdafit(pixrange,l0,xc)
 
-    wave = wave[::1]
-    wave_old = wave_old[::1]
-    flux_old = flux_old[::1]
+#--------------------- Linear interpolation and rebinning -------------------
+
+def rebin_linear(heawave,flux_old,heawave_old,o):
+    """
+    Rebin spectrum from wave_old to wave, interpolating linearly
+    """
+
+    wave = wcalib(heawave,o)[::-1]
+    wave_old = wcalib(heawave_old,o)[::-1]
+    flux_old = flux_old[::-1]
+
+    #print wave
+    #print wave_old
 
     #flux_new = np.interp(wave,wave_old,flux_old)
 
@@ -794,9 +800,99 @@ def rebin(heawave,flux_old,heawave_old,o):
     #plt.plot(flux_new)
     #plt.show()
 
-    return np.asarray(flux_new)[::1]
+    return np.asarray(flux_new)[::-1]
+    #return np.asarray(flux_new)
+
+
+#--------------------- Parabolic interpolation and rebinning -------------------
+
+def rebin2deg(heawave,flux_old,heawave_old,o):
+    """
+    Rebin spectrum from wave_old to wave, interpolating 2 degree
+    """
+
+    wave = wcalib(heawave,o)[::-1]
+    wave_old = wcalib(heawave_old,o)[::-1]
+    flux_old = flux_old[::-1]
+
+    #print wave
+    #print wave_old
+
+    #flux_new = np.interp(wave,wave_old,flux_old)
+
+    flux_new = []
+
+    for w in wave:
+        iw = min(np.searchsorted(wave_old,w),len(wave_old)-1)
+        if wave_old[iw] == w:
+            flux_new.append(flux_old[iw])
+        else:
+            try:
+                y0 = flux_old[iw-1]
+                y1 = flux_old[iw]
+                y2 = flux_old[iw+1]
+
+                x0 = wave_old[iw-1]
+                x1 = wave_old[iw]
+                x2 = wave_old[iw+1]
+
+                flux = y0*(w-x1)*(w-x2)/((x0-x1)*(x0-x2)) + \
+                       y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
+                       y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
+
+                flux_new.append(flux)
+
+            except:
+                try:
+                    y0 = flux_old[iw-2]
+                    y1 = flux_old[iw-1]
+                    y2 = flux_old[iw]
+
+                    x0 = wave_old[iw-2]
+                    x1 = wave_old[iw-1]
+                    x2 = wave_old[iw]
+
+                    flux = y0*(w-x1)*(w-x2)/((x0-x1)*(x0-x2)) + \
+                           y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
+                           y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
+
+                    flux_new.append(flux)
+
+                except:
+                    try:
+                        y0 = flux_old[iw]
+                        y1 = flux_old[iw+1]
+                        y2 = flux_old[iw+2]
 
+                        x0 = wave_old[iw]
+                        x1 = wave_old[iw+1]
+                        x2 = wave_old[iw+2]
 
+                        flux = y0*(w-x1)*(w-x2)/((x0-x1)*(x0-x2)) + \
+                               y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
+                               y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
+
+                        flux_new.append(flux)
+
+                    except:
+                        flux_new.append(flux_old[iw])
+    #plt.plot(flux_old)
+    #plt.plot(flux_new)
+    #plt.show()
+
+    return np.asarray(flux_new)[::-1]
+    #return np.asarray(flux_new)
+
+#--------- Switch between linear and parabolic interpolation and rebinning --------
+
+def rebin(heawave,flux_old,heawave_old,o):
+    """
+    Decide which rebinning to use: linear or parabolic
+    """
+    #return rebin_linear(heawave,flux_old,heawave_old,o)
+    return rebin2deg(heawave,flux_old,heawave_old,o)
+
+#--------------------- Check for keywords existence (interactive) -------------------
 
 def check_keyraw(header,filename):
     keys = [CONFIG['KEYS']['ID'] , CONFIG['KEYS']['PID'] , CONFIG['KEYS']['NODSTARE'] , CONFIG['KEYS']['EXTMODE'] , CONFIG['KEYS']['GROUPI'] , CONFIG['KEYS']['GROUPN']]
@@ -823,6 +919,7 @@ def check_keyraw(header,filename):
 
     return header
 
+#--------------------- Check for keywords existence (not interactive) -------------------
 
 def check_keywords(header,filename):
     for key in CONFIG['KEYS']:
@@ -839,27 +936,7 @@ def check_keywords(header,filename