;+
; NAME:
;  moscal
; PURPOSE:
;  Apply standard CCD image correction steps to a raw group-FITS image.
; DESCRIPTION:
;
; CATEGORY:
;  CCD data processing
;
; CALLING SEQUENCE:
;  moscal,root,outdir,start,nframes,bias,dark,flat
;
; INPUTS:
;  root    - Root of the file name (you must include . in the name).
;  outdir  - The directory to write the calibrated files to.  Don't make
;              this the same as the input directory.
;  start   - First frame number to read (integer or long).
;               Start can also be a vector of explicit frame numbers to load.
;               In this case, nframes need not be specified and in fact will
;               be ignored.
;  nframes - Number of frames to average.
;  bias    - Bias frame image name to subtract from each raw frame.
;  dark    - Dark frame image name to subtract from each raw frame.
;  flat    - Flat frame image name to subtract from each raw frame.
;
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORD INPUT PARAMETERS:
;
;   CROP     = region of original image to save, default=no cropping.
;                 [x1,x2,y1,y2]
;
;   EXCLUDE - Optional vector of image numbers that should be excluded from
;                average.  Default is to include all frames.
;
;   EXPKEY = String - FITS keyword to read to get exposure time, default = EXPTIME
;
;   OVERSCAN = column overscan region to use for frame bias level,
;                 default=no overscan subtraction.
;
; OUTPUTS:
;   The calibrated images are written to outdir.  Don't make this the same
;   as the current directory!
;
; KEYWORD OUTPUT PARAMETERS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;  2000/02/03, Written by Marc W. Buie, Lowell Observatory
;  2000/02/28, MWB, added support for frame numbers > 999.
;-
pro moscal,root,outdir,start,nframes,bias,dark,flat, $
   EXCLUDE=exclude,SCALE=in_scale,EXPKEY=expkey, $
   OVERSCAN=in_overscan,CROP=in_crop,GAINCOR=in_gaincor

   if badpar(root,   7,    0,caller='MOSCAL: (root) '   ) then return
   if badpar(outdir, 7,    0,caller='MOSCAL: (outsuf) ' ) then return
   if badpar(start,  [2,3],[0,1],caller='MOSCAL: (start) ',rank=start_rank) then return

   if badpar(exclude,[0,2,3],[0,1],caller='MOSCAL: (exclude) ',default=-1) then return

   if badpar(in_overscan,[0,2,3],[1,2],caller='MOSCAL: (overscan) ',rank=o_rank) then return
   if badpar(in_crop,    [0,2,3],[1,2],caller='MOSCAL: (crop) ',rank=c_rank) then return
   if badpar(expkey,     [0,7],      0,caller='MOSCAL: (EXPKEY) ',default='EXPTIME') then return
   if badpar(in_gaincor, [0,4,5],  [1],caller='MOSCAL: (GAINCOR) ',default=1.0) then return

   if badpar(bias,   [0,7],0,caller='MOSCAL: (bias) ',default='[[none]]') then return
   if badpar(dark,   [0,7],0,caller='MOSCAL: (dark) ',default='[[none]]') then return
   if badpar(flat,   [0,7],0,caller='MOSCAL: (flat) ',default='[[none]]') then return

   outdir = addslash(outdir)

   ; Check to see if it's a sequential list or random list.
   if start_rank eq 0 then begin
      frames=start+indgen(nframes)
      if badpar(nframes,[2,3],   0,caller='MKBIAS: (nframes) ') then return
   endif else begin
      frames=start
      nframes=n_elements(frames)
   endelse

   ; Setup the file name format string
   digits = fix(ceil(alog10(max(frames)+1))) > 3
   dig = strn(digits)
   fnfmt = '(i'+dig+'.'+dig+')'

   ; Apply the exclusion criteria to the frames list.  Then, find out how
   ;   many frames are being collected and make sure there is something to do.
   for i=0,nframes-1 do begin
      z=where(frames[i] eq exclude,count)
      if count ne 0 then frames[i]=-1
   endfor
   zg=where(frames ne -1,countg)
   if countg eq 0 then $
      message,'Error ** you have excluded all frames, nothing to do.'

   ; Make the first file name
   fname=root+string(frames[0],format=fnfmt)
   if exists(fname+'.fits') then fname=fname+'.fits'

   ; Check header of image to see if it is a multi-extension image.
   hdr=headfits(fname)
   numext=sxpar(hdr,'NEXTEND')

   if numext le 1 then begin
      print,'MOSCAL: This program only works on group fits files.'
      return
   endif

   if n_elements(in_gaincor) eq 1 then begin
      gaincor = replicate(in_gaincor,numext)
   endif else if n_elements(in_gaincor) ne numext then begin
      print,'MOSCAL: Gain correction vector has the wrong number of elements.'
      return
   endif else begin
      gaincor = in_gaincor
   endelse

   ; Setup the overscan/crop control values
   if numext eq 0 then begin
      extend=0
      numext=1
      if o_rank eq 0 then begin
         do_overscan=0
      endif else begin
         do_overscan=1
         overscan = in_overscan
      endelse
      if c_rank eq 0 then begin
         do_crop=0
      endif else begin
         do_crop=1
         crop = in_crop
      endelse
   endif else begin
      extend=1
      if o_rank eq 0 then begin
         do_overscan=0
      endif else if o_rank eq 1 then begin
         do_overscan=1
         overscan = rebin(in_overscan,n_elements(overscan),numext)
      endif else begin
         do_overscan=1
         overscan = in_overscan
      endelse
      if c_rank eq 0 then begin
         do_crop=0
      endif else if o_rank eq 1 then begin
         do_crop=1
         crop = rebin(in_crop,n_elements(crop),numext)
      endif else begin
         do_crop=1
         crop = in_crop
      endelse
   endelse

   ; Main loop over extension, done just once on "normal" frames
   for ix=1,numext do begin
      ix0=ix-1

      print,'Processing extension ',strn(ix)

      ; Read the bias frame for this image extension
      if bias ne '[[none]]' then begin
         biasim=0
         biasim=readfits(bias,exten_no=ix,/silent)
      endif

      ; Read the dark frame for this image extension
      if dark ne '[[none]]' then begin
         darkim=0
         darkim=readfits(dark,exten_no=ix,/silent)
      endif

      ; Read the flat frame for this image extension
      if flat ne '[[none]]' then begin
         flatim=0
         flatim=readfits(flat,exten_no=ix,/silent)
      endif

      ; Loop over the frames
      for i=0,nframes-1 do begin

         fname=root+string(frames[zg[i]],format=fnfmt)
         if exists(fname+'.fits') then fname=fname+'.fits'

         ; silently skip the frame if it doesn't exist.
         if exists(fname) then begin

            ; read in the new image, the dummy line before readfits is to dump
            ;   the storage for the array before its used again.
            image = 0
            if ix eq 1 then begin
               hdr=headfits(fname)
               writefits,outdir+fname,0,hdr
            endif

            image = readfits(fname,hdr,exten_no=ix,/silent)

            if do_overscan and do_crop then begin
               image = colbias(image,overscan[0,ix0],overscan[1,ix0], $
                        crop[0,ix0],crop[1,ix0],crop[2,ix0],crop[3,ix0],biasval=biasval)
               if extend then $
                  print,fname,' overscan value is ',biasval,' extension ',strn(ix) $
               else $
                  print,fname,' overscan value is ',biasval

            endif else if do_overscan then begin
               image = colbias(image,overscan[0,ix0],overscan[1,ix0],biasval=biasval)
               if extend then $
                  print,fname,' overscan value is ',biasval,' extension ',strn(ix) $
               else $
                  print,fname,' overscan value is ',biasval

            endif else if do_crop then begin
               image = image[crop[0,ix0]:crop[1,ix0],crop[2,ix0]:crop[3,ix0]]
               if extend then $
                  print,fname,' extension ',strn(ix) $
               else $
                  print,fname

            endif else begin
               if extend then $
                  print,fname,' extension ',strn(ix) $
               else $
                  print,fname

            endelse

            if bias ne '[[none]]' then $
               image = temporary(image) - biasim

            if dark ne '[[none]]' then begin
               exptime=float(sxpar(hdr,expkey))
               image = temporary(image) - darkim*exptime
            endif

            if flat ne '[[none]]' then $
               image = temporary(image)/flatim

skysclim,image,lowval,hival,meanval1,sigma1,npts=4000
            if gaincor[ix-1] ne 1.0 then $
               image = temporary(image)*gaincor[ix-1]
skysclim,image,lowval,hival,meanval2,sigma2,npts=4000
print,ix,meanval1,gaincor[ix-1],meanval2

            image = fix(temporary(image < 65535.0) - 32767.5)
            sz=size(image,/dimen)
            sxaddpar,hdr,'NAXIS1',sz[0]
            sxaddpar,hdr,'NAXIS2',sz[1]
            sxaddpar,hdr,'BSCALE',1.0,AFTER='NAXIS2'
            sxaddpar,hdr,'BZERO',32768.0,AFTER='BSCALE'
            sxaddpar,hdr,'BITPIX',16
            writefits,outdir+fname,image,hdr,/append

         endif ; frame processing clause

      endfor ; frame loop

   endfor

end