;+
;
; PROJECT:
; RHESSI
;
; NAME:
; read_xsm_4_ospex
;
; PURPOSE:
; Read XSM data from a FITS file into Spex internal variables.
; This procedure is based on read_hessi_4_ospex.
;
; CATEGORY:
; SPEX
;
; CALLING SEQUENCE:
; read_xsm_4_ospex, FILES=files, data_str=data_str, ERR_CODE=err_code,
;                   ERR_MSG=err_msg, _REF_EXTRA=_ref_extra
;
; INPUTS:
;
; INPUT KEYWORDS:
; FILES - Name of file to read.
;
; OUTPUTS:
;
; OUTPUT KEYWORDS:
; The following keywords are set to values needed by Spex.  Many are described
; in the Spex documentation.  Those not described are internal variables not
; meant to be used / accessed by the user.  Those keywords listed in the
; procedure definition line not listed here are included for Spex comptability,
; and are not used by this procedure.
;   START_TIME
;   END_TIME
;   RCOUNTS
;   ERCOUNTS
;   UT_EDGES
;   UNITS
;   AREA
;   LTIME
;   CT_EDGES
;   TITLE
;   RESPFILE
;   DETUSED
;   ATTEN_STATES
;
; OPTIONAL OUTPUTS:
;
; COMMON BLOCKS:
; read_spex_com, data, ut_data, data_file_read, read_start, read_end
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; EXAMPLE:
; Not to be used directly.  Called from within SPEX to load XSM data
;
; CALLS:
; break_file, concat_dir, get_edges, loc_file, xsm_fits2spectrum
;
; Written  16-Nov-2004 - Sandhia Bansal
;
; Modification History:
; Deleted some comments.  - Sandhia Bansal - 12/08/04
; 24-jun-2005, fixed normalization for spectrum input.  Those data
;   are already in counts
; 11-Aug-2005, Kim.  Set detused to 'XSM' by default
; 13-May-2008, Kim. Remove any energies < 1 keV from data and energy arrays
; 25-Jun-2015, Kim. Some files (e.g. 31-jul-2004 from shyama) had excess energy bins all set to [20.,20.] keV.
;   These extra bins caused a crash, so now remove bins whose width is 0. (see corresponding change in spex_xsm_fits2drm)
; 12-Oct-2017, Kim. Removal of energy bins < 1 keV and with width=0. was not being applied to ltime array. Fixed.
 
;-
;------------------------------------------------------------------------------
PRO read_xsm_4_ospex, FILES=files, $
                        data_str=data_str, $


                        ERR_CODE=err_code, $
                        ERR_MSG=err_msg, $
                        _REF_EXTRA=_ref_extra

data_str = -1

err_code = 0
err_msg = ''

respfile = ''

delta_light= 1.0

if files[0] eq '' then begin
    err_code=1
    err_msg = 'No spectrum file selected.'
    return
endif

break_file, files[0], disk, dir, fname, fext

fits = Strupcase( fext ) EQ '.FITS'

atten_states = -1

; Expect XSM spectral file to be in FITS format.
IF fits THEN BEGIN
    ; Read spectrum from the file.
    xsm_fits2spectrum, FILE=files[0], $
                   PRIMARY_HEADER=p_hdr, $
                   EXT_HEADER=sp_hdr, $
                   EXT_DATA=sp_data, $
                   ENEBAND_DATA=en_data, $
                   RESPFILE=respfile, $

                   _EXTRA=_ref_extra, $
                   ERR_MSG=err_msg, $
                   ERR_CODE=err_code, /silent

    IF err_code THEN RETURN

    ut_edges = Dblarr( 2, N_Elements(sp_data) )
    n = 0
          ; This is the new, compliant formulation for time in the RATE files.
        ; XSM file can be in one of the two FITS formats:
        ;   1:
        ;      - delta times are specified in a column or keyword named INTEGRATION_TIME
        ;      - start times are specified in a column or keyword named START_OBS
        ;   2:
        ;      - delta times are specified in a column or keyword named TIMEDEL
        ;      - start times are specified in a column or keyword named TIME
        IF (tag_exist(sp_data, 'INTEGRATION_TIME')) THEN BEGIN   ; column exists
           timedel = sp_data.integration_time
        ENDIF ELSE BEGIN
           val = fxpar(sp_hdr, 'INTEGRATION_TIME', count=n)  ; try reading from keyword
           IF (n GT 0) THEN BEGIN   ; keyword exists
              timedel = float(val)
           ENDIF ELSE BEGIN   ; try reading from TIMEDEL column/keyword
              timedel = tag_exist(sp_data, 'TIMEDEL')? $
                  sp_data.timedel : float( fxpar(sp_hdr, 'TIMEDEL'))
           ENDELSE
       endelse

        IF (tag_exist(sp_data, 'T_UTC')) THEN BEGIN

          ut_edges=anytim(sp_data.t_utc,/sec)
          ut_edges=[ut_edges, last_item(ut_edges)+last_item(timedel)]
          ut_edges= get_edges(ut_edges, /edges_2)



        ENDIF


    ct_edges = Transpose( [ [ en_data.e_min ] , [ en_data.e_max ] ] )


ENDIF ELSE BEGIN
    print, 'Invalid format for this file'
    return
ENDELSE

sp_tags = Tag_Names( sp_data )
data_tags = ['SPECTRUM', 'RATE']
ercounts_idx = -1
FOR i=0, N_Elements( sp_tags )-1L DO BEGIN
    match = Where( sp_tags[i] EQ data_tags, n_match )
    IF n_match GT 0 THEN BEGIN
        rcounts_idx = i
        unit = data_tags[ match[ 0 ] ]
    ENDIF
    ; Expect either ERATE or STAT_ERR for a column containing errors
    ematch = Where( sp_tags[i] EQ 'E' + data_tags, n_ematch ) or $
             Where (sp_tags[i] EQ 'STAT_ERR', n_ematch )
    IF n_ematch GT 0 THEN ercounts_idx = i
ENDFOR

rcounts = sp_data.( rcounts_idx )
IF (ercounts_idx[0] GT -1) THEN BEGIN
   ercounts = sp_data.(ercounts_idx)
ENDIF ELSE BEGIN
   ; Calculate errors by taking sqrt of data.  Set all errors to 0 (zero) where
   ;    data elements are negative.
   ercounts = sp_data.(rcounts_idx)
   ind = where(sp_data.(rcounts_idx) LT 0, count)
   if (count GT 0) THEN $
      ercounts[ind] = 0.0

   ercounts = sqrt(ercounts)
ENDELSE


;ercounts = ercounts_idx[0] GT -1 ? sp_data.(ercounts_idx) : sqrt(sp_data.( rcounts_idx ))
;ercounts = ercounts_idx[0] GT -1 ? sp_data.(ercounts_idx) : sqrt(sp_data.spectrum)

area = fxpar(sp_hdr, 'GEOAREA', COUNT=area_cnt)
IF area_cnt GT 0 THEN BEGIN
   area = (st2num( area, area_stat ))[0]
   IF 1 - area_stat THEN area = 1.
ENDIF ELSE $
   area = 1

; Make any 3-D arrays into 2-D - multiple detector sets must be combined for
; analysis
IF n_dimensions( rcounts ) EQ 3 THEN rcounts = Total( rcounts, 2 )

IF n_dimensions( ercounts ) EQ 3 THEN ercounts = total( ercounts, 2 )

default, ltime, rcounts*0.0  + 1.0
;The next line is an approximation, it needs refinement
;That's why this shouldn't be done here.
IF n_dimensions( ltime ) EQ 3 THEN ltime = total( ltime,2 )

nbin = n_elements( ut_edges )/2
nchan = n_elements( ct_edges )/2

; transpose the rcounts, ercounts, and ltime matrices, if necessary
delta_light = get_edges( ct_edges, /width ) # ( 1+fltarr(nbin) )
acctime = float( get_edges( ut_edges, /width)) ## ( 1+fltarr(nchan) )

IF n_elements( ltime ) EQ 1 THEN $
  ltime = ltime[0]

IF n_elements( ltime ) NE n_elements(rcounts) THEN $
  ltime = ltime ## (1+fltarr(nchan))

CASE unit OF
    'FLUX': BEGIN
        rcounts = rcounts * area * acctime * ltime * delta_light
        ercounts = ercounts * area * acctime * ltime * delta_light
    END
    'RATE': BEGIN
        rcounts = rcounts * acctime * ltime
        ercounts = ercounts * acctime * ltime

    END
    'SPECTRUM': BEGIN
;        rcounts = rcounts
;        ercounts = ercounts

    END
    ELSE: ;; It's already counts.
ENDCASE

ltime = ltime * acctime

units = ' s!u-1!n cm!u-2!n keV!u-1!n'
wchan = Lindgen( nchan )

start_time = anytim( fxpar( p_hdr, 'DATE_OBS', count=c ), /vms )
if c eq 0 then start_time = anytim( fxpar( p_hdr, 'DATE-OBS' ), /vms)
end_time = anytim( fxpar( p_hdr, 'DATE_END', count=c ), /vms )
if c eq 0 then end_time = anytim( fxpar( p_hdr, 'DATE-END' ), /vms)

title = 'XSM SPECTRUM'
default, arffile, ''
if not keyword_set(respfile) then begin
    break_file, files[0], rdisk, rdir, rfname, rfext
    s = strmid(rfname, 3, strlen(rfname))
    respfile = $
       concat_dir(concat_dir(rdisk, rdir), ('TEST_'+strmid(s,0,6)+strmid(s,8,2)+'.RMF'))
    endif
f = loc_file(respfile, COUNT=count, LOC=loc)

IF (f EQ '') THEN BEGIN
   err_msg = ['FILE ' + files[0]+':', $
                   'has RESPFILE: ' + respfile, $
                   'RESPFILE: ' + respfile + ' not found.' ]
   ;err_code = 1 - not a fatal error
   if not spex_get_nointeractive() then xmessage,err_msg else print,err_msg
ENDIF


detused = fxpar( sp_hdr, 'DETUSED', COUNT=rcount )
if rcount eq 0 then detused = 'XSM' else begin
    detused = str_replace (detused, 'SEGMENTS: ', '')
    detused = str_replace (detused, '|', ' ')
endelse

; don't keep any energies below 1. keV - chianti can't handle
; note: the drm array also has to be reduced to eliminate energies < 1 keV
; 25-jun-2015, also added requirement that width of energy bin be > 0.
wid = get_edges(ct_edges, /width)
q = where (ct_edges[0,*] gt 1. and wid gt 0.)
;e1 = q[0]
;ct_edges = ct_edges[*,e1:*]
;rcounts = rcounts[e1:*,*]
;ercounts = ercounts[e1:*,*]
ct_edges = ct_edges[*,q]
rcounts = rcounts[q,*]
ercounts = ercounts[q,*]
ltime = ltime[q,*]  ; added 12-oct-2017, Kim

data_str = { $
    START_TIME: start_time, $
    END_TIME: end_time, $
    RCOUNTS: rcounts, $
    ERCOUNTS: ercounts, $
    UT_EDGES: ut_edges, $
    UNITS: units, $
    AREA: area, $
    LTIME: ltime, $
    CT_EDGES: ct_edges, $
    TITLE: title, $
    RESPFILE: respfile, $

    detused: detused, $
    atten_states: atten_states }

END