! @(#)convolfft.prg	17.1.1.1 (ESO-DMD) 01/25/02 17:12:21
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Midas procedure convolfft.prg
!  to convolve or correlate an object image with a PSF via FFT
!
!  If the psf frame is smaller than the object frame, the PSF is inserted
!  into a frame of same size as object frame (centered) and the FFT/IMAGE done 
!  on these frames, compute the convolution real and imaginary part,
!  transform back with "FFT/INV" and swap the 4 quadrants of the result frame.
! 
!  If object + psf frame have same size, we do the FFT stuff on these frames
!  directly.
!
!  use via   @a convolfft [what] [object] [psf] [result] [del_flag]
!
!  what           CONVOLVE or CORRELATE, default = CONVOLVE
!  object         object frame, defaulted to object
!  psf            psf frame (psf centered), defaulted to psf
!  result         name of result image, defaulted to result
!  del_flag       Yes or No for temporary file deletion, defaulted to Yes
!                 Temporary files are named: mid_xp, mid_xo, mid_xr,
!		  mid_cr, mid_ci, mid_for, mid_foi, mid_fpr, mid_fpi
!
!  Note:
! Data dimensions of object frame have to be powers of 2.
! Data dimensions of psf have to be <= dimensions of object.
! If key MID$SPEC(1:5) = DEBUG, the imaginary part of inverse FFT (should be
! practically zero, for checking only) is saved in frame middummk.imagi
!
!  M. Rosa   ST-ECF	851020
!  A. Rousset  Obs. de Lyon  920129
!  K. Banse  ESO-IPG	860915, 901212, 911119, 921016, 940517, 000911
!
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
define/param p1 convolve C "Enter action, convolve or correlate: "
define/param p2 object IMA "Enter object frame: "
define/param p3 psf IMA "Enter PSF: "
define/param p4 result IMA "Enter name for result frame: "
define/param p5 Y c "Enter delete_flag: "
! 
if p1(1:3) .eq. "COR" then
   define/local ops/c/1/2 "+-"
elseif p1(1:3) .eq. "CON" then
   define/local ops/c/1/2 "-+"
else
   write/error 5
   return
endif
define/local infra/c/1/60 {p2}
define/local outfra/c/1/60 {p4}
if error(2) .lt. 2 then
   define/local verbose/c/1/1 y
else
   define/local verbose/c/1/1 n
endif
!
define/local nl/i/1/2 {{p2},naxis},0
if nl(1) .eq. 1 then
   define/local np/i/1/2 {{p2},npix(1)},1
   define/local nq/i/1/2 {{p3},npix(1)},1
   define/local sta/d/1/1 {{p2},start(1)}
   define/local ste/d/1/1 {{p2},step(1)}
   define/local maxis/i/1/1 1
else
   define/local np/i/1/2 {{p2},npix(1)},{{p2},npix(2)}
   define/local nq/i/1/2 {{p3},npix(1)},{{p3},npix(2)}
   define/local sta/d/1/2 {{p2},start(1)},{{p2},start(2)}
   define/local ste/d/1/2 {{p2},step(1)},{{p2},step(2)}
   define/local maxis/i/1/1 2
endif
define/local rmin/r/1/1 {{p3},lhcuts(3)}	!minimum of psf frame
! 
if np(1) .lt. nq(1) goto wrong_dims
if np(2) .lt. nq(2) goto wrong_dims
! 
! do FFT on object frame
if verbose .eq. "Y" write/out "Doing FFT on object frame {p2} ..."
fft/image {infra} ? mid_for mid_foi
! 
if np(1) .eq. nq(1) then
   if np(2) .eq. nq(2) then
!					do the FFT directly on psf
      if verbose .eq. "Y" write/out "Doing FFT on psf frame {p3} ..."
      fft/image {p3} ? mid_fpr mid_fpi
      goto next_step
   endif
endif
!				insert psf into big frame and do the FFT there
 nl(1) = (np(1)-nq(1))/2+1			!compute offset for insertion
 nl(2) = (np(2)-nq(2))/2+1
if np(2) .eq. 1 then
   create/image mid_xp 1,{np(1)} ? poly {rmin}
   insert/image {p3} mid_xp @{nl(1)}
else
   create/image mid_xp 2,{np(1)},{np(2)} ? poly {rmin}
   insert/image {p3} mid_xp @{nl(1)},@{nl(2)}
endif
if verbose .eq. "Y" write/out "Doing FFT on expanded psf frame mid_xp ..."
fft/image mid_xp ? mid_fpr mid_fpi
if p5(1:1) .ne. "N" delete/image mid_xp NO
!
!  multiply Fourier transforms
next_step:
if verbose .eq. "Y" write/out "Multiplying Fourier transforms ..."
compute/pixel mid_cr = (mid_for*mid_fpr){OPS(1:1)}(mid_foi*mid_fpi)
compute/pixel mid_ci = (mid_foi*mid_fpr){OPS(2:2)}(mid_for*mid_fpi)
! 
!  do inverse Fourier transform
if verbose .eq. "Y" write/out "Doing inverse FFT ..."
fft/inverse mid_cr mid_ci {outfra} middummk.imagi
if verbose .eq. "Y" write/out "Swapping result frame {outfra} ..."
nq(1) = np(1)/2-1
nq(2) = np(2)/2-1
@a swap4 {outfra} {outfra} ? @{nq(1)},@{nq(2)}
if verbose .eq. "Y" write/out "Rescale result frame {outfra} ..."
compute/image {outfra} = {outfra}*({np(1)}*{np(2)})
! 
if p5(1:1) .ne. "N" then
   delete/image mid_for NO
   delete/image mid_foi NO
   delete/image mid_fpr NO
   delete/image mid_fpi NO
   delete/image mid_cr NO
   delete/image mid_ci NO
endif
! 
if mid$spec(1:5) .eq. "DEBUG" then
   @a swap4 middummk.imagi middummk.imagi
   compute/image middummk.imagi = middummk.imagi*{np(1)}
endif
! 
if maxis .eq. 1 then
   write/descr {outfra} start {sta(1)}
   write/descr {outfra} step {ste(1)}
else
   write/descr {outfra} start {sta(1)},{sta(2)}
   write/descr {outfra} step {ste(1)},{ste(2)}
endif
write/descr {outfra} history/c/-1/80 "@a convolfft {p1} {p2} {p3} {p4} {p5} "
return
! 
wrong_dims:
write/out "PSF image has to be smaller than input frame..."