; must .run coalign first

 function shift_align,a,b,inshift=inshift,outshift=outshift,dmin=dmin,dmax=dmax

; shifts b to match a: uses coalign to get shift, grids up and
; back down to get shifted image
; inshift is alternative shift supplied
; outshift returns derived shift
; dmin, dmax are limits on range to use for correlation

 asz=size(a)
 bsz=size(b)
 if ((asz[0] ne bsz[0]) or (asz[1] ne bsz[1]) or (asz[2] ne bsz[2])) $
                then begin
    print,'Dimensions must be the same.'
    return,-1
 endif

; sh is the shift to be given to b to match a
 if keyword_set(inshift) then sh=inshift else begin
    amin=a & bmin=b
    if keyword_set(dmin) then begin
       amin=a>dmin
       bmin=b>dmin
    endif
    amax=amin & bmax=bmin
    if keyword_set(dmax) then begin
       amax=amin<dmax
       bmax=bmin<dmax
    endif
    sh=coalign(amax,bmax)
 endelse
 print,'Shift in pixels is ',sh
 if keyword_set(outshift) then outshift=sh

; OLD VERSION USING REGRIDDING - superceded by shift_image
;; regrid to a factor of 4
;
;; expand array by nex pixels
; nex=128
;
;; bgrid=congrid(b,4*asz[1],4*asz[2],cubic=-0.5)
; bgrid=rebin(b,4*asz[1],4*asz[2])
; bbig=0.0*fltarr(nex+4*asz[1],nex+4*asz[2])
; bbig[nex/4:nex/4-1+4*asz[1],nex/4:nex/4-1+4*asz[2]]=bgrid
;
; if (sh[0] ge 0.0) then xsh=fix(4*sh[0]+0.4999) else $
;                        xsh=-fix(-4*sh[0]+0.4999)
; if (sh[1] ge 0.0) then ysh=fix(4*sh[1]+0.4999) else $
;                        ysh=-fix(-4*sh[1]+0.4999)
; print,'Shifting ',xsh,ysh,' quarter pixels'
;
;; have to set /SAMPLE on binning back down or else you don't get
;; what you started with when there is no shift
;
; shift_im=rebin(bbig[nex/4-xsh:nex/4-1-xsh+4*asz[1],nex/4-ysh:nex/4-1-ysh+4*asz[2]],$
;                asz[1],asz[2],/sample)

 shift_image,b,shift_im,sh
 
 return,shift_im
 end