;+
; NAME:
;	Radio Astronomy Group Fits Read  
; PURPOSE:
;	Reads RAG spectrograms from FITS files 
;	including non-regular axes. 
;	The procedure reads also with gzip compressed files
;	  if they end with the extention ".gz"
; CALLING SEQUENCE:
;	RAGFitsRead, filename[ , image, xAxis, yAxis ]
; INPUT:
;	filename: string containing name of the file with extension
; OUTPUTS:
;	image: 2D array containing the array read.
;	xAxis, yAxis: 1D arrays containing the axes.
; KEYWORDS:
;	used for getting header informations. See "RAGFitsWrite"
;	for default information written.
;	ORIGIN, TELESCOPE, INSTRUMENT,
;	OBJECT: name of observed object
;	DATEOBS: observation date in form dd/mm/yy
;	TIMEOBS:	start time of observation ( hh:mm:ss.ddd )
;		or "No Time"
;	DATEEND, TIMEEND: as DATEOBS and TIMEOBS,
;		but for the end of the observation.
;	BZERO, BSCALE: If they are not 0 and 1, respectively, 
;		the data is reconstructed with the formula
;		with image = BZERO + BSCALE*data
;	BUNIT: string containing the units of the image pixels.
;	CTYPE1, CTYPE2: string containing the FITS type of
;		physical coordinates.
;	COMMENT, HISTORY: strings for description
;		of the image
;	NOSCALE: if set, the reconstruction with BZERO and
;		BSCALE is not done.
;	SILENT: if set, the FITS header is not printed on screen
;	DATAMIN, DATAMAX: Max and min of the data.
;	RELATIVETIME: If set, the decimal value of TIMEOBS
;		 is NOT added to the axis values. 
; SIDE EFFECT:
;	A fits file is opened.
; MODIFICATION HISTORY
;	Created: A. Csillaghy, ETHZ, October 1992
;	For IDL Sun in May 1993, A.Cs.
;	SILENT in june, 93, A.Cs
;	DATAMIN, DATAMAX; When only filename provided, the header
;		is read but not the array. Oct. 93, A.Cs
;	RELATIVETIME in august 1995, ACs
;	Read also compressed file in November 95, ACs
;       Adaptation for IDL5/SSW/Ragview in March 98 -- ACs
;-

FUNCTION RegularAxis, xtension, nEls, header 

    IF xtension EQ '1' THEN xtensionName = 'the time' $
    ELSE xtensionName = 'the frequency'

    crVal = SxPar( header, "CRVAL"+ xtension )
    crPix = SxPar( header, "CRPIX"+ xtension )
    cDelt = SxPar( header, "CDELT"+ xtension )
    Message, 'Axis corresponding to '+ xtensionName + $
	' is set regularly', /INFO, /CONT
    IF cDelt EQ 0 THEN BEGIN
      cDelt = 1
      Message, 'Warning in '+ xtensionName+$
	' : step was 0  => set to 1', /INFO, /CONT
    ENDIF
    RETURN,  LIndGen( nEls )*cDelt + crVal 

END ; Set Regular Axes


PRO SetAxes, filename, nx, ny,  header,  xAxis, yAxis, SILENT = silent

textstr = ''
FxBOpen, unit, filename, 1, binaryHeader, ERRMSG = textstr

IF textstr NE '' THEN BEGIN
    xAxis = RegularAxis( '1', nx, header )
    yAxis = RegularAxis( '2', ny, header )
ENDIF ELSE BEGIN
    IF NOT Keyword_Set( SILENT ) $
      AND header(0) NE '' THEN BEGIN
        Print, 'Header of the extension: '
        Print, binaryHeader
    ENDIF    
    FxBRead, unit, xAxis, 1, ERRMSG=textstr
    IF textstr NE '' THEN BEGIN
        Print, 'Ragfitsread: ' + textstr
        Print, 'Setting axis to pixel number'
        xAxis = FIndGen( nx )
    END
    FxBRead, unit, yAxis, 2, ERRMSG=textstr
    IF textstr NE '' THEN BEGIN
        Print, 'Ragfitsread: ' + textstr
        Print, 'Setting axis to pixel number'
        yAxis = FIndGen( ny )
    END
    
    axisSize = Size( xAxis )
    IF axisSize(0) GT 1 THEN BEGIN
        xAxis = xAxis( *, 0 )
        yAxis = yAxis(*, 1 )
    ENDIF
ENDELSE

FxBClose, unit, ERRMSG=textstr
IF textstr NE '' THEN BEGIN
    Print, 'Ragfitsread: ' + textstr
END

END                             ; set axis


PRO StoreKwrd, header, kwrdName, default,  kwrdVar

  kwrdVar = SxPar( header, kwrdName )
  IF !err LT 0 THEN BEGIN
    Message, 'Warning: the keyword '+ kwrdName + $
	' is not defined', /INFO
    kwrdVar = default
  ENDIF
  
END; store keyword

PRO RAGFitsRead, filenameP, image, xAxis, yAxis, $
 	ORIGIN = origin,  HEADER = header, $
	TELESCOPE = telescope, INSTRUMENT = instrument, $
	OBJECT= object, CONTENT = content, $
	DATEOBS = dateObs, TIMEOBS =  timeObs, $
 	DATEEND = dateEnd, TIMEEND =  timeEnd, $
	BZERO = bZero, BSCALE = bScale, BUNIT = bUnit, $
	CTYPE1 = cType1, CTYPE2 = cType2, $
	COMMENT = comment, SILENT = silent,  $
	HISTORY = history, NOSCALE = noScale, $
	DATAMIN = dataMin, DATAMAX = dataMax, $
	RELATIVETIME = relativeTime
	
  On_Error, 2

  IF N_Elements( filenameP ) EQ 0 THEN BEGIN
    Message, 'Usage: RagFitsRead, filename [ , image, xAxis, yAxis ]', $
	/INFO, /CONT 
    RETURN
  ENDIF

  fileName = fileNameP

  fileFound = FindFile( filename )
  fileCompressed = StrPos( filename, '.gz', StrLen(fileName)-3)
  isTmpFile = 0
  IF filefound(0) EQ "" OR fileCompressed(0) NE -1 THEN BEGIN
    Spawn, 'which gunzip', gunzip
    isGunzip = StrPos( gunzip, 'not found' )
    IF isGunzip(0) NE -1 THEN BEGIN
      Message, 'file is compressed but there is no gunzip command available', $
 	/INFO, /CONT 
      RETURN
    ENDIF
    IF fileCompressed(0) EQ -1 THEN BEGIN
      fileFound = FindFile( filename + '.gz' )
      IF fileFound(0) EQ "" THEN BEGIN
        Message, 'compressed file not found', /INFO, /CONT 
        RETURN
      ENDIF ELSE fileName =  fileName + '.gz'
    ENDIF
    Spawn, gunzip  + ' -c ' + filename + ' > ragfitsreadtmp.fit'
    fileName = 'ragfitsreadtmp.fit'
    isTmpFile = 1
  ENDIF
        

  IF NOT Keyword_Set( NOSCALE ) THEN noScale = 0
  IF N_Params() GT 1 THEN BEGIN
      errmsg = ''
      FxRead, filename, image, header, NOSCALE = noScale, ERRMSG = errmsg 
      IF errmsg NE '' THEN BEGIN
          Print, "Not read: " +  errmsg
          !err = 1
          RETURN
      ENDIF
  ENDIF ELSE BEGIN
    OpenR, unit, fileName, /GET_LUN
    FxHRead, unit, header
    Free_LUn, unit
  ENDELSE
  IF N_Elements( header ) EQ 0 THEN RETURN
  

  IF NOT Keyword_Set( SILENT ) THEN BEGIN
    Print, 'Header of FITS file : '
    Print, header
  ENDIF

  StoreKwrd, header, 'BITPIX', 8, bitpix
  StoreKwrd, header, 'CONTENT', '', content
  StoreKwrd, header, 'ORIGIN', '', origin
  StoreKwrd, header, 'TELESCOP', '', telescope
  StoreKwrd, header, 'INSTRUME', '', instrument
  StoreKwrd, header, 'OBJECT', '', object
  StoreKwrd, header, 'DATE-OBS', '',dateObs
  StoreKwrd, header, 'TIME-OBS', 'No Time', timeObs
  StoreKwrd, header, 'DATE-END', '', dateEnd
  StoreKwrd, header, 'TIME-END', '', timeEnd
  StoreKwrd, header, 'BZERO', 0, bZero
  StoreKwrd, header, 'BSCALE', 1, bScale
  StoreKwrd, header, 'BUNIT', '', bUnit
  StoreKwrd, header, 'CTYPE1', '', cType1
  StoreKwrd, header, 'CTYPE2', '', cType2
  StoreKwrd, header, 'COMMENT', '', comment
  StoreKwrd, header, 'HISTORY', '', history
  StoreKwrd, header, 'DATAMIN', '', dataMin
  StoreKwrd, header, 'DATAMAX', '', dataMax

  IF N_Params() GT 1 THEN BEGIN

      nx = N_Elements( image(*,0))
      ny = N_Elements( image(0,*))

      IF N_Params() GT 2 THEN BEGIN 
          SetAxes, filename, nx, ny, header, xAxis, yAxis, SILENT = silent
          IF ((timeObs NE 'No Time') OR $
              (StrCompress(timeObs) NE ' ')) AND $
            (NOT Keyword_Set(RELATIVETIME)) THEN $
            xAxis = xAxis + DecConvert( timeObs )
      ENDIF
      
      IF bitPix EQ -32   THEN  image = Float( image )
      
  ENDIF
  IF isTmpFile THEN Spawn, 'rm -f ragfitsreadtmp.fit'
  
  !err = 0 
  
END