;pro fits_read_datacube
;
;PURPOSE:
;       To read in a Fits datacube or a subregion of that datacube (in x,y, and z) and a header.
;
;INPUTS:
;       FitsFileName - The name of the file to be opened
;
;       FitsHeader   - The fits header.  It must contain the minimum amount of information necessary
;                      to write a fits file.
;
;KEYWORDS:
;	XSTART   - The first x coordinate of the subregion (0 - (Naxis1-1))
;	XSTOP    - The last x coordinate of the subregion (XStart - (Naxis1-1))
;       YSTART   - The first y coordinate of the subregion (0 - (Naxis2-1))
;       YSTOP    - The last y coordinate of the subregion (YStart - (Naxis2-1))
;       ZSTART   - The first z coordinate of the subregion (0 - (Naxis3-1))
;       ZSTOP    - The last z coordinate of the subregion (ZStart - (Naxis3-1))
;       WCOORD   - The entry in NAXIS4 that is desired (default 0 for datacube)
;                  (only supported for full NAXIS1xNAXIS2 frames)
;	MASK     - A 1-dimensional array of indices
;
;	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.
;
;       NOTE: The first and last keywords will yield bizarre results if the coordinates
;	      do not result in a cube
;
;	FIRST    - The first word of data to be returned
;	LAST     - The last word of data to be returned
;       BZEROSET: int
;         0 - Use BZero value from header (If not present, data will be signed integers)
;         1 - Force BZero = 32768 for unsigned integers
;
;RETURN:
;		Im - The subregion array
;
;		FitsHeader - The header of the .fits file
;
;CALLING SEQUENCE:
;	Fits_Read_DataCube, 'Test.fits', Image, Header
;
;AUTHOR:
;   Lance Simms, Stanford University 12/08
;****************************************************************

pro Fits_Read_Datacube, FitsFileName, Im, FitsHeader, Xstart = XStart, YStart= YStart, $
	ZStart = ZStart, XStop = XStop, YStop = YStop, ZStop = ZStop, First = First, $
	Last = Last, Offset = Offset, BZeroSet = BZeroSet, WCoord=WCoord

;Open the file and get some word from the header
if N_Elements(BZeroSet) eq 0 then BZeroSet = 0
If N_Elements(WCoord)   eq 0 then WCoord   = 0

Position = Fits_Open_DataCube(FitsFileName, FCB, FitsHeader, 'Read', BZeroSet = BZeroSet)

Naxis      = long(sxpar(FitsHeader, 'NAXIS'))
Naxis1     = long(sxpar(FitsHeader, 'NAXIS1'))
Naxis2     = long(sxpar(FitsHeader, 'NAXIS2'))
Naxis3     = long(sxpar(FitsHeader,'NAXIS3'))
Naxis4     = long(sxpar(FitsHeader,'NAXIS4'))

If Naxis3 eq 0 then Naxis3 = 1
ArrSize    = long(Naxis1*Naxis2)

If N_Elements(XStart) eq 0 then XStart = 0
If N_Elements(XStop)  eq 0 then XStop  = Naxis1-1
If N_Elements(YStart) eq 0 then YStart = 0
If N_Elements(YStop)  eq 0 then YStop  = Naxis2-1
If N_Elements(ZStart) eq 0 then ZStart = 0
If N_Elements(ZStop)  eq 0 then ZStop  = Naxis3-1
If N_Elements(Offset) eq 0 then Offset = 0

;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)

;Accomodate old FIRST and LAST keywords
If N_Elements(FIRST) ne 0 and N_Elements(LAST) ne 0 then begin
   First   = long(First)
   ZStartL = double(First / ArrSize)
   XStartL = double((First mod ArrSize) mod Naxis1)
   YStartL = double((First mod ArrSize) / Naxis1)

   Last    = Long(Last)
   ZStopL  = double(Last  / ArrSize)
   XStopL  = double((Last mod ArrSize) mod Naxis1)
   YStopL  = double((Last mod ArrSize)  / Naxis1)
EndIf

;Check Ranges
If (XStopL lt XStartL) or (XStartL lt 0) or (XStopL lt 0) or $
   (XStartL gt Naxis1-1) or (XStopL gt Naxis1-1) then begin
   print, 'X coordinates are out of range'
   stop
EndIf
If (YStopL lt YStartL) or (YStartL lt 0) or (YStopL lt 0) or $
   (YStartL gt Naxis2-1) or (YStopL gt Naxis2-1) then begin
   print, 'Y coordinates are out of range'
   stop
EndIf
If (ZStopL lt ZStartL) or (ZStartL lt 0) or (ZStopL lt 0) or $
   (ZStartL gt Naxis3-1) or (ZStopL gt Naxis3-1) then begin
   print, 'Z coordinates are out of range'
   stop
EndIf
If WCoord gt NAxis4 then begin
   print, 'W Coordinate is out of range'
   stop
endif

;Read in the subregion
Im = Fits_Read_DataCube_Region(FitFileName, FCB, FitsHeader, $
                double(XStartL), double(XStopL), $
		double(YStartL), double(YStopL), $
		double(ZStartL), double(ZStopL), $
		Offset = OffsetL, WCoord=WCoord)

;Close the File
FitsClose = Fits_Close_DataCube(FCB)

end