;+
; NAME:
;   qgemspec.pro
;
; PURPOSE:
;   Fully-automated quick reduction of GMOS spectra
;
; CALLING SEQUENCE:
;   qgemspec, files
;
; INPUTS:
;   files       - Text file or array containing list of all arcs, flats, and science exposures.
;
; OPTIONAL INPUTS:
;   nocrreject  - keyword to skip cosmic ray rejection
;   ysearch     - array of [min, max] y-coordinates where the afterglow may be.
;
; COMMENTS:
;    Automatically identifies the flats and arcs, matches files with like configurations,
;    finds the source, extracts the spectrum, gets a wavelength solution.
;    Writes out several files:
;        [obj]_[config]#.dat : ascii files (four columns:  wave, flux, skyflux, pixel)
;        specplots.ps        : plots
;        specplotsnorm.ps    : plots normalized to continuum (approximate)
;        detectedlines.dat   : list of possible spectral lines
;        arcplots_c?.fits    : arc solution plots (quick sol at left, fine sol at right)
;        [filename]z.fits    : central part of science exposure with cosmic rays removed
;        outflat_c?.fits     : normalized flatfield in each configuration
; 
;    Has been tested only with the R400, B400, and R831 gratings.  Failures 
;    with the wavelength solution for different gratings are plausible.
;    Gemini-South only (will not properly process GMOS-N).
;    Requires the GSFC IDL libraries.
;
; HISTORY:
;    v0.0 2009-11-06 (DAP) - beta release to group
;
;
; Planned features:
;    Options to tweak CR rejection
;    Checks on tracing (to try to confidently contract the aperture)
;    Sky lines (variance spectrum) on plots
;    Thoroughput normalization
;    Improved normalization and absorption line detection
;    Absorption line / redshift identification
;
;-
;------------------------------------------------------------------------------


function countlines, infile
   openr, 1, infile
   iline=''
   nlines = 0l
   while not eof(1) do begin
     readf, 1, iline
     nlines = nlines + 1
   endwhile
   close, 1
   return, nlines
end

function grabcolumn, infile, colnum, extstr=extstr, str=str, skiplines=skiplines
   nlines = countlines(infile)
   str = keyword_set(str)
   if str then arr = strarr(nlines) else arr = dblarr(nlines)
   if n_elements(skiplines) eq 0 then skiplines = 0
   iline = ''
   ndat = 0
   openr, 1, infile
   for i = 0, nlines-1 do begin
     readf, 1, iline
     if i lt skiplines then continue
     if strmid(iline, 0, 1) eq '#' then continue
     if n_elements(extstr) eq 0 then $
        idata = strtrim(strsplit(iline, /extract),2) $
     else $
        idata = strtrim(strsplit(iline, extsr, /extract),2)
     if n_elements(idata) le colnum then continue
     if str then arr[ndat] = idata[colnum] else arr[ndat] = double(idata[colnum])
     ndat = ndat + 1
   endfor  
   close, 1
   return, arr[0:ndat-1]
end

function clip, strin, maxlen
  str = string(strin)
  if n_params() eq 1 then return, strtrim(str,2)
  str = strmid(strtrim(str,2),0,maxlen)
  curlen = strlen(str)
  if curlen lt maxlen then $
    for i = 0, (maxlen-curlen-1) do $
      str = str + ' '
  return, str
end

function unique, arr
  uniquearr = arr
  n = n_elements(arr)
  nunique = 0
  for i = 0, n-1 do begin
    known = 0
    for j = 0, nunique - 1 do begin
      if arr[i] eq uniquearr[j] then known = 1
    endfor
    if known eq 0 then begin
      uniquearr[nunique] = arr[i]
      nunique = nunique + 1
    endif
  endfor
  uniquearr = uniquearr[0:nunique-1]
  return, uniquearr
end


function loadandmerge, fitsfile, h
   ; It's easier to do operations once instead of three times, so stick file mosaics together into one image.

   if n_elements(fitsfile) eq 0 then begin
      print, 'ERROR: No file specified!'
      return, 0
   endif
   dummy = mrdfits(fitsfile, 0, h, /silent) ; just to get the header
   arr1 = mrdfits(fitsfile, 1, h1, /silent, /unsigned)
   arr2 = mrdfits(fitsfile, 2, h2, /silent, /unsigned)
   arr3 = mrdfits(fitsfile, 3, h3, /silent, /unsigned)
   
   ymax = 2304
   a1xmax = 1024-1     ; The array actually goes up to 1055.  But only 1025 is visible in the file...
   a2xmin = 1024
   a2xmax = 1024*2-1
   a3xmin = 1024*2 

   arr = fltarr(1024*3-33, ymax)
   arr[0:a1xmax,*] = float(arr1[0:a1xmax,*])
   arr[a2xmin:a2xmax,*] = float(arr2[0:a1xmax,*])    
   arr[a3xmin:*,*] = float(arr3[33:a1xmax,*])

   return, arr
end

function quickzap, imdata
      nx = (size(imdata)) [1]
      ny = (size(imdata)) [2]
      gain = 1.0    ; 2.2 if not pre-corrected
      rdnoise = 4

      skyimage = fltarr(nx,ny)
      for x = 0, nx-1 do begin        ; could break this down in y too
        skyimage[x,*] = median(imdata[x,*], /even) 
      endfor
      skysubimage = imdata - skyimage
 
      sourceimage = fltarr(nx, ny)
      for x = 0, nx-1, 3 do begin
         for y = 0, ny-1 do begin
           val = median(skysubimage[x-7>1:x+7<nx-2,y], /even) 
           sourceimage[x,y] = val
           if x ne nx-1 then sourceimage[x+1,y] = val
           if x ne 0 then sourceimage[x-1,y] = val
         endfor
      endfor
      
      skysourceimage = skyimage + sourceimage	
      excessimage = imdata - skysourceimage
      sigmaimage = excessimage / sqrt(skysourceimage*gain + rdnoise^2)

      patchimage = fltarr(nx, ny)

      nbad = 0
      d0 = nx
      for x = 1, nx-1-1 do begin
         for y = 1, ny-1-1 do begin
            if sigmaimage[x,y] lt 4.0 then continue

            n1neighborblockpos = [(y-1)*d0+(x-1),(y-1)*d0+x,(y-1)*d0+(x+1), $
                                  (y  )*d0+(x-1),           (y  )*d0+(x+1), $
                                  (y+1)*d0+(x-1),(y+1)*d0+x,(y+1)*d0+(x+1)]

            neighborratio = excessimage[n1neighborblockpos] / excessimage[x,y]

            dum = where(neighborratio gt 4, cthineighbor)
            dum = where(neighborratio lt 0.25, ctloneighbor)
            dum = where(neighborratio gt 0.7 and neighborratio lt 1.3, ctsimilarneighbor)

            if ((cthineighbor ge 1 or ctloneighbor ge 1) and ctsimilarneighbor lt 2) then begin
               patchimage[x,y] = 1
               patchimage[x+1,y] = 1
               patchimage[x-1,y] = 1
               patchimage[x,y+1] = 1
               patchimage[x,y-1] = 1
               nbad = nbad + 1
            endif
         endfor
      endfor

      patches = where(patchimage gt 0)
      imdata[patches] = skysourceimage[patches]

      ;writefits, 'testpatch.fits', patchimage
      ;writefits, 'testcrremoved.fits', imdata

      return, imdata

end

function mode, l
   if n_elements(l) eq 1 then return, l
   if n_elements(l) le 3 then return, median(l, /even)
   s = l[sort(l)]
   d = s[1:n_elements(s)-1] - s[0:n_elements(s)-2]
   nd = n_elements(d)
   g = nd / 16 > 2  ; sensitive to clusters up to a little less than 1/16 of the data set
   if nd lt 12 then g = 1
   minmean = total(d)
   imean = nd / 2
   for i = 0, nd-1 do begin
       r = [i-g>0,i+g<nd-1]
       m = mean(d[r[0]:r[1]])
       if m lt minmean then begin 
          minmean = m
          imean = i
       endif
   endfor
   md = s[imean] ;+ s[imean+1])/2
   return, md
end

function solvearc, linex, brightlinesx, reflines, brightreflines, debug=debug, wav0range=wav0range, disprange=disprange, maxsep=maxsep

         ; linex = detected line list of pixel centers
         ; brightlinesx = brightest ~10 lines, sorted by brightness
         ; reflines = good reference lines, including faint ones
         ; brightreflines = bright reference lines for coarse matching

         ; Currently only coarse matching is done.  Absolutely no improvement is attempted 
         ; after the probable solution approximation is reached.

         ; This routine works by looking for trios of strong lines in the arc spectrum and in the reference list.
         ; The ratio of the spacing between (center line to furthest line)/(center line to nearer line) is calculated,
         ; which is a scale-independent quanity.  Every possible trio is matched to every other possible trio and this
         ; ratio is tabulated, then a simple routine looks for a cluster of points around a similar value.  
         ; Using lines in matched trios the dispersion and wavelength at pix=0 are then estimated.

         if n_elements(wav0range) eq 0 then wav0range = [5000,12000]
         if n_elements(disprange) eq 0 then disprange = [-2.0, -0.2]
         if n_elements(maxsep) eq 0 then maxsep = 2000.

         if disprange[1] lt disprange[0] then disprange = disprange[[1,0]]

         ndetlines = n_elements(brightlinesx)
         xdiffrat = fltarr(ndetlines, ndetlines, ndetlines)
         for i = 0, ndetlines-1 do begin
            for j = 0, ndetlines-1 do begin
               for k = 0, ndetlines-1 do begin
                  if i eq j or j eq k or i eq k then continue
                  rat = (1.0*brightlinesx[j]-brightlinesx[i])/(1.0*brightlinesx[k]-brightlinesx[i])
                  if rat lt -1 then xdiffrat[i,j,k] = rat   ; ensures that i is in the middle and bigdiff/smalldiff
               endfor
            endfor
         endfor

         nreflines = n_elements(brightreflines)
         wavdiffrat = fltarr(nreflines, nreflines, nreflines)
         for l = 0, nreflines-1 do begin
            for m = 0, nreflines-1 do begin
               for n = 0, nreflines-1 do begin
                  if l eq m or m eq n or l eq n then continue
                  rat = (brightreflines[m]-brightreflines[l])/(brightreflines[n]-brightreflines[l])
                  if rat lt -1 and abs(brightreflines[m]-brightreflines[n]) lt maxsep then $
                     wavdiffrat[l,m,n] = rat  ; ensures that l is in the middle and bigdiff/smalldiff
                                              ; and lines are within maxsep Angstroms
               endfor
            endfor
         endfor

         maxmatch = 1024
         candidatematchi = intarr(maxmatch)
         candidatematchj = intarr(maxmatch)
         candidatematchk = intarr(maxmatch)
         candidatematchl = intarr(maxmatch)
         candidatematchm = intarr(maxmatch)
         candidatematchn = intarr(maxmatch)
         candidatematchdispersion = fltarr(maxmatch)
         candidatematchwav0 = fltarr(maxmatch)      

         nmatch = 0
         for i = 0, ndetlines-1 do begin
            for j = 0, ndetlines-1 do begin
               for k = 0, ndetlines-1 do begin         ; that's right, 6 nested loops in IDL
                  if xdiffrat[i,j,k] eq 0 then continue
                  for l = 0, nreflines-1 do begin
                     for m = 0, nreflines-1 do begin
                        for n = 0, nreflines-1 do begin
                           if wavdiffrat[l,m,n] eq 0 then continue
                           if abs(((xdiffrat[i,j,k]-wavdiffrat[l,m,n])/ $
                                   (xdiffrat[i,j,k]+wavdiffrat[l,m,n]))) lt 0.02/2 then begin
                                ; possible match
                              dispersion = (brightreflines[m]-brightreflines[l]) / $
                                           (brightlinesx[j]-brightlinesx[i])
                              wav0 = brightreflines[l] - brightlinesx[i]*dispersion
                              if dispersion gt disprange[0] and dispersion lt disprange[1] and $
                                    wav0 gt wav0range[0] and wav0 lt wav0range[1] then begin
                                 candidatematchdispersion[nmatch] = dispersion
                                 candidatematchwav0[nmatch] = wav0
                                 candidatematchi[nmatch] = i
                                 candidatematchj[nmatch] = j
                                 candidatematchk[nmatch] = k
                                 candidatematchl[nmatch] = l
                                 candidatematchm[nmatch] = m
                                 candidatematchn[nmatch] = n

                                 nmatch = nmatch + 1
                                 if keyword_set(debug) then print, clip(i,3), clip(j,3), clip(k,3), clip(l,3), clip(m,3), clip(n,3), dispersion, wav0, xdiffrat[i,j,k], wavdiffrat[l,m,n]
                                 if nmatch ge maxmatch then begin 
                                    print, 'Maximum possible matches reached!! Cycling back to zero (may result in bad solution)'
                                    nmatch = 0
                                 endif
                              endif
                           endif
                        endfor
                     endfor
                  endfor
               endfor
            endfor
         endfor

         candidatematchdispersion = candidatematchdispersion[0:nmatch-1]
         candidatematchwav0 = candidatematchwav0[0:nmatch-1]


      dispersion = mode(candidatematchdispersion)
      wav0 = mode(candidatematchwav0)               ; even better would be to do this in 2D


      goodmatch = where(abs(candidatematchdispersion - dispersion)/abs(dispersion) lt 0.1 and $
                        abs(candidatematchwav0 - wav0)/abs(wav0) lt 0.05)

      ; collect all the indices together.  There are probably lots of repeats.
      xmatchindex = [candidatematchi[goodmatch], candidatematchj[goodmatch], candidatematchk[goodmatch]]
      wavmatchindex = [candidatematchl[goodmatch], candidatematchm[goodmatch], candidatematchn[goodmatch]]

      ; Find the "real" matches by majority vote (which wav-index corresponds most often to each x-index)
      xmatches = unique(xmatchindex) ; list of x indices that actually match lines, no repeats
      wavmatches = intarr(n_elements(xmatches)) - 1 ; list of wav indices corresponding to above
      for iii = 0, n_elements(xmatches)-1 do begin
         wm = where(xmatchindex eq xmatches[iii], ct)
         if ct eq 0 then continue
         wmw = wavmatchindex[wm]
         wmwhist = histogram(wmw,min=0)
         if max(wmwhist) gt 1 then $
            wavmatches[iii] = (where(wmwhist eq max(wmwhist))) [0]
      endfor

      bestmatches = where(wavmatches ge 0, ct)
      if ct gt 3 then begin
         bestxmatches = xmatches[bestmatches]
         bestwavmatches = wavmatches[bestmatches]

         ;print, 'Matched ', clip(n_elements(bestmatches)), ' lines'

         ;print, wav0, dispersion
         fit = poly_fit(brightlinesx[bestxmatches], brightreflines[bestwavmatches], 1, yfit=yfit)

         ; remove outliers one by one
         noutlier = 0
         absdev = abs(yfit - brightreflines[bestwavmatches])

         ;print, absdev
         ;print, max(absdev)

            while max(absdev) gt 3 do begin
               notworstoutlier = where(absdev lt max(absdev), ctok)
               bestxmatches = bestxmatches[notworstoutlier]
               bestwavmatches = bestwavmatches[notworstoutlier]
               fit = poly_fit(brightlinesx[bestxmatches], brightreflines[bestwavmatches], 1, yfit=yfit)
               absdev = abs(yfit - brightreflines[bestwavmatches])
               noutlier = noutlier + 1
            endwhile
            ;if noutlier gt 0 then print, 'Removed ', clip(noutlier), ' lingering outlier line matches.'


      endif else begin
         fit = [wav0, dispersion]
         bestxmatches = xmatches
         bestwavmatches = wavmatches
      endelse

      plot, brightlinesx[bestxmatches], brightreflines[bestwavmatches]-poly(brightlinesx[bestxmatches],fit), psym=1
      


      ; Now use ALL the detected lines, including faint ones, to get a refined solution.

      predwav = poly(linex, fit)
      matchx = intarr(n_elements(linex)) 
      matchwav = intarr(n_elements(linex)) 

      m = 0
      for i = 0, n_elements(linex)-1 do begin
         absdev = abs(reflines - predwav[i])
         closest = (where(absdev eq min(absdev))) [0]
         if absdev[closest] lt 5 then begin
            matchx[m] = i
            matchwav[m] = closest
            ;print, clip(linex[i]), ' @', clip(predwav[i]), ' matches ', clip(reflines[closest]), '    ', clip(absdev[closest])
            m = m + 1
         endif
      endfor

      if m lt 3 then begin
         print, 'WARNING: Insufficient line matches in full line list for precise solution.'
         return, [wav0, dispersion, 0]
      endif

      matchx = matchx[0:m-1]
      matchwav = matchwav[0:m-1]

      oldfit = fit
      fit = poly_fit(linex[matchx], reflines[matchwav], 1, yfit=yfit)
      
      noutlier = 0
      absdev = abs(yfit - reflines[matchwav])
      while max(absdev) gt 4 do begin
          notworstoutlier = where(absdev lt max(absdev), ctok)
          matchx = matchx[notworstoutlier]
          matchwav = matchwav[notworstoutlier]
          fit = poly_fit(linex[matchx], reflines[matchwav], 1, yfit=yfit)
          absdev = abs(yfit - reflines[matchwav])
          noutlier = noutlier + 1
      endwhile
      ;print, 'Flagged ', noutlier, ' additional outliers.'

      plot, linex[matchx], reflines[matchwav]-poly(linex[matchx],fit), psym=7


      if n_elements(matchx) lt 4 then begin
         return, [fit[0], fit[1], 0]
      endif

      ;print, 'First-order solution: ', fit[0], fit[1]


      ; Now fit a second-order function.

      fit2 = poly_fit(linex[matchx], reflines[matchwav], 2, yfit=yfit)

      x = 30*findgen(120)

      noutlier = 0
      absdev = abs(yfit - reflines[matchwav])
      while max(absdev) gt 0.4 do begin
          notworstoutlier = where(absdev lt max(absdev), ctok)
          matchx = matchx[notworstoutlier]
          matchwav = matchwav[notworstoutlier]
          fit2 = poly_fit(linex[matchx], reflines[matchwav], 2, yfit=yfit)
          absdev = abs(yfit - reflines[matchwav])
          noutlier = noutlier + 1
      endwhile
      ;print, 'Flagged ', noutlier, ' additional outliers.'
 
      oplot, linex[matchx], reflines[matchwav]-poly(linex[matchx],fit), psym=4
      oplot, x, poly(x,fit2)-poly(x,fit), linestyle=1


   return, [fit2[0], fit2[1], fit2[2]]
end



pro openplot, plotfilename, xs, ys
     set_plot, 'ps'
     device, filename=plotfilename
     !p.font = 0
     device, /helvetica, font_index = 17
     device, /symbol, font_index = 4
     device, /helvetica, /bold, font_index = 18
     device, /isolatin1, font_index = 19
     device, /inches, xsize=xs, ysize=ys
     device, /inches, yoffset=4.5+(5.5-ys)
     device, /color
     colors = transpose([[0,0,0],$
                         [220,220,220],$ 
                         [64, 64, 192],$ 
                         [64, 136,128],$ 
                         [64, 158, 64],$  
                         [112,112, 64],$ 
                         [192, 64, 64],$
                         [128, 64, 128]$
                    ])
      tvlct, colors
end





; --------------------------------------------------------------------------------------------------------------------



pro qgemspec, filename, nocrreject=nocrreject, ysearch=ysearch

; Quickly reduce Gemini spectroscopy data.

; Input - an unsorted list of raw files, including science frames, flats, and arcs, in any number of different grating configurations.

; Output - wavelength-calibrated spectrum of the primary target (ascii files, sav files, and plots.)

   if n_elements(nocrreject) eq 0 then crreject = 1 else crreject = 0

if n_elements(ysearch) eq 0 then ysearch=[1080,1300]
   ymax = 2304
   a1xmax = 1024-1 
   a2xmin = 1024
   a2xmax = 1024*2-1
   a3xmin = 1024*2
   a3xmax = 3038
   nx = 3039

if n_elements(filename) eq 0 then begin
  print, 'qgemspec: A quick-and-dirty Gemini spectrum reducer'
  print, 'Input:  List of Gemini files in any order, including science exposures, arcs, and flats, in one or several configurations.  Either an array or a text file.'
  print, 'Optional inputs:  /nocrreject to skip cosmic rejection'
  print, '                  ysearch=[min,max] to specify the region to search for the afterglow
  return
endif

if n_elements(filename) eq 1 then begin
   filearr = strsplit(filename,'.', /extract)
   fileroot = filearr[0]
   fileext = filearr[1]
   if fileext eq 'cat' or fileext eq 'lis' or fileext eq 'list' or fileext eq 'txt' then begin
      files = grabcolumn(filename,0,/str)
      qgemspec, files, nocrreject=nocrreject, ysearch=ysearch
      return
   endif else begin
      print, 'Error: file list not recognized or only one file found.'
   endelse
endif


nf = n_elements(filename)
obs = replicate({filename:'', type:'', grating:'', filter:'', cwave:'', config:0, inst:''}, nf)

for f = 0, nf-1 do begin
   l = mrdfits(filename[f],0,h,/silent, /unsigned)
   obs[f].filename = filename[f]
   obs[f].type = clip(sxpar(h,'OBSTYPE'))
   obs[f].grating = clip(sxpar(h,'GRATING'))
   obs[f].filter = clip(sxpar(h,'FILTER1'))
   obs[f].cwave = clip(sxpar(h,'GRWLEN'))
   obs[f].inst = clip(sxpar(h,'INSTRUME'))
endfor
filtersdiffer = n_elements(unique(obs.filter)) gt 1
if filtersdiffer then filtstr = 'filter' else filtstr = ''
print, 'Filename                 type    grating       cwave   ' + filtstr
for f = 0, nf-1 do begin
   enddir = strpos(filename[f],'/',/reverse_search)+1
   if enddir ne -1 then shortfilename = strmid(filename[f],enddir,strlen(filename[f])) else shortfilename = filename[f]
   if filtersdiffer then filtstr = clip(obs[f].filter)
   print, clip(shortfilename,25), clip(obs[f].type,8) , clip(obs[f].grating,14), clip(fix(obs[f].cwave),8), filtstr
endfor

; Remove imaging observations like acquisition spectra

wspec = where(obs.grating ne 'MIRROR' and obs.grating ne 'Mirror', ct)

if ct eq 0 then begin
   print, 'No spectra found!'
   print, 'All files appear to be in imaging mode.  Please download the spectroscopic files'
endif

obs = obs[wspec]
nf = n_elements(obs)

config = strarr(nf)
for f = 0, nf-1 do begin
   config[f] = obs[f].grating + '/' + obs[f].cwave
endfor
configlist = unique(config)
nconfig = n_elements(configlist)

for c = 1, nconfig do begin  ; note that zero means no known configuration
   print, 'Configuration #', clip(c), ': ', configlist[c-1]
   obs[where(config eq configlist[c-1])].config = c
   print, '   ', obs[where(obs.config eq c)].filename
endfor

biasfile = ''
if obs[0].inst eq 'GMOS-N' then begin
  print, 'Gemini-NORTH'
  biasfile = 'bias_n.fits'
endif 
if obs[0].inst eq 'GMOS-S' then begin
  print, 'Gemini-SOUTH'
  biasfile = 'bias_s.fits'
endif
if biasfile eq '' then begin
  print, 'Warning - not sure if this is North or South'
  print, 'Assuming South...'
  biasfile = 'bias_s.fits'
endif

biasfilefound = file_test(biasfile)
if biasfilefound eq 0 then begin
   print, 'ERROR: Bias file not found!'
   print, 'You need to rename or symlink a processed GMOS bias file as ', biasfile
   return
endif
bias = loadandmerge(biasfile)


wsciall = where(obs.type eq 'OBJECT', ctsciall)
outfiles = strarr(ctsciall)
nsciall = 0

for c = 1, nconfig do begin
   print
   print, 'Now reducing configuration #', clip(c), ' (', configlist[c-1], ')'

   obsc = obs[where(obs.config eq c)]

   wflat = where(obsc.type eq 'FLAT', ctflat)
   warc = where(obsc.type eq 'ARC', ctarc)
   wsci = where(obsc.type eq 'OBJECT', ctsci)

   print, 'Found ', clip(fix(ctflat)), ' flats, ', clip(fix(ctarc)), ' arcs, and ', clip(fix(ctsci)), ' science exposures.'


   if ctsci eq 0 then begin
      print, 'ERROR: No science exposures found!  Skipping reductions for this configuration...'
      break
   endif
   if ctarc eq 0 then begin
      print, 'WARNING: No arcs found!  Plotting in pixels and not wavelengths.  Yeah, good luck with that.'
   endif
   if ctflat eq 0 then begin
      print, 'WARNING: No flat-fields found!   Skipping flat-fielding...'
      flat = 1
   endif

   ; make the normalized flat
   if ctflat gt 0 then begin
      if ctflat gt 1 then flats = fltarr(nx,ny,ctflat)
      for nf = 0, ctflat-1 do begin
         ; bias subtraction
         flat = loadandmerge(obsc[wflat[0]].filename, flath)
         flat = flat - bias
   
         ; response division
         nx = (size(flat)) [1]
         ny = (size(flat)) [2]

         for x = 0, nx-1 do begin
            xmed = 1.0 * median(flat[x,*], /even);this is ignoring all the zeros but for now we're just doing a tentative thing
            flat[x,*] = flat[x,*] / xmed
         endfor     

         ; reverse gain-correct the flat so that later exposures are automatically gain-corrected
         flat[0:a1xmax,*] = flat[0:a1xmax,*] / 2.372
         flat[a2xmin:a2xmax,*] = flat[a2xmin:a2xmax,*] / 2.076
         flat[a3xmin:a3xmax,*] = flat[a3xmin:a3xmax,*] / 2.097

         if ctflat gt 1 then flats[*,*,nf] = flat
      endfor
      if ctflat gt 1 then begin      
         for x = 0, nx-1 do begin
            for y = 0, ny-1 do begin
               flat[x,y] = median(flats[x,y,*], /even)
            endfor
         endfor
      endif
      mwrfits, flat, 'outflat_c'+clip(c)+'.fits', flath, /silent, /create
   endif else begin
      flat = 1
   endelse



   ; now get the science spectrum

   objspec = fltarr(ctsci, nx)
   skyspec = fltarr(ctsci, nx)

   for nsci = 0, ctsci-1 do begin

      if ctsci gt 1 then begin
         print, 'Exposure #', clip(nsci+1), '/', clip(ctsci)
      endif

      sci = loadandmerge(obsc[wsci[nsci]].filename, h)
      sci = sci - bias
      sci = sci / flat  ; includes gain correction

      object = clip(sxpar(h,'OBJECT'))
      pa = sxpar(h,'PA')
      utstart = sxpar(h,'UTSTART')
      utend = sxpar(h,'UTEND')
      slit = sxpar(h,'MASKNAME')
      exptime = sxpar(h,'EXPTIME')
      grating = sxpar(h,'GRATING')
      cwave = sxpar(h,'CENTWAVE')
      airmass = sxpar(h,'AIRMASS')

      print, 'Processing ', obsc[wsci[nsci]].filename

      if (crreject gt 0) then begin
         print, 'Removing cosmic rays...'
         zap1 = quickzap(sci[0:a1xmax,900:1500-1])
         sci[0:a1xmax,900:1500-1] = zap1
         zap2 = quickzap(sci[a2xmin:a2xmax,900:1500-1])
         sci[a2xmin:a2xmax,900:1500-1] = zap2
         zap3 = quickzap(sci[a3xmin:a3xmax,900:1500-1])
         sci[a3xmin:a3xmax,900:1500-1] = zap3

         in2dspecname = obsc[wsci[nsci]].filename
         pos = strpos(in2dspecname,'.')
         if pos gt 1 then begin
            froot = strmid(in2dspecname,0,pos)
            ext = strmid(in2dspecname,pos+1)
            mwrfits, sci, froot+'z'+'.'+ext, h, /silent, /create
         endif
         print, '...done.'
      endif

      ; find the object

      ysearchmin = ysearch[0]
      ysearchmax = ysearch[1]
      wavsum = fltarr(ysearchmax-ysearchmin)
      for y = ysearchmin, ysearchmax-1 do begin  ; this is very CR-nonrobust
         wavsum[y-ysearchmin] = total(sci[a2xmin:a2xmax,y])   ; looking in a2 will not be good for high redshift
      endfor
      kernel = [0.2,1.0,1.0,0.2]
      wavsumconv = convol(wavsum, kernel, /edge_truncate)

      centersigma = max(wavsum)/stdev(wavsum)
      if centersigma lt 3 then begin
         print, 'WARNING:  Cannot find the source!   Best center is only ', centersigma, '-sigma above background.'
         print, 'Continuing anyway, but you will likely want to rerun after specifying ysearch.'
      endif

      avcenty = ysearchmin + (where(wavsumconv eq max(wavsumconv))) [0] ; max val after convolution


      profmax = max(wavsum)
      baseline = median(wavsum, /even)
      for y = avcenty, ysearchmax-1 do begin
         if wavsum[y-ysearchmin] lt baseline + (profmax-baseline)*0.2 then begin
            apmax = y
            break
         endif
      endfor
      for y = avcenty, ysearchmin-1, -1 do begin
         if wavsum[y-ysearchmin] lt baseline + (profmax-baseline)*0.2 then begin
            apmin = y
            break
         endif
      endfor

      if (apmax-avcenty) gt 6 then apmax = avcenty + 5
      if (avcenty-apmin) gt 6 then apmin = avcenty - 5

      print, 'Found object at y=', strtrim(avcenty,2), ' (+', clip(apmax-avcenty), ' -', clip(avcenty-apmin), ')'


      ; trace the object across the chips
   
      ; march across the image, left to right, in blocks of 8, summing up horizontally for every y to get 
      ; the variable spatial profile
      inc = 8 ; this should probably be an adjustable parameter
      centy = intarr(a3xmax/inc)
      centx = intarr(a3xmax/inc)
      scansize = 20 ; how far away from the overall center to look for trace variation
      wavsum = fltarr(scansize*2+1)
      lastgoodtracei = 0
      ntraceloss = 0
      ntracefound = 0
      lastgoodtracex = -1
      lastgoodtracey = avcenty
      for i = 3, (a3xmax/inc)-1 do begin   ; start not quite at the left edge in case this is really far to red
         x1 = i*inc
         x2 = x1+inc-1 < a3xmax
         for y = avcenty-scansize, avcenty+scansize do begin
            wavsum[y-(avcenty-scansize)] = total(sci[x1:x2,y])  ; sum across 8 pix within the block
         endfor
         wavsumconv = convol(wavsum, kernel, /edge_truncate)
         centyi = avcenty-scansize + (where(wavsumconv eq max(wavsumconv))) [0]
         centxi = i*inc + inc/2

         if abs(centyi-lastgoodtracey) gt 3 and i ne 0 then begin
            ntraceloss = ntraceloss + 1
            ntracefound = 0
            if ntraceloss eq 5 then begin ; trace bad for last five blocks: officially lost, mark the position
               tracelossi = i - 5
               print, 'Lost trace at x=', clip(lastgoodtracex)
            endif
         endif else begin
            centy[i] = centyi
            centx[i] = centxi
            ntracefound = ntracefound + 1
         endelse

         if ntraceloss eq 0 then begin
            lastgoodtracex = centxi
            lastgoodtracey = centyi
         endif
            
         if ntracefound eq 3 then begin
            if ntraceloss ge 5 then print, 'Refound trace'
            ntraceloss = 0    ; good 3 times in a row: it's officially back, so all is cool
         endif
      endfor

      centy = centy[where(centy gt 0)]
      centx = centx[where(centy gt 0)]

           ; lastgoodtracei = lastgoodtracei - 1 ; remove one more to be safe
           ; print, 'Lost the trace at x=',clip(lastgoodtracei*inc)
           ; centx = centx[0:tracelossi]
           ; centy = centy[0:tracelossi]
           ; if centx[tracelossi] lt 2800 then $
           ;; print, 'This could be either to poor conditions (faint source, bright moon) or the burst could be at high-z.'
           ; break


      ; fit a polynomial to the profile, excluding possible bad traces

      goodtrace = where(abs(centy-avcenty) lt 6)

      tracepar = poly_fit(centx[goodtrace], centy[goodtrace],2)
      print, 'Trace function: '+ clip(tracepar[0]) + ' + ' + clip(tracepar[1]) + 'x + ' + clip(tracepar[2]) +'x^2'

      tracey = round(poly(findgen(nx), tracepar))
 
      ;av1centy = median(centy[0:a2xmin/inc-1], /even)
      ;av2centy = median(centy[a2xmin/inc:a3xmin/inc-1], /even)
      ;av3centy = median(centy[a3xmin/inc-1:2800/inc], /even) ; don't go quite to the end

      ;mad1 = total(abs(centy[0:a2xmin/inc-1] - av1centy))/(a2xmin/inc)
      ;mad2 = total(abs(centy[a2xmin/inc:a3xmin/inc-1] - av2centy))/(a3xmin/inc-a2xmin/inc)
      ;mad3 = total(abs(centy[a3xmin/inc-1:2800/inc] - av3centy))/(2800/inc-a3xmin/inc)
      ;print, 'Trace deviations: ', clip(mad1) + ', ' + clip(mad2) + ', ' + clip(mad3)

      ;if mad3 gt 14 then begin
      ;   print, 'Having difficulty tracing the object in right chip.  The S/N is probably poor.  Use of this'
      ;   print, 'routine is not advised for low-signal spectra.'
      ;   if mad2 lt 10 then av3centy = av2centy
      ;endif
      ;if mad2 gt 8 and (av3centy - av2centy) gt 6 then begin
      ;   print, 'Having difficulty finding the object in middle chip.  The S/N may be poor or the source may be'
      ;   print, 'at very high redshift.   Using right chip position.'
      ;   av2centy = av3centy
      ;endif
      ;if mad1 gt 8 and (av2centy - av1centy) gt 12 then begin
      ;   print, 'Having difficulty finding the object in left chip.  The S/N may be poor or the source may be'
      ;   print, 'at high redshift.   Using middle chip position'
      ;   av1centy = av2centy
      ;endif

      print, 'Extracting...'

      ; extract the object and the sky
      fobj = fltarr(nx)
      fbak1 = fltarr(nx)
      fbak2 = fltarr(nx)
      apradp = apmax-avcenty
      apradm = avcenty-apmin
      apwidth = apradp+apradm + 1
      skyrad = 11  
      skywidth = 11   ;note confusing terminology here.  ap is -apradm to apradp, sky is skyrad to skyrad+skywidth
      for x = 0, nx-1 do begin
         ;if x lt a2xmin then cy = av1centy
         ;if x ge a2xmin and x lt a3xmin then cy = av2centy
         ;if x ge a3xmin then cy = av3centy
         cy = tracey[x]
         apy1 = cy-apradm
         apy2 = cy+apradp
         s1y1 = cy-skyrad-(skywidth-1)
         s1y2 = cy-skyrad
         s2y1 = cy+skyrad
         s2y2 = cy+skyrad+(skywidth-1)
         fobj[x]  = total(sci[x,apy1:apy2])
         fbak1[x] = total(sci[x,s1y1:s1y2])
         fbak2[x] = total(sci[x,s2y1:s2y2])
      endfor

      bakdiff = fbak1 - fbak2
      avbakdiff = median(bakdiff, /even)
      baksigma = stdev(bakdiff)
      f1cosmic = where(fbak1 gt fbak2 + avbakdiff + 6*baksigma, ctf1c)
      f2cosmic = where(fbak2 gt fbak1 - avbakdiff + 6*baksigma, ctf2c)
      if ctf1c gt 0 then fbak1[f1cosmic] = fbak2[f1cosmic] + avbakdiff
      if ctf2c gt 0 then fbak2[f2cosmic] = fbak1[f2cosmic] - avbakdiff
      baksigma = stdev(bakdiff)
      f1cosmic = where(fbak1 gt fbak2 + avbakdiff + 5*baksigma, ctf1c)
      f2cosmic = where(fbak2 gt fbak1 - avbakdiff + 5*baksigma, ctf2c)
      if ctf1c gt 0 then fbak1[f1cosmic] = fbak2[f1cosmic] + avbakdiff
      if ctf2c gt 0 then fbak2[f2cosmic] = fbak1[f2cosmic] - avbakdiff

      ;res1 = poly_fit(centx[0:a2xmin/inc-1]         , centy[0:a2xmin/inc-1], 2)
      ;res2 = poly_fit(centx[a2xmin/inc:a3xmin/inc-1], centy[a2xmin/inc:a3xmin/inc-1], 2)
      ;res3 = poly_fit(centx[a2xmin/inc-1:*]         , centy[a2xmin/inc-1:*], 2)

      objspec[nsci,*] = fobj - (((fbak1+fbak2)/2)*apwidth/skywidth)
      skyspec[nsci,*] = ((fbak1+fbak2)/2)*apwidth/skywidth

      print, 'Extraction complete.'

   endfor

   !p.multi = [0, 2, 3]
   psopen, 'arcplots_c'+clip(c)+'.ps'

   ; Figure out the wavelength calibration

   if ctarc gt 0 then begin

      arc = loadandmerge(obsc[warc[0]].filename, arch) ; only examine the first arc if there is more than one
      arc = arc - bias
      arc = arc / flat

      grating = sxpar(arch,'GRATING') ; should be same as original grating since same config
      grat = strmid(grating,0,4)
      cwav = fix(sxpar(arch,'CENTWAVE'))

      ; extract arc spectrum
      print, 'Extracting arc spectrum.'
      farc = fltarr(nx)
      for x = 0, nx-1 do begin
         ;if x lt a2xmin then cy = av1centy                       ; Note that this is doing the lazy thing 
         ;if x ge a2xmin and x lt a3xmin then cy = av2centy       ; and using whatever aperture was last
         ;if x ge a3xmin then cy = av3centy                       ; extracted.  To be marginally more precise
         cy = tracey[x]                                           ; one could remember all apertures and loop.
         apy1 = cy-apradm                                       
         apy2 = cy+apradp
         farc[x]  = total(arc[x,apy1:apy2])/(apy2-apy1+1.0)
      endfor
      ; do a really rough background continuum subtraction
      fbakconv = convol(farc, [0.1,0.2,0.2,0.2,0.1], /edge_truncate)
      for x = 0, nx-1 do begin
         farc[x] = farc[x] - min(fbakconv[x-25>1:x+25<nx-1])
      endfor

      farcx = findgen(n_elements(farc))
      arckernel = [0.2,0.4,0.9,1,0.9,0.3,0.2]
      arckernel = arckernel / total(arckernel)
      arcconv = convol(farc, arckernel, /edge_truncate)
      ; identify lines
      arcthreshold = 55 ; this is a parameter that should probably be adjustable at top level
      linex = fltarr(1000) ; up to 1000 lines...
      linestrength = fltarr(1000)
      l = 0
      for x = 5, nx-5 do begin
         if arcconv[x] lt arcthreshold then continue
         if arcconv[x] gt arcconv[x-1] and arcconv[x] gt arcconv[x+1] then begin
            ; get a more precise centroid
            ; linex[l] = x ;

            linex[l] = total(farc[x-3:x+3] * farcx[x-3:x+3])/total(farc[x-3:x+3])
            linestrength[l] = arcconv[x]
            l = l + 1
         endif
      endfor
      if l gt 9 then begin
         linex = linex[0:l-1]   
         print, 'Found ', strtrim(l,2), ' arc lines.'
      endif else begin
         print, 'Only found ', strtrim(l, 2), ' arc lines!   Not enough to calibrate reliably...'
         break  ; should just rerun with a lower threshold
      endelse

      brightarcfile = 'CuAr_GMOS_bright.dat'
      allarcfile = 'CuAr_GMOS.dat'
      checkarclist = file_test(brightarcfile)
      checkallarclist = file_test(allarcfile)
      if checkarclist eq 0 or checkallarclist eq 0 then begin
         print, 'Cannot find arc line list!'
         print, 'Must be titled ', brightarcfile, ' and ', allarcfile, '.'
         return
      endif
      brightreflines = grabcolumn(brightarcfile,0)
      reflines = grabcolumn(allarcfile,0)

      if grat eq 'R400' then begin
        maxsep = 1024*1.34
        dispersionrange = [-1.25, -1.45]
        guesswav0 = [10157.8, 10110.5, 10036.1] + 10*(cwav-800)
        guessdispersion = [-1.395, -1.376, -1.353]
        wav0range1 = [7950,8300] + 10*(cwav-600)
        wav0range2 = [7900,8250] + 10*(cwav-600)
        wav0range3 = [7850,8200] + 10*(cwav-600)
      endif else if grat eq 'B600' then begin
        maxsep = 1024*0.95
        dispersionrange = [-0.8, -1.0]
        guesswav0 = [7440.761, 7409.889, 7360.568] + 10*(cwav-600)
        guessdispersion = [-0.9345, -0.9232, -0.9087]
        wav0range1 = [7250,7650] + 10*(cwav-600)
        wav0range2 = [7200,7600] + 10*(cwav-600)
        wav0range3 = [7150,7550] + 10*(cwav-600)
      endif else if grat eq 'R600' then begin         ; remaining gratings have not been tested
        maxsep = 1024*2*0.47
        dispersionrange = [-0.8, -1.05]
        guesswav0 = replicate(cwav*10 + 2860/2, 3)
        guessdispersion = replicate(-2*0.47,3)
        wav0range1 = [guesswav0[0]-500*abs(guessdispersion[0]),guesswav0[0]+500*abs(guessdispersion[0])]
        wav0range2 = [guesswav0[1]-500*abs(guessdispersion[1]),guesswav0[1]+500*abs(guessdispersion[1])]
        wav0range3 = [guesswav0[2]-500*abs(guessdispersion[2]),guesswav0[2]+500*abs(guessdispersion[2])]
      endif else if grat eq 'R150' then begin
        maxsep = 1024*2*1.74
        dispersionrange = [-3, -4]
        guesswav0 = replicate(cwav*10 + 10710/2, 3)
        guessdispersion = replicate(-2*1.74,3)
        wav0range1 = [guesswav0[0]-500*abs(guessdispersion[0]),guesswav0[0]+500*abs(guessdispersion[0])]
        wav0range2 = [guesswav0[1]-500*abs(guessdispersion[1]),guesswav0[1]+500*abs(guessdispersion[1])]
        wav0range3 = [guesswav0[2]-500*abs(guessdispersion[2]),guesswav0[2]+500*abs(guessdispersion[2])]
      endif else if grat eq 'B1200' then begin
        maxsep = 1024*2*0.23
        dispersionrange = [-0.6, -0.75]
        guesswav0 = replicate(cwav*10 + 1530/2, 3)
        guessdispersion = replicate(-2*0.34,3)
        wav0range1 = [guesswav0[0]-500*abs(guessdispersion[0]),guesswav0[0]+500*abs(guessdispersion[0])]
        wav0range2 = [guesswav0[1]-500*abs(guessdispersion[1]),guesswav0[1]+500*abs(guessdispersion[1])]
        wav0range3 = [guesswav0[2]-500*abs(guessdispersion[2]),guesswav0[2]+500*abs(guessdispersion[2])]
      endif else if grat eq 'R831' then begin
        maxsep = 1024*2*0.34
        dispersionrange = [-0.65, -0.75]
        guesswav0 = replicate(cwav*10 + 2070/2, 3)
        guessdispersion = replicate(-2*0.23,3)
        wav0range1 = [guesswav0[0]-500*abs(guessdispersion[0]),guesswav0[0]+500*abs(guessdispersion[0])]
        wav0range2 = [guesswav0[1]-500*abs(guessdispersion[1]),guesswav0[1]+500*abs(guessdispersion[1])]
        wav0range3 = [guesswav0[2]-500*abs(guessdispersion[2]),guesswav0[2]+500*abs(guessdispersion[2])]
      endif else begin
        print, 'WARNING: Unrecognized grating.   It is very unlikely the wavelength solution will be accurate.'
        maxsep = 1024*4.0
        dispersionrange = [-0.2, -2]
        guessdispersion = replicate(-1,3)
        wav0range1 = [5000,12000]
        wav0range2 = [5000,12000]
        wav0range3 = [5000,12000]
      endelse
      
      ; left chip
      a1line_w = where(linex le a2xmin, ctla1)
      if ctla1 lt 3 then begin
         print, 'WARNING: Only ', clip(ctla1), ' arc line(s) on left chip!  This is not enough to calibrate.'
         print, '         Using approximate archival solution.'
         a1sol = [guesswav0[0], guessdispersion[0],0]
      endif else begin
         print, clip(ctla1), ' arc lines on left chip.  Determining wavelength solution...'
         a1line_s = reverse(sort(linestrength[a1line_w]))
         if ctla1 gt 9 then a1line_s = a1line_s[0:9]
         a1brightlinesx = linex[a1line_w[a1line_s]]      
         
         a1sol = solvearc(linex[a1line_w], a1brightlinesx, reflines, brightreflines, wav0range=wav0range1, disprange=dispersionrange, maxsep=maxsep)
         a1wavrange = poly([0,a1xmax],a1sol)
         cdisp1 = a1sol[1] + 2 * a1sol[2] * (a1xmax/2)
         print, 'Wavelength range for left chip:   ', clip(a1wavrange[0]), '--', clip(a1wavrange[1]), $
                ' (central dispersion ', clip(cdisp1), ' A/pix)' 
         ;print, 'Wavelength solution for left chip:   lambda = ', clip(a1sol[0]), ' - ', clip(-a1sol[1]), '*pix'
      endelse


      ; middle chip
      a2line_w = where(linex ge a2xmin and linex lt a3xmin, ctla2)
      if ctla2 lt 3 then begin
         print, 'WARNING: Only ', ctla2, ' line(s) on middle chip!  This is not enough to calibrate.'
         print, '         Using approximate archival solution.'
         a2sol = [guesswav0[1], guessdispersion[1],0]
      endif else begin
         print, clip(ctla2), ' arc lines on middle chip.  Determining wavelength solution...'
         a2line_s = reverse(sort(linestrength[a2line_w]))
         if ctla2 gt 9 then a2line_s = a2line_s[0:9]
         a2brightlinesx = linex[a2line_w[a2line_s]]      
         
         a2sol = solvearc(linex[a2line_w], a2brightlinesx, reflines, brightreflines, wav0range=wav0range2, disprange=dispersionrange, maxsep=maxsep)
      endelse
      a2wavrange = poly([a2xmin,a2xmax],a2sol)
      cdisp2 = a2sol[1] + 2 * a2sol[2] * (a2xmin + (a2xmax-a2xmin)/2)
      print, 'Wavelength range for middle chip: ', clip(a2wavrange[0]), '--', clip(a2wavrange[1]), $
                ' (central dispersion ', clip(cdisp2), ' A/pix)' 
      ;print, 'Wavelength solution for middle chip: lambda = ', clip(a2sol[0]), ' - ', clip(-a2sol[1]), '*pix'


      ; right chip
      a3line_w = where(linex ge a3xmin, ctla3)
      if ctla3 lt 3 then begin
         print, 'WARNING: Only ', ctla3, ' line(s) on right chip!  This is not enough to calibrate.'
         print, '         Using approximate archival solution.'
         a3sol = [guesswav0[2], guessdispersion[2],0]
      endif else begin
         print, clip(ctla3), ' arc lines on right chip.  Determining wavelength solution...'
         a3line_s = reverse(sort(linestrength[a3line_w]))
         if ctla3 gt 9 then a3line_s = a3line_s[0:9]
         a3brightlinesx = linex[a3line_w[a3line_s]]      
         
         a3sol = solvearc(linex[a3line_w], a3brightlinesx, reflines, brightreflines, wav0range=wav0range3, disprange=dispersionrange, maxsep=maxsep)
      endelse
      a3wavrange = poly([a3xmin,a3xmax],a3sol)
      cdisp3 = a3sol[1] + 2 * a3sol[2] * (a3xmin + (a3xmax-a3xmin)/2)
      print, 'Wavelength range for right chip:  ', clip(a3wavrange[0]), '--', clip(a3wavrange[1]), $
                ' (central dispersion ', clip(cdisp3), ' A/pix)' 
      ;print, 'Wavelength solution for right chip:  lambda = ', clip(a3sol[0]), ' - ', clip(-a3sol[1]), '*pix'

      ;print, 'Central dispersions: ', cdisp1, cdisp2, cdisp3

      if abs(cdisp2-cdisp1) gt 0.08 then print, 'WARNING: wavelength solutions appear to be inconsistent'
      if abs(cdisp3-cdisp2) gt 0.08 then print, 'WARNING: wavelength solutions appear to be inconsistent'

      ;if abs(a2sol[1]-a1sol[1]) gt 0.08 then print, 'WARNING: wavelength solutions appear to be inconsistent'
      ;if abs(a3sol[1]-a2sol[1]) gt 0.08 then print, 'WARNING: wavelength solutions appear to be inconsistent'

      pix = findgen(nx)
      wav = fltarr(nx)

      wav[0:a1xmax] = poly(pix[0:a1xmax],a1sol)
      wav[a2xmin:a2xmax] = poly(pix[a2xmin:a2xmax],a2sol)
      wav[a3xmin:nx-1] = poly(pix[a3xmin:nx-1],a3sol)

   endif

   psclose
   !p.multi = 0


   ; Remove points due to known bad columns: 1025, 1229, 1230, 1569, 1979, 2755-57, first and last (all-1)

   x = indgen(nx)
   badcol = [1,1025, 1229, 1230, 1569, 1979, 2755, 2756, 2757, n_elements(x)]-1
   x[badcol] = -1
   nbadcol = n_elements(badcol)
   goodcol = where(x ge 0)
   x = x[goodcol]
   wav = wav[goodcol]
   objspec = objspec[*,goodcol]
   skyspec = skyspec[*,goodcol]

   ; Save to disk

   for nsci = 0, ctsci-1 do begin
      outfile =  strlowcase(object)+'_'+clip(grating,4)+'c'+clip(fix(cwave))+'_'+clip(nsci+1)+'.dat'
      outfiles[nsciall] = outfile
      nsciall = nsciall + 1
      openw, 1, outfile
      printf, 1, '# OBJECT   = ', object
      printf, 1, '# GRATING  = ', grating
      printf, 1, '# CWAVE    = ', cwave
      printf, 1, '# NEXP     = ', nsci
      for i = 0, n_elements(wav)-1 do begin
          printf, 1, wav[i], objspec[nsci,i], skyspec[nsci,i], x[i]
      endfor
      close, 1
      print, 'Spectrum saved to ', outfile
   endfor

endfor

print
print, 'Reduction of 2D frames complete.   Proceeding to analysis.'
print

;save, /variables, filename='save.sav'
;resumehere:
;restore, 'save.sav'

ctsciall = nsciall
ctsci = ctsciall ; overwrite these parameters (lazy)

prevconfig = ''

nc = nx-nbadcol
allwav   = fltarr(ctsci,nc)
allspec  = fltarr(ctsci,nc)
allsky   = fltarr(ctsci,nc)
allnorm  = fltarr(ctsci,nc)
object = strarr(ctsci)
cwave = strarr(ctsci)
grating = strarr(ctsci)
nexp = strarr(ctsci)

for nsci = 0, ctsci-1 do begin

   infile = outfiles[nsci]

   allwav[nsci,*] =  grabcolumn(infile, 0)
   allspec[nsci,*] = grabcolumn(infile, 1)
   allsky[nsci,*] = grabcolumn(infile, 2)
   iline = ''
   openr, 1, infile
   for i = 0, 5 do begin
     readf, 1, iline
     if strmid(iline, 0, 1) eq '#' then begin
        idata = strtrim(strsplit(iline, /extract),2)
        if n_elements(idata) lt 4 then continue
        if idata[1] eq 'OBJECT'  then object[nsci] = idata[3]
        if idata[1] eq 'GRATING' then grating[nsci] = idata[3]
        if idata[1] eq 'CWAVE'   then cwave[nsci] = idata[3]
        if idata[1] eq 'NEXP'    then nexp[nsci] = idata[3]
     endif
   endfor
   close, 1

endfor

   ; Line detection

   maxlines = 200
   linelists = fltarr(ctsci,maxlines+1)
   nliness = intarr(ctsci)
   alllines = fltarr(maxlines,nc) 

   for nsci = 0, ctsci-1 do begin
      wav  = (allwav[nsci,*]) [*]
      spec = (allspec[nsci,*]) [*]
      normfitx = 5+5*indgen(nc/5-2)
      normfitv = fltarr(nc/5-2)
      for b = 0, n_elements(normfitx)-1 do begin
        nearby = spec[normfitx[b]-100>0:normfitx[b]+100<nc-1]
        goodregions = where(nearby gt 0)
        for i = 0, 5 do begin
          av = median(nearby[goodregions], /even)
          sd = stdev(nearby[goodregions])
          goodregions = where(nearby gt av - 3.3*sd and nearby lt av + 4*sd)
        endfor
        normfitv[b] = median(nearby[goodregions], /even)
      endfor
      fit = poly_fit(wav[normfitx], normfitv, 6)
      yfit = poly(wav, fit)
      normspec = spec/yfit
      allnorm[nsci,*] = normspec

      kernel = [0.02,0.04,0.1,0.5,0.8,0.9,0.95,0.98,1.0,1.0,1.0,0.98,0.95,0.9,0.8,0.5,0.1,0.04,0.02]
      kernel = kernel / total(kernel)
      convnormspec = convol(normspec, kernel)
      ;plot, wav, convnormspec
      lineregions = where(convnormspec lt 0.75 or convnormspec gt 2.0, ct)
      linelist = fltarr(maxlines+1)
      nline = 0
      if ct gt 0 then begin
         lineregionwav = wav[lineregions]
         for i = 0, n_elements(lineregions)-1 do begin   ; any sane programming language could do this block
            if i eq 0 then begin                         ; in 1/10 the code that IDL demands ><
               starti = i
            endif else begin
               if lineregions[i] ne lineregions[i-1]+1 then starti = i
            endelse
            endthelinehere = 0
            if i eq (n_elements(lineregions)-1) then begin 
               endthelinehere = 1
            endif else begin
               if lineregions[i] ne lineregions[i+1]-1 then endthelinehere = 1
            endelse
            if endthelinehere eq 1 then begin   
               endi = i
               thiswavregion = starti + intarr(endi-starti+1)
               linecenter = lineregionwav[min(i)]
               linelist[nline] = linecenter
               nline = nline + 1
            endif
            if nline ge maxlines then break
         endfor
         linelists[nsci,*] = linelist
         linelist = linelist[0:nline-1]
         nliness[nsci] = nline
      endif
   endfor

   ; Line matching / printout

   linelists = transpose(linelists[*,0:max(nliness)])
   linelists[where(linelists eq 0)] = 99999
   compthreshold = 10

   ijk = intarr(ctsci)
   offset = indgen(ctsci)*((size(linelists)) [1])
   ctlines = 0
   ctmultlines = 0
   openw, 1, 'detectedlines.dat'
   printf, 1, 'Possible lines:'
   if ctsci eq 1 then begin
      for a = 0, n_elements(linelists)-1 do printf, 1, clip(linelists[a])
   endif else begin
    for a = 0, n_elements(linelists)-1 do begin
      match = intarr(ctsci)
      aaa = ijk + offset
      unfinished = where(linelists[aaa] ne 99999, ct)
      if ct eq 0 then break
      currentlist = (where(linelists[aaa] eq max(linelists[aaa[unfinished]]) )) [0]
      match[currentlist] = 1
      otherlist = where(indgen(ctsci) ne currentlist)
      for o = 0, n_elements(otherlist)-1 do begin
        if abs(linelists[aaa[otherlist[o]]]-linelists[aaa[currentlist]]) lt compthreshold then begin
           match[otherlist[o]] = 1
        endif
      endfor
      outstr = ''
      for c = 0, ctsci-1 do begin
        if match[c] eq 1 then begin
          outstr = outstr + clip(linelists[aaa[c]],10)
          ijk[c] = ijk[c] + 1
        endif else begin
          outstr = outstr + clip('',10)
        endelse
      endfor

      if total(match) ge 1 then ctlines = ctlines + 1
      if total(match) ge 2 then ctmultlines= ctmultlines + 1
      printf, 1, outstr
    endfor
   endelse
   close, 1
   print, 'Detected ', clip(ctlines), ' line candidates (', clip(ctmultlines), ' multiple detections)'
   print, 'List of possible lines written to detectedlines.dat'

   ; Plotting

   ; open the plot
     plotfilename = 'specplots.ps'
     xs = 6
     ys = 4.5
     initdev = !d.name
     openplot, plotfilename, xs, ys


   prevgrating = ''
   col = 2

   for nsci = 0, ctsci-1 do begin

      if nsci eq 0 or grating[nsci] ne prevgrating then begin
         sspec = allspec[nsci,sort(allspec[nsci,*])]
         yr = [0,1.05*sspec[fix((n_elements(sspec)-1)*0.99)]] ; set max to 5% above the 99 percentile flux
         xr = [max(allwav),min(allwav)]

         plot, [0], [0], xrange=xr, xstyle=1, yrange=yr, ystyle=1, xtitle='Wavelength', ytitle='Flux (counts)', $
                         title = object[nsci] + ' ('+clip(grating[nsci],4)+')', /nodata
         
         ; telluric lines
         ybox = [yr[0]+yr[1]/1000.,yr[0]+yr[1]/1000.,yr[1]*0.999,yr[1]*0.999,yr[0]+yr[1]/1000.]
         if min(xr) lt 7594 and max(xr) gt 7621 then $
            polyfill, [7594, 7621+35, 7621+35, 7594, 7594], ybox, color=1
         if min(xr) lt 6884 and max(xr) gt 6867 then $
            polyfill, [6867, 6884+10, 6884+10, 6867, 6867], ybox, color=1
         plot, [0], [0], xrange=xr, xstyle=1, yrange=yr, ystyle=1, /nodata, /noerase

         refscale = nsci

         ;xyouts, 0.2, 0.86, /norm, object[nsci]
         ;xyouts, 0.2, 0.8, /norm, strmid(grating[nsci],0,4)

      endif else begin
         col = col ;; + 1
      endelse

         ;other lines ???

      ;
      rescale = total(allspec[refscale,*])/total(allspec[nsci,*])
      oplot, allwav[nsci,*], allspec[nsci,*]*rescale, color=col
      col = col + 1

      if ctsci eq 1 then begin
         spec = (allspec[nsci,*]) [*]
         smoothspec = convol(spec, [0.075,0.2,0.45,0.2,0.075], /edge_truncate)
         oplot, allwav[nsci,5:nc-5], smoothspec[5:nc-5], psym=10
         break
      endif

      prevgrating = grating[nsci]
      if nsci ne ctsci-1 then nextgrating = grating[nsci+1]
      if nsci ne ctsci-1 then nextcwave = cwave[nsci+1]
      if nsci eq ctsci-1 or grating[nsci] ne nextgrating or cwave[nsci] ne nextcwave then begin
         avspec = fltarr(nc)
         for x = 0, nc-1 do begin
            avspec[x] = total(allspec[where(grating eq grating[nsci] and cwave eq cwave[nsci]),x])
         endfor
         avspec = avspec * total(allspec[refscale,*])/total(avspec)
         smoothspec = convol(avspec, [0.075,0.2,0.45,0.2,0.075], /edge_truncate)
         oplot, allwav[nsci,5:nc-5], smoothspec[5:nc-5]
      endif
   endfor

   psclose
   print, 'Reduced spectra plotted to ', plotfilename 

     plotfilename = 'specplotsnorm.ps'
     xs = 6
     ys = 4.5
     initdev = !d.name
     openplot, plotfilename, xs, ys

     col = 2

     sspec = allnorm[sort(allnorm)]
     yr = [0,min([1.15*sspec[fix((n_elements(sspec)-1)*0.992)],1.3])] ; set max to 15% above the 99.2 percentile flux
     xr = [max(allwav),min(allwav)]
     plot, [0], [0], xrange=xr, xstyle=1, yrange=yr, ystyle=1, xtitle='Wavelength', ytitle='Normalized flux', $
                     title = object[0], /nodata
     ybox = [yr[0]+yr[1]/1000.,yr[0]+yr[1]/1000.,yr[1]*0.999,yr[1]*0.999,yr[0]+yr[1]/1000.]
     if min(xr) lt 7594 and max(xr) gt 7621 then $
        polyfill, [7594, 7621+35, 7621+35, 7594, 7594], ybox, color=1
     if min(xr) lt 6884 and max(xr) gt 6867 then $
        polyfill, [6867, 6884+10, 6884+10, 6867, 6867], ybox, color=1
     plot, [0], [0], xrange=xr, xstyle=1, yrange=yr, ystyle=1, /nodata, /noerase

     for nsci = 0, ctsci-1 do begin
        spec = (allnorm[nsci,*]) [*]
        smoothspec = convol(spec, [0.075,0.2,0.45,0.2,0.075], /edge_truncate)
        oplot, allwav[nsci, 5:nc-5], smoothspec[5:nc-5], color=col, psym=10
        col = col + 1
     endfor

     minwavall = fix(min(allwav))+1
     maxwavall = fix(max(allwav))-1
     nwav = fix(maxwavall-minwavall)
     interpwav = findgen(nwav) + minwavall ; rebin to 1 A-pix
     interpfluxall = fltarr(ctsci,nwav)
   
     for nsci = 0, ctsci-1 do begin
        minwav = min(allwav[nsci,*])
        maxwav = max(allwav[nsci,*])
        interpok = where(interpwav gt minwav and interpwav lt maxwav)
        interpfluxall[nsci,interpok] = interpol(allnorm[nsci,*], allwav[nsci,*], interpwav[interpok])
     endfor
     
     interpflux = fltarr(nwav)
     for x = 0, nwav-1 do begin
        fluxknown = where(interpfluxall[*,x] ne 0, ct)
        if ct gt 0 then interpflux[x] = median(interpfluxall[fluxknown,x], /even)
     endfor

     kernel = [0.1,0.2,0.6,0.9,1.0,1.0,1.0,0.9,0.6,0.2,0.1] ; ideally use aperture size
     kernel = kernel/total(kernel)
     interpfluxsmooth = convol(interpflux, kernel, /edge_truncate)

     oplot, interpwav, interpfluxsmooth

   


   psclose
   print, 'Normalized spectra plotted to ', plotfilename


set_plot, initdev
print
print, 'All done.'
print

end