pro sxt_interp,index,data,ntimes,unc_data,satpix,verbose=verbose
; --------------------------------
; Setup/initialization
; --------------------------------
sz = size(data)
if (sz(sz(0)+1) eq 1) or (n_elements(data(0,0,*)) lt 2) or (sz(0) ne 3) then begin
print,'**** Warning in sxt_interp: ***' & tbeep
if sz(sz(0)+1) eq 1 then print,' Data must be decompressed'
if n_elements(data(0,0,*)) lt 2 then $
print,' Require at least 2 images for interpolation'
if (sz(0) ne 3) then print,' Data must be a 3-d array'
help,data
return
endif
nntimes = anytim2ints(ntimes) ; Convert to structure format
tim0 = int2secarr(index) ; Get array of times
tim1 = int2secarr(nntimes, index(0)) ; Times when images are desired
j0 = sort(tim0) ; Time sorted
tim0 = tim0(j0)
nimg2 = n_elements(tim1) ; Number of new images
N_x = n_elements(data(*,0,0)) ; Number of x pixels
N_y = n_elements(data(0,*,0)) ; Number of y pixels
sz(3) = nimg2
new_data = make_array(size=sz) ; Set up output array
; Round off if the final data is not floating:
if sz(sz(0)+1) le 3 then round_off0 = 0.5 else round_off0 = 0.
; Set up new satpix and unc_data arrays
if n_params() ge 4 then begin
sz1 = size(unc_data) & sz1(3) = nimg2
unc_data2 = make_array(size=sz1)
; Round off if the unc_data is not floating:
if sz1(sz1(0)+1) le 3 then round_off1 = 0.5 else round_off1 = 0.
endif
if n_params() ge 5 then begin
sz2 = size(satpix) & sz2(3) = nimg2
satpix2 = make_array(size=sz2)
sz2_1 = sz2(0:4) ; Set up for satpix2
sz2_1(0) = 2 & sz2_1(3) = sz2_1(4) ; Preserve the type
endif
nk = intarr(2,nimg2) ; Indices of arrays which will
;- be interpolated
ss = intarr(nimg2) ; Index of the shorter index
; --------------------------------
; Find image pairs to interpolated
; --------------------------------
for i=0,nimg2-1 do begin
jj = where(tim0 ge tim1(i),nj)
if (nj gt 1) and (nj ne n_elements(tim0)) then begin
nk(0,i) = j0(min(jj)-1)
nk(1,i) = j0(min(jj))
endif else if (nj eq n_elements(tim0)) then begin
nk(*,i) = j0([0,1]) ; Take the first two images
print,'** sxt_interp warning: Extrapolating for time = ',fmt_tim(nntimes(i)),' (',strtrim(i,2),')'
endif else begin ; nj = 0
nk(1,i) = j0(n_elements(j0)-1) ; Take the last two images
nk(0,i) = j0(n_elements(j0)-2)
if tim1(i) gt max(tim0) then $
print,'** sxt_interp warning: Extrapolating for time = ',fmt_tim(nntimes(i)),' (',strtrim(i,2),')'
endelse
; Choose the shorter exposure for the exposure normalization:
texp = gt_expdur(index(nk(*,i))) ; Get the exposure durations
xx = min(texp,xmin) & ss(i) = nk(xmin,i) ; Choose the shorter exposure to normalize
; Make sure the images to be interpolated have the same resolution and filters
res = gt_res( index(nk(*,i))) ; Check the resolutions
filb = gt_filtb(index(nk(*,i))) ; Check the filters
if (res(0) ne res(1)) or (filb(0) ne filb(1)) then begin
print,'** sxt_interp warning: Can not interpolate for time = ', $
fmt_tim(nntimes(i))
if res(0) ne res(1) then print,'Resolutions of the adjacent images are different',format='(15x,a)'
if filb(0) ne filb(1) then print,'Filters (B) of the adjacent images are different',format='(15x,a)'
nk(1-xmin,i) = nk(xmin,i)
endif
endfor ; Find image pairs to interpolate
new_index = index(ss) ; Create the new index array
; --------------------------------
; Call interp_arr to do the interpolation
; --------------------------------
dt = gt_expdur(index) ;jmm, 29-sep-95
new_dt = gt_expdur(new_index) ;jmm, 29-sep-95
for i = 0, nimg2-1 do begin ; Do the interpolation
if keyword_set(verbose) then $
print, fmt_tim(index(nk(0, i))), '\\',fmt_tim(nntimes(i)),'\\', $
fmt_tim(index(nk(1, i))), '\\'
; New uncertainty calculation used in Interp_arr, Jmm, 29-sep-95
img1 = data(*, *, nk(0, i))/dt(nk(0, i))
img2 = data(*, *, nk(1, i))/dt(nk(1, i))
if n_params() ge 4 then begin ;need uncertainties
uimg1 = unc_data(*, *, nk(0, i))/dt(nk(0, i))
uimg2 = unc_data(*, *, nk(1, i))/dt(nk(1, i))
ddata = new_dt(i) * interp_arr(img1, img2, index(nk(0, i)), index(nk(1, i)), $
nntimes(i), uimg1, uimg2, udata2)
unc_data2(*, *, i) = new_dt(i)*udata2+round_off1
endif else ddata = new_dt(i) * interp_arr(img1, img2, index(nk(0, i)), $
index(nk(1, i)), nntimes(i))
; End of new stuff, Jmm, 29-sep-95
ex2int, anytim2ex(nntimes(i)), msod, ds79
new_index(i).gen.time = msod
new_index(i).gen.day = ds79
; Assign to the new array (Could go negative if extrapolated)
new_data(0, 0, i) = (ddata + round_off0) > 0
; Make any pixel that is saturated in either image saturated in the final image
if n_params() ge 5 then begin
jj = where((satpix(*, *, nk(0, i)) gt 0) or (satpix(*, *, nk(1, i)) GT 0), nj)
if nj gt 0 then begin
ddata = make_array(size = sz2_1)
ddata(jj) = 1b ; saturated pixels
satpix2(0, 0, i) = ddata
endif
endif ; n_params() ge 5 (satpix)
endfor
; --------------------------------
; Copy new_index and new_data to index and data: (Clobber input!)
; --------------------------------
data = 0 ; clear memory
data = new_data
new_data = 0 ; clear memory
index = new_index
if n_params() ge 4 then begin
unc_data = unc_data2
unc_data2 = 0 ; clear memory
endif
if n_params() ge 5 then begin
satpix = satpix2 ; Saturated pixels
satpix2 = 0 ; clear memory
endif
end