PRO fg_extract, index, data, index_out, data_out,      $
                x0=x0, y0=y0, x1=x1, y1=y1,            $
                nx1=nx1, ny1=ny1,                      $
                loud=loud, verbose=verbose,            $
                version=version, name=name

;+
; NAME:
;	FG_EXTRACT
;
; PURPOSE:
;      Extract a sub-image from each image in the DATA array.
;
; CATEGORY:
;	FITS processing
;
; SAMPLE CALLING SEQUENCE:dex_out.crval1 = index_out.ycen
;	fg_extract, index, data, index_out, data_out,                       $
;                 [x0=x0], [y0=y0], [x1=x1], [y1=y1], [nx1=nx1], [ny1=ny1], $
;                 [/quiet],                                                 $
;                 [version=version, name=name ]
;
; INPUTS:
;	INDEX	  - the SSW index structure for each image
;	DATA	  - a 2- or 3-D data array (NX x NY [x N])
;
; OUTPUTS:
;       INDEX_OUT - the SSW header index structure for the corrected
;                   data. Currently this is just INDEX, unmodified.
;       DATA_OUT - The extracted image array (NX1 x NY1 [x N]).
;
; OPTIONAL INPUT KEYWORD PARAMETERS: 
;       Any two from the set {X0, X1, NX1}
;       Any two from the set {Y0, Y1, NY1}
;       LOUD       - Not set to minimize the messages
;	VERBOSE    - Set for a few more messages; overrides /loud
;
; OUTPUT KEYWORD PARAMETERS: 
;	VERSION    - The current version number of the program
;       NAME       - The name of this routine

;
; PROCEDURE:
;	1. Determine if sub-image parameters are acceptable.
;       2. Extract sub-images from data array.
;       3. Return data_out array of sub-images.
;
; MODIFICATION HISTORY:
;V0.1  Created 26-Apr-2007. T.E. Berger, LMSAL.
;V1.0  Heavily modified, added array capability. 10-Mar-2009, TEB, LMSAL.
;-
;-----------------------------------------------------------------------------
;
prognam = 'FG_EXTRACT.PRO'
progver = '1.0'

if KEYWORD_SET(loud) then loud = 1 else loud = 0
;Constants
xccd0 = 0 
yccd0 = 0
xccdmax = 4095
yccdmax = 2047

;Original image information
sz = SIZE(data)
if sz[0] lt 2 or sz[0] gt 3 then begin
   MESSAGE,'Image data must be a 2-D image or a 3-D array of 2-D images. Returning...'
   RETURN
end
nx0 = sz[1]
ny0 = sz[2]
if sz[0] eq 3 then nim = sz[3] else nim = 1
if loud then MESSAGE,/INFO,'Original image size = '+STRTRIM(nx0,2)+' x '+STRTRIM(ny0,2)+' pixels'

xlo = 0
ylo = 0
xhi = 0
yhi = 0
nx = 0
ny = 0
;Keyword check
if N_ELEMENTS(x0) gt 0 then begin
   xlo = 1
   if (x0 lt xccd0) || (x0 gt xccdmax) then begin
      MESSAGE,'Subimage X0 parameter < 0 or greater than CCD size. Returning...'
      RETURN
   end
end else x0=0

if N_ELEMENTS(y0) gt 0 then begin
   ylo = 1
   if (y0 lt yccd0) || (y0 gt yccdmax) then begin
      MESSAGE,'Subimage Y0 parameter < 0 or greater than CCD size. Returning...'
      RETURN
   end
end else y0=0

if KEYWORD_SET(x1) then begin
   xhi = 1
   if (x1 lt xccd0) || (x1 gt xccdmax) then begin
      MESSAGE,'Subimage X1 parameter < 0 or greater than CCD size. Returning...'
      RETURN
   end
end else x1 = nx0-1

if KEYWORD_SET(y1) then begin
   yhi = 1
   if (y1 lt yccd0) || (y1 gt yccdmax) then begin
      MESSAGE,'Subimage Y1 parameter < 0 or greater than CCD size. Returning...'
      RETURN
   end
end else y1 = ny0-1

if KEYWORD_SET(nx1) then begin
   nx = 1
   if (nx1 gt xccdmax) or (nx1 gt nx0) then begin
      MESSAGE,'Subimage size NX1 > original image or CCD size. Returning...'
      RETURN
   end
   if xhi and xlo then begin
      dx = x1 - x0 + 1
      if dx ne nx1 then begin
         MESSAGE,'Subimage size NX1 incompatible with X0 and X1 specification. Returning...'
         RETURN
      end
   end
end

if KEYWORD_SET(ny1) then begin
   ny = 1
   if (ny1 gt yccdmax) or (ny1 gt ny0) then begin
      MESSAGE,'Subimage size NY1 > original image or CCD size. Returning...'
      RETURN
   end
   if yhi and ylo then begin
      dy = y1 - y0 + 1
      if dy ne ny1 then begin
         MESSAGE,'Subimage size NY1 incompatible with Y0 and Y1 specification. Returning...'
         RETURN
      end
   end
end

;To extract an image you need one of the following combinations for
;each axis:
;
;    x0 & x1
;    x0 & nx
;    x1 & nx

if xlo && xhi then nx1 = x1 - x0 + 1 else $
if xlo && nx then x1 = x0 + nx1 - 1 else $
if xhi && nx then x0 = x1 - nx1 + 1 else $
   begin
   MESSAGE,/INFO,'Insufficient number of keywords supplied. Returning original image...'
   data_out = data
   RETURN
end

if ylo && yhi then ny1 = y1 - y0 + 1 else $
if ylo && ny then y1 = y0 + ny1 - 1 else $
if yhi && ny then y0 = y1 - ny1 + 1 else $
   begin
   MESSAGE,/INFO,'Insufficient number of keywords supplied. Returning original image...'
   data_out = data
   RETURN
end

case datatype(data) of
   'BYT': data_out = BYTARR(nx1,ny1,nim)
   'INT': data_out = INTARR(nx1,ny1,nim)
   'FLO': data_out = FLTARR(nx1,ny1,nim)
   'DOU': data_out = DBLARR(nx1,ny1,nim)
   else: begin
      MESSAGE,/INFO,'Only BYTE, INT, FLOAT, or DOUBLE arrays allowed. Returning original data'
      data_out = data
      RETURN
   end      
end 

index_out = index
for i=0,nim-1 do begin

   data_out[*,*,i] = data[x0:x1,y0:y1,i]

;Modify keywords
;Should we modify FGCCDIX0, etc.? I don't think so. Keeps original CCD readout reference in header.
   index_out[i].naxis1 = nx1
   index_out[i].naxis2 = ny1
   dxc = (x0 + x1)/2. - (nx0/2.)
   dyc = (y0 + y1)/2. - (ny0/2.)
   index_out[i].crpix1 = FLOAT(nx1)/2.0
   index_out[i].crpix2 = FLOAT(ny1)/2.0
   index_out[i].xcen = index[i].xcen + dxc*index[i].cdelt1
   index_out[i].crval1 = index_out[i].xcen
   index_out[i].ycen = index[i].ycen + dyc*index[i].cdelt2
   index_out[i].crval2 = index_out[i].ycen
   
end


RETURN
END