pro rd_ar, infil, ss, index, data, roadmap, ref_image=refimage, $
             korrel=korrel, quiet=quiet, latitude=lat, longitude=lon, $
             date=date, size_ffi=size_ffi,ar_size=ar_size,missing=missing, $
             compress=compress,nocorrections=nocorrections, $
             norm_fact=norm_fact,sfd=sfd


  if n_params() LT 2 then begin
     print,'RD_AR: rd_ar,infil,ss,index,data,roadmap    - OR - '
     print,'       rd_ar,index,data'
     return
  endif

  ; Do we need to read in the data?
  if n_params() EQ 2 then noread = 1 else noread = 0

  if noread then begin
     roadmap = infil
     infil = 0  ; save some memory
  endif

  ARCSEC_PER_PIX = gt_pix_size()
  if n_elements(roadmap) LE 0 and (NOT noread) then rd_roadmap,infil,roadmap
  if n_elements(ss) LE 0 and (NOT noread) then ss=indgen(n_elements(roadmap))
  if (NOT noread) then pftype = gt_pfi_ffi(roadmap(ss)) MOD 2 $
  else pftype = gt_pfi_ffi(roadmap) MOD 2
  pfi = where(pftype EQ 0)
  ffi = where(pftype EQ 1)

  if noread then begin
     if n_elements(ss) EQ 0 or n_elements(roadmap) EQ 0 then $
        message,'data and/or index suppied is empty'
     message,/info,'Using supplied data and index - not reading data'
     if pfi(0) GE 0 then begin
        data_pfi = ss(*,*,pfi)
        index_pfi = roadmap(pfi)
     endif
     if ffi(0) GE 0 then begin
        indata_ffi = ss(*,*,ffi)
        inindex_ffi = roadmap(ffi)
     endif
     ss = 0
     roadmap=0
  endif

  if pfi(0) GE 0 then begin
     if (NOT noread) then begin  ; Read in the PFI data
        if NOT keyword_set(quiet) then message,'Reading PFI data',/info
        dset = mk_dset_str(infil,ss(pfi))
        rd_sda, infil, dset, index_pfi, data_pfi
        data_pfi = ass_or(data_pfi,index_pfi,subs,/only_full)
        index_pfi = index_pfi(subs)
     endif

     sdpfi = size(data_pfi)

     if sdpfi(0) LE 1 then pfi=-1 $   ; No data selected
     else begin
        if sdpfi(n_elements(sdpfi)-2) EQ 1 then begin ; decompress byte data
           compressed = where(gt_comp(index_pfi) EQ 0)
           if compressed(0) GE 0 then begin  ; Uncompress the data
              uncdata = sxt_decomp(data_pfi(*,*,compressed))
              data_pfi = fix(data_pfi)
              data_pfi(*,*,compressed) = uncdata
              uncdata = 0  ; save the memory
             ; message,/info,'Decompressed ' + $
             ;               strcompress(string(n_elements(compressed))) + $
             ;               ' PFI images.'
           endif
        endif
        if NOT keyword_set(nocorrections) then begin
           data_pfi = dark_sub(index_pfi,temporary(data_pfi))
           data_pfi = leak_sub(index_pfi,temporary(data_pfi))
        endif 

        data_pfi = align_ar(data_pfi,index_pfi,ref_image=ref_image, $
                            ref_sunc=refsc,ref_xy=refxy, ref_ll=refll, $
                             ref_day=ref_day,missing=missing,ar_size=ar_size, $
                             center_lat=lat,center_lon=lon,center_date=date)
        sdpfi = size(data_pfi)
        ;if keyword_set(compress) and sdpfi(sdpfi(0)+1) NE 1 then $
        ;  data_pfi=sxt_comp(temporary(data_pfi))
     endelse
  endif 
  if pfi(0) EQ -1 then  begin
     message,'No PFI data',/info
     if NOT (n_elements(lat) GT 0 and n_elements(lon) GT 0) then $
        message,'Lat/Lon must be specifed in the absence of PFI data'
     if NOT (n_elements(date) GT 0) then $
        message,'Date must be specified in the absence of PFI data'
  endif

  if ffi(0) NE -1 then begin
     if NOT keyword_set(quiet) and (NOT noread) then $
        message,'Reading FFI data',/info
     nimage_ffi = n_elements(ffi)
     if pfi(0) NE -1 then begin   ; PFI data?
        nimage_pfi = n_elements(data_pfi(0,0,*))
        nx = n_elements(data_pfi(*,0,0))
        ny = n_elements(data_pfi(0,*,0))
     endif $
     else begin
        nimage_pfi = 0
        if n_elements(ar_size) EQ 2 then begin
           nx = ar_size(0)
           ny = ar_size(1)
        endif $
        else if n_elements(size_ffi) EQ 1 then begin
           nx = size_ffi
           ny = size_ffi
        endif $
        else begin
           nx = 128
           ny = 128
        endelse
     endelse
     sdpfi = size(data_pfi) & tdpfi = sdpfi(sdpfi(0)+1)
     if nimage_pfi GT 0 then $
        data = make_array(nx,ny,nimage_ffi+nimage_pfi,type=tdpfi) $
     else $
        data = make_array(nx,ny,nimage_ffi+nimage_pfi, $
                          type=2)
     if pfi(0) NE -1 then begin
        index = [index_pfi]
        data(*,*,0:nimage_pfi-1) = data_pfi
        data_pfi = 0    ; save some memory
     endif 
     ccount_ffi = 0L

     if (NOT noread) then dset = mk_dset_str(infil,ss(ffi))
     for image=0,nimage_ffi-1 do begin

        ; Read one FFI image at a time to save memory
        if (NOT noread) then rd_sda, infil, dset(image), index_ffi, data_ffi $
        else begin 
           data_ffi = indata_ffi(*,*,image)
           index_ffi = inindex_ffi(image)
        endelse
        
        dffisize = size(data_ffi)
        ; decompress byte data
        if dffisize(n_elements(dffisize)-2) EQ 1 then begin
           if gt_comp(index_ffi(0)) EQ 0 then begin
              if keyword_set(sfd) then data_ffi = sfd_decomp(data_ffi) $
              else data_ffi = sxt_decomp(data_ffi) 
              ccount_ffi = ccount_ffi + 1L
           endif $
           else data_ffi = fix(data_ffi)
        endif
        if NOT keyword_set(nocorrections) then begin
           data_ffi = dark_sub(index_ffi,temporary(data_ffi))
           data_ffi = leak_sub(index_ffi,temporary(data_ffi))
        endif ;$
        ;else begin
        ;   if (pfi(0) EQ -1) and (image EQ 0) and $
        ;      (NOT keyword_set(compress)) then begin
        ;         data=make_array(nx,ny,nimage_ffi,type=dffisize(dffisize(0)+1))
        ;   endif
        ;endelse

        nxf = n_elements(data_ffi(*,0,0))
        nyf = n_elements(data_ffi(0,*,0))
        if pfi(0) EQ -1 AND image EQ 0 then index=[index_ffi] $
        else index=[index,index_ffi(0)]

        if image EQ 0 then begin
           rd_pnt,index_ffi(0),index_ffi(0),pdata,flag=flag
           if flag then message,'Could not read PNT file'
           rd_fem,index_ffi(0),index_ffi(0),femdata,/near
        endif

        if int2secarr(index_ffi(0),pdata(0)) LT 0 OR $
           int2secarr(index_ffi(0),pdata(n_elements(pdata)-1)) GT 0 $
        then begin
           rd_pnt,index_ffi(0),index_ffi(0),pdata,flag=flag
           if flag then message,'Could not read PNT file'
        endif
        if int2secarr(index_ffi(0),femdata(0)) LT 0 OR $
           int2secarr(index_ffi(0),femdata(n_elements(femdata)-1)) GT 0 $
        then begin
           rd_fem,index_ffi(0),index_ffi(0),femdata,/near
        endif
        ii = tim2dset(pdata,index_ffi(0))
;        sc = gt_hxa(pdata(ii),/sxtpix,fm=fm)  ; Sun center
;        if fm(0) then message,/info,'WARNING: HXA fiducial mark'
;        sc = hxa_suncenter(pdata(ii),code=hxacode,fem=femdata)
        sc = get_suncenter(pdata,index=index_ffi(0),code=hxacode,fem=femdata)
        sc=[sc(0,0),sc(1,0)]  ; Sun center
        if hxacode(0) LE 1 then $
           message,'WARNING: HXA Data Problem: ' + $
                   strcompress(string(hxacode(0))) + ' ' + $
                   fmt_tim(pdata(ii)),/info

        if pfi(0) EQ -1 AND image EQ 0 then begin  ; Set up reference coords
           refsc = sc
           refll = [lat,lon]
           if n_elements(ref_image) NE 1 then ref_image=0
           day = gt_day(index_ffi(0))
           time = gt_time(index_ffi(0))
           ref_day = day+time/8.64d7
           ex2int,anytim2ex(date),ms,ds
           dday = ref_day - (ds+ms/8.64d7)
           refll(1) = refll(1) + diff_rot(dday,lat)
        endif

        time = gt_time(index_ffi(0))  ; Compute differential rotation
        day = gt_day(index_ffi(0))
        dday = (day+time/8.64d7-ref_day) ;difference in days from AR reference
        lat = refll(0)
        lon = refll(1) + diff_rot(dday,lat)
        xy = hel2arcmin(lat,lon,date=fmt_tim(index_ffi(0)))*60.0/ARCSEC_PER_PIX

        ; Extract a sub-image of the FFI data
        res = 2^gt_res(index_ffi(0))
        x = xy(0)+sc(0)
        y = xy(1)+sc(1)
        x0 = fix(x)-nx/2 & x1 = fix(x)+nx/2-((nx MOD 2) EQ 0)
        y0 = fix(y)-ny/2 & y1 = fix(y)+ny/2-((ny MOD 2) EQ 0)
        xoffset = -(x-(x0+x1)/2.0)
        yoffset = -(y-(y0+y1)/2.0)

        if x0 LT 0 OR y0 LT 0 OR x1 GE nxf*res OR y1 GT nyf*res then $
           message,'Can not get a full FFI sub-image '+fmt_tim(index_ffi(0)),/info
        if x0 LT 0 then begin
           xoffset=xoffset-x0 & x1=nx-1 & x0=0
        endif
        if y0 LT 0 then begin 
           yoffset=yoffset-y0 & y1=ny-1 & y0=0 
        endif
        if x1 GE nxf*res then begin 
           xoffset=xoffset+(nxf*res-1-x1) & x0=nxf*res-nx
           x1=nxf*res-1 
        endif
        if y1 GE nyf*res then begin 
           yoffset=yoffset+(nyf*res-1-y1) & y0=nyf*res-ny
           y1=nyf*res-1 
        endif
        ; Rebin to full resolution and align to sub-pixel acc.
;        if res NE 1 then $
;           timage = rebin(data_ffi(*,*,0),nxf*res,nyf*res,/sample) $
;        else timage = data_ffi(*,*,0)
        if res NE 1 then begin
           if keyword_set(sfd) and NOT keyword_set(norm_fact) then $
              norm_fact=1
           timage = mod_res(index_ffi(0),data_ffi(*,*,0),outres=0, $
                            norm_fact=norm_fact) 
        endif $
        else timage = data_ffi(*,*,0)
        timage = timage(x0:x1,y0:y1)
        stimage = size(timage)
        ;if keyword_set(compress) and stimage(stimage(0)+1) NE 1 then $
        ;   data(*,*,nimage_pfi+image) = sxt_comp(poly_2d(temporary(timage), $
        ;                                [-xoffset,0,1,0],[-yoffset,1,0,0],1)) $
        ;else $
           data(*,*,nimage_pfi+image) = poly_2d(temporary(timage), $
                                        [-xoffset,0,1,0],[-yoffset,1,0,0],1)
     endfor
    ; if ccount_ffi GT 0 then $
    ;    message,/info,'Decompressed ' + $
    ;                  strcompress(string(ccount_ffi)) + $
    ;                  ' FFI images.'
     data_ffi = 0  ; Save some memory
     indata_ffi = 0 ; Save some memory
     time = 0.0d0
     for i=0,n_elements(index)-1 do begin
        ex2int, anytim2ex(index(i)), ms, ds
        time = [time,ds+ms/8.64d7]
     endfor
     s = sort(time(1:n_elements(time)-1))  ; Sort on the index time
     index = index(s)
     data = data(*,*,s)

  endif $
  else begin
     message,'No FFI data',/info
     data = data_pfi
     index = index_pfi
  endelse


; Do a correlation to attempt to get rid of leftover jitter.  Not too useful

  if keyword_set(korrel) then begin
      nimage = n_elements(data(0,0,*))
      nx = n_elements(data(*,0,0))
      ny = n_elements(data(0,*,0))
      xshift = fltarr(nimage)
      yshift = fltarr(nimage)
      for image=0,nimage-2 do begin

         ii = (image+1+ref_image) MOD nimage       ; First image
         ii1 = (image+ref_image) MOD nimage        ; Previous image

         xrange = [nx/4,3*nx/4]
         yrange = [ny/4,3*ny/4]
         ref = reform(float(data(xrange(0):xrange(1),yrange(0):yrange(1),ii1)))
         img = reform(float(data(xrange(0):xrange(1),yrange(0):yrange(1),ii)))

         nxk = n_elements(ref(*,0))  ; This is the size of the mosaic'ed images
         nyk = n_elements(ref(0,*))

         MAX_SHIFT = 6  ; Should be an even positive number GT 0
                        ; Actual max is MAX_SHIFT/2.0

         ; Size of subimage in korrel:
         sizes = [nxk-MAX_SHIFT,nyk-MAX_SHIFT]
         ; This centers the subimage:
         offset = fix(0.5*([nxk,nyk]-sizes))
         xy = korrel(img,ref,sizes,offset)  ; Do the correlation

         if NOT keyword_set(quiet) then $
            print,"korrel output: ",xy(0),xy(1)," image ",ii," of ",nimage

         xshift(ii) = xy(0,0)
         yshift(ii) = xy(1,0)

      endfor

      ; Now apply the cumulative korrel shifts
      xs = 0.0
      ys = 0.0
      for image=0,nimage-2 do begin
         ii = (image+1+ref_image) MOD nimage
         xs = xs + xshift(ii)  ; Compute the cumulative shift
         ys = ys + yshift(ii)

         data(*,*,ii) = poly_2d(reform(data(*,*,ii)), $
                                   [xs,0,1,0],[ys,1,0,0],1,nx,ny)

      endfor

  endif

  if keyword_set(compress) then begin
     if keyword_set(sfd) then begin
        pftype = gt_pfi_ffi(index) MOD 2
        pfi = where(pftype EQ 0)
        ffi = where(pftype EQ 1)
        if ffi(0) GE 0 then data(*,*,ffi)=sfd_comp(data(*,*,ffi))
        if pfi(0) GE 0 then data(*,*,pfi)=sxt_comp(data(*,*,pfi))
     endif $
     else data = sxt_comp(data)
  endif

  if noread then begin
     infil = temporary(index)
     ss = temporary(data)
  endif

en