;pro Fits_Read_DataCubeRegion
;
;PURPOSE:
;	To return a 3-d region defined by xstart - xstop, ystart - ystop, for
;	the time interval zstart - zstop
;
;	The read will be handled in three different ways:
;
;	1 - Full Frame  - No loop needed to read in all data since
;			  it is contiguous
;	2 - Sub-Frame   - One loop needed since all bytes are
;			  contiguous until a new frame is requested
;	3 - Sub-Region  - Two loops are needed since the bytes for
;			  the rows, columns and frames are all not contiguous
;
;INPUTS:
;	XStart - The beginning column pixel number
;	XStop  - The end column pixel number
;	YStart - The beginning row pixel number
;	YStop  - The end row pixel number
;	ZStart - The beginning frame number
;       WCoord - The entry in naxis4 that is desired (default 0 for datacube)
;                (only supported for full NAXIS1xNAXIS2 frames)
;	ZStop  - The end frame number
; 	Naxis1 - The Number of Columns in the full frame
;	Naxis2 - The Number of Rows in the full frame
;
;COMMON:
;	FCB - the file control block of the (already opened fits file)
;	      FCB.bitpix = Number of Bits per pixel
;	      FCB.bzero  = Zero point of Bytes
;	      FCB.start_data = the position in bytes of the start of data
;
;	NOTE: The ieee functions are used to make sure the bytes are in
;	      the correct order
;
;KEYWORDS:
;       OFFSET - An offset in data words corresponding to the location
;                where data actually starts.  This value is intended
;                to correct files that were incorrectly written with
;                a 1440 or 2880 word offset between where the
;                header ends and where the data begins.
;
;
;OUTPUTS:
;	CubeRegion :
;	      An (XStop - XStart, YStop - YStart, ZStop - ZStart)
;	      array holding the data
;
function Fits_Read_DataCube_Region, FitsFileName, FCB, FitsHeader,$
	      XStart, XStop, YStart, YStop, ZStart, ZStop, Offset = Offset, $
              WCoord=WCoord

If N_Elements(Offset) eq 0 then Offset = 0
If N_Elements(WCoord) eq 0 then WCoord = 0

;Test to see if file is already open.  If not, open it with fits_open_datacube
If (Size(Fcb))[0] eq 0 then $
   Pos = Fits_Open_DataCube(FitsFileName, FCB, FitsHeader, 'Read')

;Make sure that gigantic cubes can be handled by making the pointer a double
XStartL = double(XStart) & XStopL = double(XStop)
YStartL = double(YStart) & YStopL = double(YStop)
ZStartL = double(ZStart) & ZStopL = double(ZStop)
OffsetL = double(Offset)

;Get the full frame dimensions
Naxis1 = fcb.axis(0)
Naxis2 = fcb.axis(1)
Naxis3 = fcb.axis(2)
Naxis4 = fcb.axis(3)

;Form a 3-d Array based on the Input Coordinates
XSize = XStopL - XStartL + 1
YSize = YStopL - YStartL + 1
ZSize = ZStopL - ZStartL + 1

;Get the number of bits per pixel
BitPix = fcb.bitpix(0)
Bzero = fcb.bzero(0)

;Check to see if we need to swap bytes
Bswap = 1 - Is_Ieee_Big()

case BitPix of
   8 : begin
     IDL_type = 1
     DataCubeRegion = BytArr(XSize, YSize, ZSize)
   end
   16: begin
     if BZero ne 0 then begin
       IDL_type = 12
       DataCubeRegion = UIntArr(XSize, YSize, ZSize)
     endif else begin
       IDL_type = 2
       DataCubeRegion = IntArr(XSize, YSize, ZSize)
     endelse
   end
   32: begin
     IDL_type = 3
     DataCubeRegion = LonArr(XSize, YSize, ZSize)
   end
  -32: begin
     IDL_type = 4
     DataCubeRegion = FltArr(XSize,YSize,ZSize)
   end
  -64: begin
     IDL_type = 5
     DataCubeRegion = DblArr(XSize,YSize,ZSize)
   end
endcase

BytesPerWord = (abs(bitpix)/8)
Position = fcb.start_data
ArrSize  = Naxis1*double(Naxis2)
CubeSize = double(Naxis3)*ArrSize
 
;Determine wheter we can get read in an entire array without skipping any
;bytes, i.e., an entire contiguous frame

;CASE 3 --------------------------------------------------------Sub Region
If not(Naxis1 eq XSize) then begin

  DataRow = make_array(dim = [XSize], type = idl_type, /nozero)

  ;Loop through the datacube and obtain the pixel values
  For Frame = ZStartL, ZStopL do begin
    for Row = YStartL, YStopL do begin

	;Move to the right place in the Data
	point_lun, FCB.unit, Position + BytesPerWord * ( $
			(Frame*ArrSize)+$   	; Distance Along Z
			(Row  * Naxis1) +$   	; Distance Along Y
			(XStartL)+$	    	; Distance Along X
			OffsetL)
	;Read in one row and put it into the datacube that will be returned
	readu, FCB.unit, DataRow

	if bswap then ieee_to_host, DataRow
  	If bitpix eq 16 then begin
	  DataCubeRegion(*, Row - YStartL, Frame - ZStartL) = $
		       UInt(DataRow) + UInt(BZero)
	Endif else begin
          DataCubeRegion(*, Row - YStartL, Frame - ZStartL) = $
                      DataRow + bzero
	EndElse
     endfor
  Endfor

;CASE 2 ---------------------------------------------------------Sub Frame
Endif Else If (Naxis1 eq XSize) and not(Naxis2 eq YSize) then Begin

  DataSubFrame = make_array(dim = [XSize*YSize], type = idl_type, /nozero)

  For Frame = ZStartL, ZStopL do begin
     ;Move to the right place in the Data
     point_lun, FCB.unit, Position + BytesPerWord * ( $
                      (Frame  * ArrSize)+$       ; Distance Along Z
                      (YStartL * Naxis1)+$      ; Distance Along Y
 		      OffsetL)
     ;Read in one row and put it into the datacube that will be returned
     readu, FCB.unit, DataSubFrame
     if bswap then ieee_to_host, DataSubFrame

     If bitpix eq 16 then begin
        DataCubeRegion(*, *, Frame - ZStartL) = $
                       UInt(DataSubFrame) + UInt(Bzero)
        Endif else begin
        DataCubeRegion(*, *, Frame - ZStartL) = $
                      DataSubFrame + bzero
     EndElse

   Endfor


;CASE 1 ---------------------------------------------------------Full Frame
Endif Else Begin
  ;Full Frame should be much faster if we don't go row by row and move the
  ;pointer in the readu function

  DataFrame = make_array(dim = [ArrSize], type = idl_type, /nozero)
  For Z = ZStartL, ZStopL do begin
    ;Move to the right place in the Data
    point_lun, FCB.unit, Position + BytesPerWord*(WCoord*CubeSize+Z*ArrSize+OffsetL)
    readu, FCB.unit, DataFrame
    if bswap then ieee_to_host, DataFrame

    If bitpix eq 16 then begin
      DataCubeRegion(*,*,Z-ZStartL) = $
	reform(UInt(DataFrame), XSize, YSize) + UInt(Bzero)
    EndIf else begin
      DataCubeRegion(*,*,Z-ZStartL) = $
	reform(DataFrame, XSize, YSize) + bzero
    EndElse
  EndFor

EndElse

return, DataCubeRegion

end