PRO xrt_mreadfits, files, info, data, noscale=noscale,                    $
                   strtemplate=strtemplate, nodata=nodata, header=header, $
                   quiet=quiet, outsize=outsize, maxx=maxx, maxy=maxy,    $
                   add_standard=add_standard, comsep=comsep,              $
                   hissep=hissep, comments=comments, history=history,     $
                   ccnts=ccnts, hcnts=hcnts, nocomments=nocomments,       $
                   nohistory=nohistory, all_keywords=all_keywords,        $
                   nofill=nofill, silent=silent, wcs=wcs, _extra=_extra

; =========================================================================
;+
;
; PROJECT:
;       Solar-B / XRT
;
; NAME:
;       XRT_MREADFITS
;
; CATEGORY:
;       File I/O
;
; PURPOSE:
;       This routine is simply a wrapper to <mreadfits.pro>.
;       The object is to read a combination of images of different sizes,
;       without rescaling any of them. Instead, the output data-cube 
;       dimensions are set to the maximum size, and smaller images are 
;       inset to their lower-left corner. All header tags should still be 
;       accurate for the actual image, ignoring any buffer pixels.
;
; CALLING SEQUENCE:
;
;       XRT_MREADFITS, files, info, data, 
;                      [,/noscale],       [,strtemplate=strtemplate] 
;                      [,/nodata]         [,header=header] 
;                      [,/quiet]          [,outsize=outsize] 
;                      [,maxx=maxx]       [,maxy=maxy]
;                      [,/add_standard]   [,/comsep]
;                      [,/hissep]         [,comments=comments] 
;                      [,history=history] [,ccnts=ccnts] 
;                      [,hcnts=hcnts]     [,/nocomments]
;                      [,/nohistory]      [,/all_keywords]
;                      [,/nofill]         [,/silent] 
;                      [,wcs=wcs]         [,_extra=_extra]
;
; INPUTS:
;       (See <mreadfits.pro>. 
;       All arguments map directly into the call to that routine.)
;
; KEYWORDS:
;       (See <mreadfits.pro>. 
;       All arguments map directly into the call to that routine.)
;
; OUTPUTS:
;       (See <mreadfits.pro>. 
;       All arguments map directly into the call to that routine.)
;
; EXAMPLES:
;
;       Simple case: get header and data from an array of FITS filenames.
;       IDL> xrt_mreadfits, filenames, header, data
;
; COMMON BLOCKS:
;       None.
;
; NOTES:
;
;       1) This wrapper routine should be used to read XRT FITS files,
;          rather than using <mreadfits.pro> directly. This ensures
;          that the data and the header tags remain consistent.
;
; CONTACT:
;
;       Comments, feedback, and bug reports regarding this routine may be
;       directed to this email address:
;                xrt_manager ~at~ head.cfa.harvard.edu
;
; MODIFICATION HISTORY:
;
progver = 'v2006.Nov.30' ;--- (MW) Written.
;
;
;-
; =========================================================================


  if n_elements(files) eq 1 then begin

    mreadfits, files, info, data, noscale=noscale,                    $
               strtemplate=strtemplate, nodata=nodata, header=header, $
               quiet=quiet, outsize=outsize, maxx=maxx, maxy=maxy,    $
               add_standard=add_standard, comsep=comsep,              $
               hissep=hissep, comments=comments, history=history,     $
               ccnts=ccnts, hcnts=hcnts, nocomments=nocomments,       $
               nohistory=nohistory, all_keywords=all_keywords,        $
               nofill=nofill, silent=silent, wcs=wcs, _extra=_extra

  endif else begin

    mreadfits, files, info0, /nodata, /quiet


    uniq_naxis1 = info0[uniq(info0.naxis1, sort(info0.naxis1))].naxis1
    n1 = n_elements(uniq_naxis1)

    uniq_naxis2 = info0[uniq(info0.naxis2, sort(info0.naxis2))].naxis2
    n2 = n_elements(uniq_naxis2)

    ss_stack = 0L
    info = info0[0]
    data = intarr(8,8)

    for ii = 0,(n1-1) do begin
      for jj = 0,(n2-1) do begin

        ss = where((info0.naxis1 eq uniq_naxis1[ii]) and $
                   (info0.naxis2 eq uniq_naxis2[jj]), count)

        if (count gt 0) then begin

          mreadfits, files[ss], info1, data1, noscale=noscale,           $
                  strtemplate=strtemplate, nodata=nodata, header=header, $
                  quiet=quiet, outsize=outsize, maxx=maxx, maxy=maxy,    $
                  add_standard=add_standard, comsep=comsep,              $
                  hissep=hissep, comments=comments, history=history,     $
                  ccnts=ccnts, hcnts=hcnts, nocomments=nocomments,       $
                  nohistory=nohistory, all_keywords=all_keywords,        $
                  nofill=nofill, silent=silent, wcs=wcs, _extra=_extra

          ss_stack = [ss_stack, ss]
          info = concat_struct(temporary(info), info1)

          if (not keyword_set(nodata)) then begin

            Nx_out = n_elements(data[*,0,0])
            Ny_out = n_elements(data[0,*,0])
            Nz_out = n_elements(data[0,0,*])
            Nx1    = n_elements(data1[*,0,0])
            Ny1    = n_elements(data1[0,*,0])
            Nz1    = n_elements(data1[0,0,*])

            ;=== First do x-dimension.
            ;=== If data (x) is bigger:
            if (Nx_out gt Nx1) then begin
              aa0 = Nx_out - Nx1
              data1 = [temporary(data1), bytarr(aa0,Ny1,Nz1)]
              Nx1 = n_elements(data1[*,0,0])

            endif else begin
              ;=== Else if data1 (x) is bigger:
              if (Nx_out lt Nx1) then begin
                aa0 = Nx1 - Nx_out
                data = [temporary(data), bytarr(aa0,Ny_out,Nz_out)]
                Nx_out = n_elements(data[*,0,0])

              ;=== Else data (x) and data1 (x) are same size.
              endif else begin
              endelse
            endelse

            ;=== Now do y-dimension.
            ;=== If data (y) is bigger:
            if (Ny_out gt Ny1) then begin
              aa0 = Ny_out - Ny1
              data1 = [[temporary(data1)], [bytarr(Nx1,aa0,Nz1)]]
              Ny1 = n_elements(data1[0,*,0])

            endif else begin
              ;=== Else if data1 (y) is bigger:
              if (Ny_out lt Ny1) then begin
                aa0 = Ny1 - Ny_out
                data = [[temporary(data)], [bytarr(Nx_out,aa0,Nz_out)]]
                Ny_out = n_elements(data[0,*,0])

              ;=== Else data (y) and data1 (y) are same size.
              endif else begin
              endelse
            endelse

            ;=== Concatenate data1 onto data.
            data = [[[temporary(data)]], [[data1]]]

          endif

        endif

      endfor
    endfor


    ss_stack = ss_stack[1:*]
    info = temporary(info[1:*])
    if (not keyword_set(nodata)) then data = temporary(data[*,*,1:*])


    ss_sort = sort(ss_stack)

    info = temporary(info[ss_sort])
    if (not keyword_set(nodata)) then data = temporary(data[*,*,ss_sort])

  endelse


END ;======================================================================