pro read_sot, sotfiles, index, data, _extra=_extra, no_patch=no_patch, $
   debug=debug , silent=silent, quiet=quiet
; 
;+ 
;   Name: read_sot
;
;   Purpose: read sot FITS files -> "index,data" w/header-index tweaks if any
;
;   Input Paramters:
;      sotfiles - list of one or more sot files to read
;
;   Output Parameters:
;        index - output header index vector (adjusted if required)
;       data - 2D/3D output array
;
;   Keyword Parameters:
;      _extra - passed to mreadfits.pro via inheritance
;     no_patch - if set, bypass attempts to find&fix header problems
;
;   History:
;      9-March-2006 - S.L.Freeland
;      1-Aug-2006 - S.L.Freeland - name change fpp->sot
;      4-oct-2006 - S.L.Freeland - turn off auto-patch for post launch data
;                   use mrdfits for SP / 4 dimensional data
;     16-oct-2006 - S.L.Freeland - nindex=ndata for 4 dimen
;      9-apr-2007 - S.L.Freeland - allow 3D SP4D reads
;     15-jun-2007 - S.L.Freeland - _extra->mreadfits for naxis>2 cases also
;     19-Feb-2010 - Y. Katsukawa - allow reading 3 iamges data (e.g. FGIVDG)
;     22-oct-2014 - S.L.Freeland - add SP3D (hao) support 
;
;   Restrictions:
;      cannot mix #axis on single call e.g., cannot mix simple FG + SP4D
;      but multi-files of either type are ok on a given call
;- 
debug=keyword_set(debug)
silent=keyword_set(silent) or keyword_set(quiet)
loud=1-silent

sot_types=str2arr('AC,CL,CR,FG,FGRAW,FGSIQUV,SP4D,SRAW') ; obs-depend tweaks? 

axis_gt2=where(strmatch(sotfiles,'*SP?D*') or $ ; SP4D/SP3D (latter added 22-oct-2014)
               strmatch(sotfiles,'*FG*IQUV*') or $
               strmatch(sotfiles,'*FG*IV*'),spcnt)
if spcnt gt 0 and n_params() gt 2 then begin 
      nfiles=spcnt
      if spcnt gt 1 then $
         if loud then box_message,['File list includes '+ strtrim(spcnt,2) + $
              ' NAXIS>2 file'+(['','s'])(spcnt gt 1) ]
      mreadfits,sotfiles(axis_gt2),index,_extra=_extra ; header->structure
      naxis=gt_tagval(index,/naxis)
      nind=n_elements(index)
      case 1 of 
         nind eq nfiles:  
         nind eq (2*nfiles): index=index(indgen(nfiles)*2)
         nind eq (3*nfiles): index=index(indgen(nfiles)*3)
         nind eq (4*nfiles): index=index(indgen(nfiles)*4)
         else: box_message,'Warning: may be mismatch #index:#images
      endcase
      aaxis=all_vals(naxis)
      if n_elements(aaxis) gt 1 then begin 
          box_message,'You cannot mix NAXIS files in one call (ex IV + QUIV)'
          return
      endif
      data=mrdfits(sotfiles(axis_gt2(0)),silent=silent) ; use mrdfits
      if n_elements(sotfiles) gt 1 then begin ; multi 3+D data merge...
         sdata=size(data)
         osize=last_nelem(sdata*nfiles)
         sdata(0)=sdata(0)+1
         sdata=[sdata(0:sdata(0)-1),nfiles,data_chk(data,/type),osize]
         temp=temporary(data)
         data=make_array(size=sdata) ; make output array (multi-4D files ok)
         data(0,0,0,0,0)=temp
         ndim=sdata(0)
         for i=1,nfiles-1 do begin
            temp=mrdfits(sotfiles(axis_gt2(i)),silent=silent) ;
            data(0,0,$
               ([0,i])(ndim eq 3),([0,i])(ndim eq 4),([0,i])(ndim eq 5))=temp 
            if debug then stop,'image '+ string(i)
         endfor
      endif
      if debug then stop,'>2 axis case'
      return ; !!!!! Early exit for NAXIS>2 !!!!!
endif 

mreadfits,sotfiles,index,data,_extra=_extra,nodata=n_params() lt 3 
no_patch=keyword_set(no_patch) or ssw_deltat(index(0).date_obs, ref='1-aug-2006') gt 0
patchit=1-keyword_set(no_patch) 



; check for "required" keywords and tweak if possible
if not required_tags(index,'crpix1,cdelt1,xcen,ycen,cunits1,crota1')  $
   and patchit then begin
   box_message,'Missing tags? - attempting patch

;  ------ plate scale ----
   if not tag_exist(index(0),'CDELT1') and $
          tag_exist(index(0),'XSCALE') then begin 
      index=add_tag(index,index.xscale,'CDELT1')
      index=add_tag(index,index.yscale,'CDELT2')
      index=add_tag(index,0l,'CROTA1')
      update_history,index,/caller,"Missing CDELTx"
   endif 

;  ----- pointing ---
   ncrpix=['CRPIX1','CRPIX2']
   hcen  =['XCEN','YCEN']
   for i=0,1 do begin
      if not tag_exist(index(0),ncrpix(i)) and $
         tag_exist(index(0),hcen(i))  then begin
         crpixn=comp_fits_crpix(gt_tagval(index,hcen(i)),$
             gt_tagval(index,'cdelt'+strtrim(i+1,2)), $
             gt_tagval(index,'naxis'+strtrim(i+1,2)))
         index=add_tag(temporary(index),crpixn,ncrpix(i))
         update_history,index,/caller,"Missing "+ncrpix(i)
      endif
   endfor
endif

return
end