read_cube_fits.pro

; FUNCTION READ_CUBE_FITS
;
; Reads in a fits file generated by cube, and places the header in a
; structure. Scales the fits file, undoes the log from cube in tipsy.
;
; Inputs:
;
;    INFILE: The fits data cube to be read.
;
;    HEADER: A name for the header structure created.
;
;    MRDH: Array of strings describing the header as read by MRDFITS.
;
;
; Return value: The data cube stored as a 3d array
;
;
function read_cube_fits, infile, header, mrdh, kpcunit=kpcunit


if not keyword_set(infile) then begin
    print,"Syntax:"
    print,""
    print,"read_cube_fits(infile, [header, mrdh])"
    print,""
    return,0.
endif

array = mrdfits(infile, 0, mrdh, /silent)
;array = double(10^(double(array)*1.0))
header = {cubeheader, $
          naxis: 0, $
          naxis1: 0, $
          naxis2: 0, $
          naxis3: 0, $
          crval1: 0.0, $
          cdelt1: 0.0, $
          crpix1: 0.0, $
          crval2: 0.0, $
          cdelt2: 0.0, $
          crpix2: 0.0, $
          crval3: 0.0, $
          cdelt3: 0.0, $
          crpix3: 0.0, $
          bscale: 0.0, $
          bzero: 0.0, $
          blank: 0.0, $
          kpcunit: 0.0, $
          munit: 0.0 $
         }

header.naxis = sxpar(mrdh, 'NAXIS')
header.naxis1 = sxpar(mrdh, 'NAXIS1')
header.naxis2 = sxpar(mrdh, 'NAXIS2')
header.naxis3 = sxpar(mrdh, 'NAXIS3')

header.crval1 = sxpar(mrdh, 'CRVAL1')
header.cdelt1 = sxpar(mrdh, 'CDELT1')
header.crpix1 = sxpar(mrdh, 'CRPIX1')

header.crval2 = sxpar(mrdh, 'CRVAL2')
header.cdelt2 = sxpar(mrdh, 'CDELT2')
header.crpix2 = sxpar(mrdh, 'CRPIX2')

header.crval3 = sxpar(mrdh, 'CRVAL3')
header.cdelt3 = sxpar(mrdh, 'CDELT3')
header.crpix3 = sxpar(mrdh, 'CRPIX3')

header.bscale = sxpar(mrdh, 'BSCALE')
header.bzero = sxpar(mrdh, 'BZERO')
header.blank = sxpar(mrdh, 'BLANK')

; scale the array
array = array * header.bscale + header.bzero
array = 10.0^(float(array))

; set blanks to 0.
vals = where(array eq 10^(float(header.blank * header.bscale + header.bzero)))
if (vals[0] ne -1) then array[vals] = 0.0

if keyword_set(kpcunit) then begin
    header.crval1 = header.crval1 * kpcunit
    header.crval2 = header.crval2 * kpcunit
    header.cdelt1 = header.cdelt1 * kpcunit
    header.cdelt2 = header.cdelt2 * kpcunit
endif

return,array

end