pro vmake

   psf = iudf_psf('IUDF_CH1_v0.9_psf.fits.gz', 53.12, -27.81, spsf=spsf, map=map)

   tv,bytscl(psf,min=-5e-5,max=5e-5)

end

; reconstruct effective IRAC psf at the specified ra, dec (decimal degrees)
function iudf_psf, fpsf, ra, dec, spsf=spsf, map=map

    ; multi-extension fits file: first extension is the super psf
    ; second extension is the weight and PA map
    ; 4d map = [nx,ny,2,naor+1]
    ; first two dimensions specify a rectangular grid containing 2 vectors on each location
    ; - a weight vector and rotation angle vector, with one entry per aor
    ; - the first element of each vector contains the ra and dec on that grid position
    if not keyword_set(spsf) then spsf = readfits(fpsf, hpsf, exten=1,/silent)
    if not keyword_set(map) then map = readfits(fpsf, exten=2,/silent)

    szpsf = size(spsf,/dim)
    sz = size(map,/dim)
    len = sz[3]-1
    xy0 = (szpsf-1)/2.
    rscl=1.0

    cra = map[*,*,0,0]
    cdec = map[*,*,1,0]
    ix = value_locate(cra[*,0], ra)   ; could interpolate for finer grid result
    iy = value_locate(cdec[0,*], dec)
    psf = fltarr(size(spsf,/dim))

    ; if all weights zero then return
    if total(map[ix,iy,0,1:*]) eq 0 then return, 0.0

    for i=0,len-1 do begin
        aor_w = map[ix,iy,0,i+1]
        aor_pa = map[ix,iy,1,i+1]
        if aor_w[0] le 0 then continue  ; skip if no weight
        psf += rot(spsf,-aor_pa[0],rscl[0],xy0[0],xy0[1],cubic=-0.4,missing=0.0)*aor_w[0]/rscl^2
    end
    totcor = 1.0/total(psf)

    return, psf
end