pro read_sotsp, files, index,data, _extra=_extra, scan_info=scan_info, phead=phead, slitmap=slitmap, mapxy=mapxy
;
;+
;   Name: read_sotsp
;
;   Purpose: ssw style api for SOT-SP level2 and level 2.1
;
;   Input Parameters:
;      files - SOT-SP Level 2 (merlin/binary extended fits)
;
;   Output Parameters:
;      index,data - ssw compatible "index,data", 1:1 sp maps
;
;   Keyword Parameters:
;      _extra - keywords -> sotsp_l2struct2data via inheritance (e.g., /XYCOORD)
;                        -> mrdfits for "non-merlin"
;      scan_info (output) -> structure containg 1D vectors
;                           (TIMES, SLITPOS, and SCATT)
;      phead (output) -> primary FITS header(structure)
;      slitmap (output) -> for "new" version ("level 2.01") of interpolated
;                          level 2 data, provides an array indicating which
;                          scan columns have been interpolated (output of
;                          sotsp_rasterize.pro)
;      mapxy (output) -> for "new" version ("level 2.01") of interpolated
;                        level 2 data, indicates (x,y) in arcsec of center of
;                        interpolated image(s)
;      silent (switch) -> inherit via mrdfits
; 
;   Calling Sequence:
;      IDL> read_sotsp, l2file, index, data , scan_info=scan_info
;
;   Calling Example:
;      IDL> read_sotsp,l2,index,data,scan_info=si
;      IDL> help,data,si,index,/str
; DATA            FLOAT     = Array[4156, 256, 36] <<<DATA 
; ** Structure <2663bd4>, 3 tags, length=66944, data length=66944, refs=1:
;   TIMES           STRING    Array[4156]  <<< times per slit position 
;   SLITPOS         FLOAT     Array[4156]  <<< slit pos 
;   SCATTPRO        FLOAT     Array[112]   <<< scattered light profile
; ** Structure <2930804>, 54 tags, length=624, data length=620, refs=1:
;   SIMPLE          INT              1      <<< index.simple
;   BITPIX          LONG                 8  <<< index... etc 
;   (...etc; more tags here...)
;   FITSTAGS        STRING    'SFLD'        <<< ssw/info added info
;   CUNIT1          STRING    'kG'
;   CUNIT2          STRING    'kG'
;   FTYPE           STRING    'Magnetic Field'
;   INSTRUME        STRING    'SOT/SP Magnetic Field' <<< ssw superset added
;   TELESCOP        STRING    'HINODE'
;   CDELT1          FLOAT          0.148565
;   CDELT2          FLOAT          0.159991
;   DATE_OBS        STRING    ' 1-Dec-2006 11:00:47.000'
;   DATE_END        STRING    ' 1-Dec-2006 16:59:52.000'
;   
;
;   History:
;      1-Jun-2008 - S.L.Freeland - ssw interface to merlin/bin extended FITS
;     25-aug-2016 - S.L.Freeland - train to read Ted/Phil style SP "movie"
;                                  cubes, add PHEAD output
;     26-Apr-2018 - PGShirts - added two if/then catches to reform 'dat' if
;                              n_dimensions(dat) eq 1 (this is the
;                              sit-and-stare case).  NOTE: as each column is
;                              read in a for-loop, the read for a
;                              sit-and-stare observation can take a long time.
;                              --perhaps this could be vectorized in the
;                              future to improve the sit-and-stare obs ingest
;                              speed.
;     14-Jul-2020 - M.DeRosa - added capability to read in level 2.1
;                              (disambiguated and reprojected) data if it
;                              happened to be retrieved by sotsp_getdata.pro. 
;     12-Mar-2025 - M.DeRosa - added ability to read in "new" interpolated data
;                              files, and distinguish between the two formats,
;                              and also added slitmap keyword
;-
;

if not file_exist(files[0]) then begin
  box_message,'Require existing & local SOTSP Level2 filename'
  return
endif

mreadfits,files[0],phead
next=get_fits_nextend(files[0]) ; get #nextesions

if required_tags(phead,'invcode') then begin  ;  assume MERLIN output

  ;  read lev2 header/images into the standard SSWIDL index/data arrays
  l2s=merlin_read_fits2(files[0],_extra=_extra)
  sotsp_l2struct2data, l2s, index, data , _extra=_extra, scan_info=scan_info

  ;  read in lev2.1 (disambiguated data) if available
  lev21file=str_replace(files[0],'level2hao','level2.1hao')
  lev21file=str_replace(lev21file,'.fits','_L2.1.fits')
  if file_test(lev21file) then begin

    lev21s=ssw_binfits2struct(lev21file,_extra=_extra)
    lev21s=join_struct({header:headfits(lev21file)},temporary(lev21s))
    sotsp_l21struct2data,lev21s,index21,data21,_extra=_extra,scan_info=scan_info

    ;  are all lev2 header (=index) keywords in the lev2.1 headers?
    lev2tags=tag_names(index[0])
    flag=0
    for ii=0,n_elements(lev2tags)-1 do begin
      if ~tag_exist(index21[0],lev2tags[ii]) then begin
        print,'  tag '+lev2tags[ii]+' not in lev2.1 header'
        flag++
      endif
    endfor
    if flag gt 0 then begin
      box_message,'Found '+strtrim(flag,2)+' tags from lev2 header that are not in lev2.1 header, returning only lev2 data'
      return
    endif
    
    ;  next, extend lev2 headers to include new lev2.1 keywords
    lev21tags=tag_names(index21[0])
    for ii=0,n_elements(lev21tags)-1 do $
      if ~tag_exist(index,lev21tags[ii]) then $
        index=add_tag(temporary(index),index21[0].(ii),lev21tags[ii])
    
    ;  now should be able to concatenate index structures and data arrays
    index=merge_struct(index,index21)
    data=[[[temporary(data)]],[[temporary(data21)]]]
    
  endif
  
endif else begin  ;  non-MERLIN data file
  dat=mrdfits(files[0],1,head,_extra=_extra)
  headstruct=fitshead2struct(head)
  is_repeating=headstruct.naxis eq 3

  ;  check whether we have an interpolated lev2 data file (old or new)
  is_oldinterp=0b
  is_newinterp=0b
  if required_tags(headstruct,'LVL_NUM') then begin
    if abs(headstruct.lvl_num-1.5) lt 1e-5 then is_oldinterp=1b 
    if abs(headstruct.lvl_num-2.01) lt 1e-5 then is_newinterp=1b 
  endif
  if required_tags(headstruct,'DATA_LEV') then begin
    if abs(headstruct.data_lev-1.5) lt 1e-5 then is_oldinterp=1b 
    if abs(headstruct.data_lev-2.01) lt 1e-5 then is_newinterp=1b 
  endif

  ;  adjust number of extensions if new-style file
  if is_newinterp then nimext=next-2 else nimext=next

  ;  create index array
  index=replicate(fitshead2struct(head),nimext)

  ; in case of single slit position (where extension is a 1-d array)
  if n_dimensions(dat) eq 1 then dat=reform(dat,[1,n_elements(dat)],/overwr)
  
  ;  create data array
  data=make_array([size(dat,/dim),nimext],type=size(dat,/type))

  ;  loop through extensions
  for e=0,nimext-1 do begin
    delvarx,dat,head
    dat=mrdfits(files[0],e+1,head,_extra=_extra)
    
    ; in case of single slit position (where extension is a 1-d array)
    if n_dimensions(dat) eq 1 then dat=reform(dat,[1,n_elements(dat)],/overwr)

    ;  populate data array
    if is_repeating then data[0,0,0,e]=dat else data[0,0,e]=dat
    index[e]=fitshead2struct(head)
    
  endfor
  
  ;  fill slitmap keyword for output (new files only)
  if is_newinterp then begin
    delvarx,dat,head
    slitmap=mrdfits(files[0],e+1,head,_extra=_extra)
    mapxy=mrdfits(files[0],e+2,head,_extra=_extra)
  endif
        
endelse

return
end