Modifications to mktelluric

lcogtgemini/__init__.py

@@ -351,11 +351,20 @@ def get_chipedges(data):
         try:
             w = np.where(data == 0.0)[0]
 
+
             # Note we also grow the chip gap by 8 pixels on each side
             # Chip 1
             chip_edges = []
             left_chipedge = 10
             morechips = True
+
+            # This is necessary after speccombine.
+            # If we have only the red side of a spectrum,
+            # the first ~1000 elements of the combined spectrum are 0
+            if len(w) > 1000 and w[1000] == 1000:
+                    blue_cut = np.min(np.where(data > 0))
+                    left_chipedge=blue_cut + 10
+
             while morechips:
                 try:
                     right_chipedge = np.min(w[w > (left_chipedge + 25)]) - 10
@@ -367,6 +376,7 @@ def get_chipedges(data):
                 left_chipedge = np.max(w[w < (right_chipedge + 200)]) + 10
         except:
             chip_edges = []
+
         return chip_edges
 
 
@@ -581,13 +591,21 @@ def mktelluric(filename):
     spec[chip_gaps] = interpolate.splev(waves[chip_gaps], intpr)
 
     not_telluric = telluric_mask(waves)
+
     # Smooth the spectrum so that the spline doesn't go as crazy
     # Use the Savitzky-Golay filter to presevere the edges of the
     # absorption features (both atomospheric and intrinsic to the star)
     sgspec = signal.savgol_filter(spec, 31, 3)
     # Get the number of data points to set the smoothing criteria for the 
     # spline
     m = not_telluric.sum()
+
+    # If we have only the red side of a spectrum,
+    # the first ~1000 elements of the combined spectrum and noise are 0,
+    # so replace zeros of noise with 1e10 so that weights, w, for the following interpolation are 0
+    if len(noise) > 100 and np.sum(noise[:100]==0.0)==100:
+    	noise[noise<1e-10]=1e10
+
     intpr = interpolate.splrep(waves[not_telluric], sgspec[not_telluric],
                                w=1 / noise[not_telluric], k=2, s=20 * m)
 
@@ -597,15 +615,18 @@ def mktelluric(filename):
     # Extrapolate the ends linearly
     # Blue side
     w = np.logical_and(waves > 3420, waves < 3600)
-    bluefit = np.poly1d(np.polyfit(waves[w], spec[w], 1))
-    bluewaves = waves < 3420
-    smoothedspec[bluewaves] = bluefit(waves[bluewaves])
+    if len(waves[w]) > 0:
+    	bluefit = np.poly1d(np.polyfit(waves[w], spec[w], 1))
+    	bluewaves = waves < 3420
+    	smoothedspec[bluewaves] = bluefit(waves[bluewaves])
 
     # Red side
     w = np.logical_and(waves > 8410, waves < 8800)
-    redfit = np.poly1d(np.polyfit(waves[w], spec[w], 1))
-    redwaves = waves > 8800
-    smoothedspec[redwaves] = redfit(waves[redwaves])
+    if len(waves[w]) > 0:
+    	redfit = np.poly1d(np.polyfit(waves[w], spec[w], 1))
+    	redwaves = waves > 8800
+    	smoothedspec[redwaves] = redfit(waves[redwaves])
+
     smoothedspec[not_telluric] = spec[not_telluric]
     # Divide the original and the telluric corrected spectra to
     # get the correction factor
@@ -636,7 +657,12 @@ def telluric(filename, outfile):
     # to find wavelength shift of standard star.
     w = np.logical_and(waves > 7550., waves < 8410.)
     tw = np.logical_and(telwave > 7550., telwave < 8410.)
-    p = fitxcor(waves[w], spec[w], telwave[tw], telspec[tw])
+    # If the main telluric bands overlap the spectrum, refine the wavelength solution of
+    # the correction using cross correlation.
+    if w.sum() > 5 and tw.sum() > 5:
+	    p = fitxcor(waves[w], spec[w], telwave[tw], telspec[tw])
+    else:
+    	p = [1.0, 0.0]
     # shift and stretch standard star spectrum to match science
     # spectrum.
     telcorr = np.interp(waves, p[0] * telwave + p[1], telspec)
@@ -894,6 +920,7 @@ def makemasterflat(flatfiles, rawpath, plot=True):
 
         y = data / np.median(data)
 
+
         fitme_x = x[chip_edges[0][0]:chip_edges[0][1]]
         fitme_x = np.append(fitme_x, x[chip_edges[1][0]:chip_edges[1][1]])
         fitme_x = np.append(fitme_x, x[chip_edges[2][0]:chip_edges[2][1]])