;+
; NAME:
;       ISKYSELECT
;
; PURPOSE:
;       Select sky apertures for ISKYSUBTRACT2D.
;
; CALLING SEQUENCE:
;       iskyselect, skylist=, skyfile=, datapath=, skyapfile=, $
;          skyfirst=, skylast=, norder_sky=, navg_sky=, /skyremember
;
; INPUTS:
;       skylist - input list of FITS files
;       skyfile - scalar or vector text file containing SKYLIST as
;                 generated by ISPEC_MAKELISTS; if SKYFILE is a vector
;                 and SKYAPFILE is defined then they must have the
;                 same number of elements
;
; OPTIONAL INPUTS:
;       datapath   - I/O path
;       skyapfile  - output file name to write the selected sky
;                    apertures (default 'iskyselect.out'); also see
;                    SKYFILE notes
;       skyfirst   - scalar integer index number; only select apertures
;                    for SKYLIST[SKYFIRST:SKYLAST]; useful for testing
;                    on a subset of objects
;       skylast    - see SKYFIRST
;       norder_sky - default order of the sky fit along the spatial
;                    dimension (default 2)
;       navg_sky   - number of columns to average (around center of
;                    chip) for plot
;
; KEYWORD PARAMETERS:
;       skyremember - initialize sky windows to those selected for
;                     previous object 
;
; OUTPUTS:
;       See SKYAPFILE.
;
; OPTIONAL OUTPUTS:
;
; PROCEDURES USED:
;       READCOL, IFORAGE(), SPLOG, RD2DSPEC(), DJS_MEDIAN(), LEGEND,
;       DJS_PLOT, DJS_OPLOT, DJS_ITERSTAT, POLY_ITER, ICLEANUP
;
; COMMENTS:
;       This routine should be run after initial reductions (CCDPROC)
;       to generate the sky aperture file passed to ISKYSUBTRACT2D.  
;
;       If both SKYLIST and SKYFILE are given then use SKYLIST.  At
;       least one of SKYLIST or SKYFILE must be passed.
;
; EXAMPLES:
;
; MODIFICATION HISTORY:
;       J. Moustakas, 2004 Jun 23, U of A, nearly all of this code was 
;          written by Andy Marble; I have excised the interactive
;          selection of sky apertures from ISKYSUBTRACT2D into this
;          routine
;       jm04jul04uofa - added NORDER_SKY optional input
;
; Copyright (C) 2004, John Moustakas
; 
; This program is free software; you can redistribute it and/or modify 
; it under the terms of the GNU General Public License as published by 
; the Free Software Foundation; either version 2 of the License, or
; (at your option) any later version. 
; 
; This program is distributed in the hope that it will be useful, but 
; WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
; General Public License for more details. 
;-

pro iskyselect, skylist=skylist, skyfile=skyfile, datapath=datapath, $
  skyapfile=skyapfile, skyfirst=skyfirst, skylast=skylast, $
  norder_sky=norder_sky, navg_sky=navg_sky, skyremember=skyremember

    nskylist = n_elements(skylist)
    nskyfile = n_elements(skyfile)

    if (nskylist eq 0L) and (nskyfile eq 0L) then begin
       print, 'Syntax - iskyselect, skylist=, skyfile=, datapath=, skyapfile=, $'
       print, '   skyfirst=, skylast=, norder_sky=, navg_sky=, /skyremember'
       return
    endif

    if n_elements(datapath) eq 0L then datapath = cwd()
    if n_elements(skyapfile) eq 0L then skyapfile = 'iskyselect.out'

    if n_elements(norder_sky) eq 0L then norder_sky = 2L
    if n_elements(navg_sky) eq 0L then navg_sky = 20L
    if odd(navg_sky) then navg_sky = navg_sky + 1L

    if (nskylist eq 0L) then begin ; if SKYLIST is not given then read SKYFILE
       
; if SKYFILE is a vector then call this routine recurvsively    
       
       if (nskyfile gt 1L) then begin

          if (nskyfile ne n_elements(skyapfile)) then begin
             splog, 'SKYFILE and SKYAPFILE must have the same number of elements.'
             return
          endif

          for k = 0L, nskyfile-1L do iskyselect, skyfile=skyfile[k], datapath=datapath, $
            skyapfile=skyapfile[k], skyfirst=skyfirst, skylast=skylast, norder_sky=norder_sky, $
            navg_sky=navg_sky, skyremember=skyremember
          return
          
       endif

       if (file_test(datapath+skyfile,/regular) eq 0L) then begin
          splog, 'Sky list file '+datapath+skyfile+' not found.'
          return
       endif else begin
          splog, 'Reading sky list file '+datapath+skyfile+'.'
          readcol, datapath+skyfile, skylist, /silent, format='A', $
            comment='#', delimiter=' '
       endelse

    endif

; optionally crop the list of objects    
    
    if (n_elements(skyfirst) eq 0L) then skyfirst = 0L else skyfirst = skyfirst > 0L
    if (n_elements(skylast) eq 0L) then skylast = n_elements(skylist)-1L else $
      skylast = skylast < (n_elements(skylist)-1L)

    if (skyfirst gt n_elements(skylist)-1L) then begin
       splog, 'SKYFIRST exceeds the number of objects in SKYLIST!'
       return
    endif
    
    skylist = skylist[skyfirst:skylast]
    nimage = n_elements(skylist)

; forage header information

    forage = iforage(datapath+skylist)
    if (size(forage,/type) ne 8L) then begin
       splog, 'Error in IFORAGE().'
       return
    endif

    skyfileobject = strcompress(forage.object,/remove_all)
    len1 = strcompress(max(strlen(skylist)),/remove_all)
    len2 = strcompress(max(strlen(skyfileobject)),/remove_all)
    
; do not overwrite SKYAPFILE    
    
    if file_test(datapath+skyapfile,/regular) then begin
       splog, 'Overwrite existing file '+datapath+skyapfile+' [Y/N]? ', format='(A,$ )'
       cc = get_kbrd(1)
       if strupcase(cc) ne 'Y' then return
    endif

    openw, lun, datapath+skyapfile, /get_lun
    printf, lun, '# Sky apertures file list (ISKYSELECT).'
    printf, lun, '# Generated by ISKYSELECT ('+im_today()+').'
    printf, lun, '# Do not modify this line: REFZERO=1 (all apertures are zero-referenced).'
    printf, lun, '# Columns: FILE, OBJECT, SKY POSITION(s) 1-N, SKY APERTURE(s) 1-N, NORDER_SKY'
    printf, lun, '#'

; loop on each object
    
    window, 0, xsize=550, ysize=550

    for i = 0L, nimage-1L do begin

       cube = rd2dspec(skylist[i],datapath=datapath)
       image = cube.image
       mask = cube.mask
       sigmamap = cube.sigmamap
       invvar = (1.0/sigmamap^2.0) ; * (mask eq 0)
       header = *cube.header
       
       ncols = cube.naxis1
       nrows = cube.naxis2
       rowaxis = findgen(nrows)
       colaxis = findgen(ncols)

       file = strmid(skylist[i],0,strpos(skylist[i],'.fits'))
       
; initialize arrays used below

       if (n_elements(selected) eq 0L) then selected = lonarr(nrows)
       if (n_elements(skyremember) eq 0L) then selected = lonarr(nrows)

       text = "'b'"
       text2 = "'q'"
       text3 = "'f'"

       mdn1 = djs_median(image[0:ncols*0.1,*], 1)
       mdn2 = djs_median(image[0.45*ncols:0.55*ncols,*], 1)
       mdn3 = djs_median(image[0.9*ncols:ncols-1,*], 1)
       djs_iterstat, mdn2, mean=mn, sigma=sg

       select:

       vals = [mdn1, mdn2, mdn3]

       yrange1 = [min(vals)-0.15*(max(vals)-min(vals)), $
         max(vals)+0.15*(max(vals)-min(vals))]
       yrange2 = mn + sg * [-5.0, 5.0]

       plot, rowaxis, mdn2, psym=10, thick=2, xthick=2, ythick=2, $
         yrange=yrange2, ystyle=1, xstyle=3, charsize=1.5, charthick=2, $
         xtitle='row (pixels)', position=[0.2,0.1,0.95,0.35], /normal
       oplot, rowaxis, mdn1, psym=10, thick=2, line=1
       oplot, rowaxis, mdn3, psym=10, thick=2, line=1

       plot, rowaxis, mdn2, /noerase, psym=10, thick=2, xthick=2, ythick=2, $
         yrange=yrange1, ystyle=1, xstyle=3, charsize=1.5, charthick=2, $
         xtickname=replicate(' ', 30), $
         title='column median (near chip center & ends)', $
         position=[0.2,0.35,0.95,0.925], /normal
       oplot, rowaxis, mdn1, psym=10, thick=2, line=1
       oplot, rowaxis, mdn3, psym=10, thick=2, line=1
       legend, /top, /left, box=0, file+ ' ('+skyfileobject[i]+')', charsize=1.5, charthick=2

       xyouts, 0.05, 0.5, align=0.5, orientation=90, charthick=2, charsize=1.5, 'counts', /normal

       whsel = where(selected eq 1L, count)
       if count gt 0L then djs_oplot, psym=4, thick=2, color='green', $
         whsel, whsel*0+(min(vals)-0.065*(max(vals)-min(vals))), symsize=0.5
       
       done = 0L
       nsky = 0L

       while not done do begin

;         junk = dialog_message('Sky window selection: Position cursor and press '+$
;           'any mouse button (once to mark, twice to quit).',/information)
          
          print, format='(" ISKYSELECT: Sky window selection: Position cursor and press' + $
            ' any mouse button (once to mark, twice to quit).   ",A1,$)', string(13b)
          cursor, xmouse1, ymouse, /down
          xmouse1 = (round(xmouse1) > 0L) < (nrows-1)
          djs_oplot, psym=4, color='green', symsize=0.5, xmouse1*[1,1], $
            (min(vals)-0.065*(max(vals)-min(vals)))*[1,1]
          
          print, format='(" Sky window selection: Reposition and press: lef' + $
            't-select, right-unselect, middle-reverse selection.",A1,$)', string(13b)
          
          cursor, xmouse2, ymouse, /down
          xmouse2 = (round(xmouse2) > 0L) < (nrows-1)
          if xmouse2 ne xmouse1 then begin
             wh = where(rowaxis ge min([xmouse1, xmouse2]) and $
               rowaxis le max([xmouse1, xmouse2]), count)
             if count gt 0 then begin
                if !mouse.button eq 1L then selected[wh] = 1L $
                else if !mouse.button eq 4L then selected[wh] = 0L $
                     else selected[wh] = 1L - selected[wh]
             endif
             whsel = where(selected eq 1L, count, comp=whnotsel)
             if count lt nrows then oplot, psym=4, thick=2, color=0, $
               whnotsel, whnotsel*0+(min(vals)-0.065*(max(vals)-min(vals))), symsize=0.5
             if count gt 0 then djs_oplot, psym=4, thick=2, color='green', $
               whsel, whsel*0+(min(vals)-0.065*(max(vals)-min(vals))), symsize=0.5

          endif else begin

             if selected[xmouse1] ne 1L then $
               djs_oplot, psym=4, color='black', symsize=0.5, xmouse1*[1,1], $
                 (min(vals)-0.065*(max(vals)-min(vals)))*[1,1]

             if total(selected) gt 0L then begin

                done = 1L 

                skyrows = where(selected eq 1L, count)
                
                temp = [skyrows, max(skyrows)+2] ; pad the array
                ndx = [0, where(temp[1:count]-temp[0:count-1] ne 1L)+1]
                skyposition=temp[ndx[0:n_elements(ndx)-2]]
                skyaperture=temp[ndx[1:n_elements(ndx)-1]-1] - $
                  temp[ndx[0:n_elements(ndx)-2]]+1
                nsky = n_elements(skyposition)
                
             endif else print, format='(" At least one sky window MUST be selected!' + $
               string('',f='(A60)')+'",a1,$)', string(13b) 
             
          endelse 

       endwhile 
             
; overplot sky windows

       skyrows = skyrows[uniq(skyrows,sort(skyrows))] ; don't count rows twice
       nskyrows = n_elements(skyrows)

       skyimage = image[*,skyrows]
       skyinvvar = invvar[*,skyrows]

       first = 1L
       promptme = 1L
       changeorder:
       if not first then begin

          if promptme then splog, 'Select '+text+' to go back, '+text2+' to quit, or ENTER to continue.'+$
            string('',format='(A50)')
          print, format='(" ISKYSELECT: Specify NORDER_SKY = '+string(norder_sky,format='(i0)')+$
            ' to refit.'+'",A1,$)', string(13b)

          key = read_key(1) - 48
          if key eq 50 or key eq 18 then goto, select
          if key eq -38 then goto, done 
          if key eq 65 or key eq 33 then goto, quit
          if key ge 1 and key le 9 then norder_sky = key else goto, changeorder

          promptme = 0L

       endif 
                
       vals = total(image[(ncols-navg_sky)/2.:(ncols+navg_sky)/2.,*], 1) / navg_sky
       djs_iterstat, vals[skyrows], sigrej=4, mask=msk
       ok = where(msk)
       poly_iter, skyrows[ok], vals[skyrows[ok]], norder_sky-1, 3.0, polyfit, coeff=coeff
       if n_elements(coeff) gt 1L then polyfit = poly(rowaxis, coeff) else $
         polyfit = rowaxis*0+coeff[0]
       
       erase

       plot, rowaxis, /nodata, yrange=[0,1], ystyle=5, xstyle=7, thick=2, xthick=2, ythick=2, $
         position=[0.2,0.1,0.95,0.925], /normal
       for k=0,n_elements(skyposition)-1 do begin
          x = [skyposition[k], skyposition[k] + skyaperture[k]]
          polyfill, [x[0], x[1], x[1], x[0], x[0]], $
            [0, 0, 1, 1, 0], color=djs_icolor('dark gray')             
       endfor
       
       yrange1 = [min(vals)-0.10*(max(vals)-min(vals)), $
         max(vals)+0.15*(max(vals)-min(vals))]
       yrange2 = [min(vals)-0.1*(max(vals[skyrows[ok]])-min(vals[skyrows[ok]])), $
         max(vals[skyrows[ok]])+0.1*(max(vals[skyrows[ok]])-min(vals[skyrows[ok]]))]

       plot, rowaxis, vals, /noerase, psym=10, thick=2, xthick=2, ythick=2, $
         yrange=yrange1, ystyle=1, xstyle=3, charsize=1.5, charthick=2, $
         title='average of columns '+strcompress(fix((ncols-navg_sky)/2.),/remove_all)+ $
         '-'+strcompress(fix((ncols+navg_sky)/2.),/remove_all), $
         pos=[.2,.35,.95,.925], /normal, xtickname=replicate(' ', 30)
       legend, /top, /left, box=0, file+ ' ('+skyfileobject[i]+')', charsize=1.5, charthick=2
       djs_oplot, rowaxis, polyfit, color='black', thick=4
       djs_oplot, rowaxis, polyfit, color='green', thick=2

       plot, rowaxis, vals, /noerase, psym=10, thick=2, xthick=2, ythick=2, $
         yrange=yrange2, ystyle=1, xstyle=3, charsize=1.5, charthick=2, $
         xtitle='row (pixels)', position=[.2,.1,.95,.35], /normal
       xyouts, 0.05, 0.5, align=.5, orientation=90, charsize=1.5, charthick=2, 'counts', /normal
       djs_oplot, rowaxis, polyfit, color='black', thick=4
       djs_oplot, rowaxis, polyfit, color='green', thick=2
    
       first = 0L
       goto, changeorder
       
       done:

; write out the results

       nsky = n_elements(skyposition)
       printf, lun, f='(x,A'+len1+',2x,A'+len2+','+ $
         strcompress(nsky*2,/remove_all)+'f9.3,2x,i1)', $
         skylist[i], skyfileobject[i], skyposition*forage[i].cd2_2, $
         skyaperture*forage[i].cd2_2, norder_sky

; clean up memory

       icleanup, cube
       
    endfor 
    
    quit:

    free_lun, lun
    
    print

return
end
