;+
; NAME:
;	READPHARO
;
; PURPOSE:
;	Read PHARO image file into a 1024x1024 array.
;
; EXPLANATION:
;	PHARO files are saved as either 512x512x4 (difference) or 512x512x8
;	(endpoint read) cubes.  In the case of endpoint read cubes, will
;	subtract first read from second, returning the difference.  This 
;	program is intended to replace READFITS.pro for PHARO images.
;
; CALLING SEQUENCE:
;	Result = READPHARO(Filename, Header, exten_no=, /silent)
;
; INPUTS:
;       Filename = Scalar string containing the name of the FITS file
;                 (including extension) to be read.   If the filename has
;                  a *.gz extension, it will be treated as a gzip compressed
;                  file.   If it has a .Z extension, it will be treated as a
;                  Unix compressed file.
;
; OUTPUTS:
;       Result = FITS data array constructed from designated record.
;                If the specified file was not found, then Result = -1
;
; OPTIONAL OUTPUT:
;       Header = String array containing the header from the FITS file.
;
; OPTIONAL INPUT KEYWORDS:
;       EXTEN_NO - non-negative scalar integer specifying the FITS extension to
;               read.  For example, specify EXTEN = 1 or /EXTEN to read the
;               first FITS extension.
;
;       /SILENT - Normally, READFITS will display the size the array at the
;                 terminal.  The SILENT keyword will suppress this
;
; PROCEDURES USED:
;	Functions: READFITS()
;	Procedures: SXADDPAR, SXDELPAR
;
; MODIFICATION HISTORY:
;	Create 05/07 by A. Bouchez, Caltech Optical Observatories
;-
FUNCTION        READPHARO, filename, header, exten_no=exten_no, silent=silent

  dat = READFITS(filename, header, exten_no=exten_no, silent=silent)

  sz = SIZE(dat)
  if sz[0] eq 3 then begin

    if not KEYWORD_SET(silent) then $
      MESSAGE, /info, 'Converting to ' + $
        STRING(sz[1:2]*2,f='(I0.0," by ",I0.0)')

    case sz[4] of
      2: im = INTARR(sz[1]*2, sz[2]*2)
      3: im = LONARR(sz[1]*2, sz[2]*2)
      4: im = FLTARR(sz[1]*2, sz[2]*2)
      5: im = DBLARR(sz[1]*2, sz[2]*2)
    endcase

    im[ sz[1]:sz[1]*2-1,     0:sz[2]-1  ] = dat[*,*,0]
    im[     0:sz[1]-1  ,     0:sz[2]-1  ] = dat[*,*,1]
    im[     0:sz[1]-1  , sz[2]:sz[2]*2-1] = dat[*,*,2]
    im[ sz[1]:sz[1]*2-1, sz[2]:sz[2]*2-1] = dat[*,*,3]

    if sz[3] eq 8 then begin
      unbias = im
      unbias[sz[1]:sz[1]*2-1,    0:sz[2]-1  ] = dat[*,*,4]
      unbias[    0:sz[1]-1  ,    0:sz[2]-1  ] = dat[*,*,5]
      unbias[    0:sz[1]-1  ,sz[2]:sz[2]*2-1] = dat[*,*,6]
      unbias[sz[1]:sz[1]*2-1,sz[2]:sz[2]*2-1] = dat[*,*,7]
      im = unbias - im
    endif

    if N_PARAMS() gt 1 then begin
      SXADDPAR, header, 'NAXIS', 2
      SXADDPAR, header, 'NAXIS1', sz[1]*2
      SXADDPAR, header, 'NAXIS2', sz[2]*2
      SXDELPAR, header, 'NAXIS3'
    endif
  endif else begin

    MESSAGE, /info, 'Dimensions do not match those of a PHARO image cube!'
    im = -1
  endelse

RETURN,im
END