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