;+
; NAME:
;       EIS_DATA__READFITS
;
; PURPOSE:
;       EIS_DATA__READFITS reads EIS data files and stores data and other
;       information in the object fields of an EIS_DATA class object.
;
; CATEGORY:
;       Hansteen/Wikstol Data analysis SW
;
; CALLING SEQUENCE:
;       An object of class EIS_DATA is declared with this statement:
;                   data = obj_new('eis_data')
;
;       EIS_DATA__READFITS reads EIS data fits files and stores information
;       into the object, with this statement:
;                     data-> read, filename, [hdr = hdr]
;
;       If called from within an eis_fits_data procedure, the
;       EIS_DATA__READFITS statement looks like this:
;                     self-> read, file [, hdr = hdr]
;
; INPUTS:
;       file:  Name of fits file (string) containing EIS data (and headers).
;
; KEYWORD PARAMETERS:
;       hdr: Set this keyword to obtain also the header
;       information in the form of an object of type EIS_FITS_HDR.
;
; OUTPUTS:
;       Reads EIS data fits files and fills the parameters of the EIS_DATA
;       object with information.
;
; CALLS:
;       FITS_HDR__DEFINE
;       FITS_HDR__READ
;       FXBOPEN
;       FXBREAD
;       FXBTFORM
;       FXBTDIM
;
; COMMON BLOCKS:
;
;
; PROCEDURE:
;
;
;
; RESTRICTIONS:
;
;
; MODIFICATION HISTORY:
;       14-Jan-2003: Oivind Wikstol.
;       20-Apr-2004: Oivind Wikstol - Cleaned up
;       17-Jun-2004: Oivind Wikstol - Added self.filename
;       31-Oct-2006: Harry Warren   - Modified to read compressed fits files
;       29-Sep-2007: A. Gardini - Freed pointers.
;       10-Oct-2007: A. Gardini - Added Doppler vel.
;       11-Oct-2007: A. Gardini - Alternative reading of the aux
;                                         variables (working but commented)
;        4-Dec-2007: A. Gardini - Corrected "noaaux" typo.
;        8-Jan-2008: A. Gardini - Commented FF references.
;        2-Feb-2008: Mike Marsh - Added photon calibration status, if true 
;                                 set units to photon events
;        8-Feb-2008: A. Gardini - Added check on input file.
;       20-May-2008: A. Gardini - Added retain calibration status.
;        3-Jul-2008: A. Gardini - Added init keyword.
;       16-Feb-2009: A. Gardini - Fixed units of level 1 slot data.
;       25-Jun-2009: A. Gardini - Fixed the photon unit tag.
;       26-Jun-2009: A. Gardini - Roll back of the previous change.
;       08-Apr-2010: T. Fredvik - Added slit tilt and wavelength
;                                 correction
;       22-Oct-2019: P. Young   - Added 'rwp' to calstat.
;  $Id: eis_data__readfits.pro 11677 2019-12-05 21:00:25Z pyoung $
;
;-
pro eis_data::readfits, file, hdr = hdr, error = error, noaux=noaux, init=init
  if n_elements(init) eq 1 then return

  if n_params() lt 1 then begin
    message, 'No datafile specified!', /informational
    message, 'usage: data->read, file, hdr = hdr', $
              /informational
    return
  endif
;
  if n_elements(noaux) eq 0 then noaux=0
;
  self.title='Intensity'

; open and construct header object
  if ptr_valid(self.hdr) then begin
    heap_free,self.hdr
  end
  self.hdr = ptr_new(obj_new('fits_hdr',file))
  self.filename = file

; open the file using readfits (HPW 31-OCT-2006)
  bte_data = readfits(file,headr,exten_no=1,/silent)
  if bte_data[0] eq -1 then begin
    error=-1
    return
  endif

; extract information from fits headr:
  self.nwin = fxpar(headr, 'nwin')  ; number of line windows
  self.ys = fxpar(headr, 'yws')     ; y window start pixel (cam)
  self.yw = fxpar(headr, 'yw')      ; y window width (pixels) (cam)
  self.nexp = fxpar(headr, 'nexp')  ; number of exposures
  self.nexp_prp = fxpar(headr, 'nexp_prp')
;  self.nraster = self.nexp/self.nexp_prp
  self.nraster = fxpar(headr, 'nraster')
  self.slit_id = fxpar(headr, 'slit_id')
  self.slit_ind = fxpar(headr, 'slit_ind')

; set calibration status
  self.calstat.dc = fxpar(headr,'cal_dc')
  self.calstat.hp = fxpar(headr,'cal_hp')
  self.calstat.dp = fxpar(headr,'cal_dp')
  self.calstat.wp = fxpar(headr,'cal_wp')
  self.calstat.rwp = fxpar(headr,'cal_rwp')
  self.calstat.cr = fxpar(headr,'cal_cr')
  self.calstat.abs = fxpar(headr,'cal_abs')
  self.calstat.sens = fxpar(headr,'cal_sens')
  self.calstat.phot = fxpar(headr,'cal_phot')	
  self.calstat.retain = fxpar(headr,'cal_reta')	
; self.calstat.ff = fxpar(headr,'cal_ff')
; self.calstat.wvl = fxpar(headr,'cal_wvl')
; check if absolute calibration and photon calibration has been done and set unit if so
  if (self->getcalstat()).abs and not (self->getcalstat()).phot then begin
    if (self.slit_ind eq 0 or self.slit_ind eq 2) then $
      self->setunit,['erg/cm!e2!N/s/sr/'+string("305B),'',''] $
    else self->setunit,['erg/cm!e2!N/s/sr','','']
  endif
  if (self->getcalstat()).abs and (self->getcalstat()).phot then $
      self->setunit,['Photon-Events']

; read parameters related to fits reformat:
  self.fitslev = fxpar(headr, 'data_lev')
  date_rf0 = fxpar(headr, 'date_rf0')
  self-> setfits_reformat, 'date_rf0', date_rf0
  orig_rf0 = fxpar(headr, 'orig_rf0')
  self-> setfits_reformat, 'orig_rf0', orig_rf0
  ver_rf0 = fxpar(headr, 'ver_rf0')
  self-> setfits_reformat, 'ver_rf0', ver_rf0
  date_rf1 = fxpar(headr, 'date_rf1')
  self-> setfits_reformat, 'date_rf1', date_rf1
  orig_rf1 = fxpar(headr, 'orig_rf1')
  self-> setfits_reformat, 'orig_rf1', orig_rf1
  ver_rf1 = fxpar(headr, 'ver_rf1')
  self-> setfits_reformat, 'ver_rf1', ver_rf1

  ; read fine mirror step size
  fmirr_ss = fxpar(headr, 'fmir_ss')
  if fmirr_ss eq 0 then begin
    self.obsmode = 'sit-and-stare'
    self.sit_and_stare = 1
;    self.nraster = 1
  endif else begin
    self.obsmode = 'scanning'
    self.sit_and_stare = 0
  endelse
  self.xw = 0                     ; preliminary xsize of line windows

  ; get information about data types in the BTE columns
  fxbtform, headr, tbcol, idltype, format, numval, maxval
  nauxcol=(size(tbcol))[1]-self.nwin

  ; read each binary table extension column for line windows
  for icol = 0, self.nwin-1 do begin
    istr = strtrim(string(icol+1), 2)
    ;extract line width (x) from 'TDIM' keyword
    tdim = fix(fxbtdim(fxpar(headr, 'TDIM' + istr)))
    self.xw[icol] = tdim[0]

    ; line window start pixel:
    self.xs[icol] = fxpar(headr, 'tdetx' + istr)

    ;; set up size of line window pointer
    ptr_free,self.w[icol]
    self.w[icol] = ptr_new(intarr(self.xw[icol],self.yw[icol],self.nexp))
    ;; read data using tbget (HPW 31-OCT-2006)
    d = tbget(headr,bte_data,icol+1)
    *self.w[icol] = reform(d,self.xw[icol],self.yw[icol],self.nexp)
  endfor
  if noaux then goto,finished

  ; read auxiliary data (first 12 columns)
  ptr_free,self.aux_data.ti_1
  ptr_free,self.aux_data.ti_2
  ptr_free,self.aux_data.mhc_dur
  ptr_free,self.aux_data.exp_dur
  ptr_free,self.aux_data.fmirr
  ptr_free,self.aux_data.hslstat
  ptr_free,self.aux_data.xrt_flfl
  ptr_free,self.aux_data.xrt_fl_x
  ptr_free,self.aux_data.xrt_fl_y
  ptr_free,self.aux_data.aec_hepc
  ptr_free,self.aux_data.aec_lepc
  ptr_free,self.aux_data.mhcfmsg

  icol = self.nwin
  d = tbget(headr,bte_data,icol+1)
  self.aux_data.ti_1 = ptr_new(ulonarr(self.nexp))
  *self.aux_data.ti_1 = ulong(d)
  d = tbget(headr,bte_data,icol+2)
  self.aux_data.ti_2 = ptr_new(ulonarr(self.nexp))
  *self.aux_data.ti_2 = ulong(d)
  d = tbget(headr,bte_data,icol+3)
  self.aux_data.mhc_dur = ptr_new(ulonarr(self.nexp))
  *self.aux_data.mhc_dur = ulong(d)
  d = tbget(headr,bte_data,icol+4)
  self.aux_data.exp_dur = ptr_new(uintarr(self.nexp))
  *self.aux_data.exp_dur = ulong(d)
  self.exp = ptr_new(fltarr(self.nexp))
  *self.exp = ulong(d)/1.e2  ; exposure time in seconds (exp_dur is in units of 10ms)
  d = tbget(headr,bte_data,icol+5)
  self.aux_data.fmirr = ptr_new(uintarr(self.nexp))
  *self.aux_data.fmirr = d
  d = tbget(headr,bte_data,icol+6)
  self.aux_data.hslstat = ptr_new(uintarr(self.nexp))
  *self.aux_data.hslstat = d
  d = tbget(headr,bte_data,icol+7)
  self.aux_data.xrt_flfl = ptr_new(bytarr(self.nexp))
  *self.aux_data.xrt_flfl = d
  d = tbget(headr,bte_data,icol+8)
  self.aux_data.xrt_fl_x = ptr_new(bytarr(self.nexp))
  *self.aux_data.xrt_fl_x = d
  d = tbget(headr,bte_data,icol+9)
  self.aux_data.xrt_fl_y = ptr_new(bytarr(self.nexp))
  *self.aux_data.xrt_fl_y = d
  d = tbget(headr,bte_data,icol+10)
  self.aux_data.aec_hepc = ptr_new(ulonarr(self.nexp))
  *self.aux_data.aec_hepc = d
  d = tbget(headr,bte_data,icol+11)
  self.aux_data.aec_lepc = ptr_new(ulonarr(self.nexp))
  *self.aux_data.aec_lepc = d
  d = tbget(headr,bte_data,icol+12)
  self.aux_data.mhcfmsg = ptr_new(uintarr(self.nexp))
  *self.aux_data.mhcfmsg = d
  if nauxcol eq 12 then begin
    message,'old style format fits file, no xcen, ycen, ccd or mhc temperatures',/info
    goto,finished
  endif

  ptr_free,self.aux_data.xcen
  ptr_free,self.aux_data.ycen
  d = tbget(headr,bte_data,icol+13)
  self.aux_data.xcen = ptr_new(fltarr(self.nexp))
  *self.aux_data.xcen = d
  d = tbget(headr,bte_data,icol+14)
  self.aux_data.ycen = ptr_new(fltarr(self.nexp))
  *self.aux_data.ycen = d
  if nauxcol eq 14 then begin
    message,'old style format fits file, no ccd or mhc temperatures',/info
    goto,finished
  endif

  ptr_free,self.aux_data.ccda_temp
  ptr_free,self.aux_data.ccdb_temp
  d = tbget(headr,bte_data,icol+15)
  self.aux_data.ccda_temp = ptr_new(fltarr(self.nexp))
  *self.aux_data.ccda_temp = d
  d = tbget(headr,bte_data,icol+16)
  self.aux_data.ccdb_temp = ptr_new(fltarr(self.nexp))
  *self.aux_data.ccdb_temp = d
  if nauxcol eq 16 then begin
    message,'old style format fits file, no mhc temperatures',/info
    goto,finished
  endif

  iadd=0
  if nauxcol eq 18 or nauxcol eq 21 OR nauxcol EQ 23 then begin
    ptr_free,self.aux_data.mhc_hz_t10
    ptr_free,self.aux_data.mhc_hz_t15
    d = tbget(headr,bte_data,icol+17)
    self.aux_data.mhc_hz_t10 = ptr_new(fltarr(self.nexp))
    *self.aux_data.mhc_hz_t10 = d
    d = tbget(headr,bte_data,icol+18)
    self.aux_data.mhc_hz_t15 = ptr_new(fltarr(self.nexp))
    *self.aux_data.mhc_hz_t15 = d
    iadd=2
    message,'new style format fits file, mhc temperatures',/info
  endif
  if nauxcol eq 19 or nauxcol eq 21 OR nauxcol EQ 23 then begin
    ptr_free,self.aux_data.v_sat
    ptr_free,self.aux_data.v_sun
    ptr_free,self.aux_data.v_eth
    d = tbget(headr,bte_data,icol+17+iadd)
    self.aux_data.v_sat = ptr_new(fltarr(self.nexp))
    *self.aux_data.v_sat = d
    d = tbget(headr,bte_data,icol+18+iadd)
    self.aux_data.v_sun = ptr_new(fltarr(self.nexp))
    *self.aux_data.v_sun = d
    d = tbget(headr,bte_data,icol+19+iadd)
    self.aux_data.v_eth = ptr_new(fltarr(self.nexp))
    *self.aux_data.v_eth = d
    iadd += 3
    message,'new style format fits file, Doppler vel.',/info
  endif
  
  IF nauxcol EQ 23 THEN BEGIN
     ptr_free,self.aux_data.hkpixcorrtilt
     ptr_free,self.aux_data.hkpixcorrtime
     d = tbget(headr,bte_data,icol+17+iadd)
     self.aux_data.hkpixcorrtilt = ptr_new(fltarr(self.nexp))
     *self.aux_data.hkpixcorrtilt = d
     d = tbget(headr,bte_data,icol+18+iadd)
     self.aux_data.hkpixcorrtime = ptr_new(fltarr(self.nexp))
     *self.aux_data.hkpixcorrtime = d
     message,'new style format fits file, wavelength correction w/HK temperatures.',/info 
  ENDIF
  
; alternatively
;
;  if nauxcol gt 16 then begin
;    ptr_free,self.aux_data.mhc_hz_t10
;    ptr_free,self.aux_data.mhc_hz_t15
;    ptr_free,self.aux_data.v_sat
;    ptr_free,self.aux_data.v_sun
;    ptr_free,self.aux_data.v_eth
;    for iaux = 17,nauxcol do begin
;      ttypenr = self.nwin + iaux
;      ttypestr = 'TTYPE' + strtrim(string(ttypenr),2)
;      column_name = fxpar(headr, ttypestr)
;      d = tbget(headr,bte_data,icol+iaux)
;      case column_name of
;        'MHC_HZ_T10': begin
;         self.aux_data.mhc_hz_t10 = ptr_new(fltarr(self.nexp))
;         *self.aux_data.mhc_hz_t10 = d
;         end
;        'MHC_HZ_T15': begin
;         self.aux_data.mhc_hz_t15 = ptr_new(fltarr(self.nexp))
;         *self.aux_data.mhc_hz_t15 = d
;         end
;        'V_SAT   ': begin
;         self.aux_data.v_sat = ptr_new(fltarr(self.nexp))
;         *self.aux_data.v_sat = d
;         end
;        'V_SUN   ': begin
;         self.aux_data.v_sun = ptr_new(fltarr(self.nexp))
;         *self.aux_data.v_sun = d
;         end
;        'V_ETH   ': begin
;         self.aux_data.v_eth = ptr_new(fltarr(self.nexp))
;         *self.aux_data.v_eth = d
;         end
;        else : print,'Column name mismatch'
;      endcase
;    endfor
;  endif

  finished:

end