pro sonic3d,name,xdim,ydim,tdim
;
; Subsonic filter in 3 dimensions (k_x,k_y,w)
; The whole datacube is filtered at once, so that the computer
; memory must be approx. 10x larger than the datacube size.
;
; Input parameters:
;
;      name : I/O file name (string)
;      xdim : x-dimension of the input images
;      ydim : y-dimension of the input images
;      tdim : t-dimension of the input images
;  time_step: mean time distance of the images
;      v_ph : maximum phase velocity
;
;  I/O file : unformatted binary file containing the 3D x,y,t
;	 INTEGER array (datacube). Input overwritten by output!!
;	 Recommended normalisation: Iph = 10000
;
; CAUTION!   number of images and frame dimensions MUST BE EVEN
;		(not checked)
;
; History: 
;	   Based on the J.Hirzberger's procedure "wave.pro"
;	   Adaptation made on 6,Nov,1995 by M.Sobotka and J.A.Bonet.
;	   I/O changes and 'temporary' command use made on May 21 1997
;	   by M.S.


;PARAMETERS TO BE EDITED-------------------------------------------------
;

      scale=0.070  ;arcsec/pixel

;
;-------------------------------------------------------------------------

;
;  reading input parameters
;
      print,''
      read,'Mean time distance [s] : ', time_step
      read,'Max. Phase Velocity [km/s] : ', v_ph
      read,'Shall we apodize?: YES(=1), NO(=0): ',ap
      print,'Filter cutoff:' 
      read,'abrupt (=0), smoothed WEDGE-MODE(=1), CONSTANT-MODE(=2) ',cut
      if (cut gt 2) or (cut lt 0) then begin
	print,'UNEXPECTED VALUE. Program terminated.' 
	goto,fin
      endif				;last chance to escape...
      if ap ne 0 or cut ne 0 then begin
	read,'length of the damped edge (% of total length)=',perc
      endif

;
;  reading input 3D array imeq - file organized x,y,t, but
;				 imeq organized t,x,y !
     imeq=fltarr(tdim,xdim,ydim)
     openr,1,name	
     as=assoc(1,intarr(xdim,ydim))	; input is integer !
     for i=0,tdim-1 do imeq(i,*,*)=float(as(i))	
     close,1
     as=0
;
;  defining unit in the Fourier domain
;
     print,'Spatial Resolution =',scale,' arcsec/pixel'
     print,''
     kx_step=1./(scale*725.*xdim)
     ky_step=1./(scale*725.*ydim)
     w_step=1./(time_step*tdim)
;
;  optional apodization with a cosine window
;
     if ap ne 0 or cut ne 0 then begin
       smooth_x=fix(xdim*perc/100.) ; width of the edge in x-dimension
       smooth_y=fix(ydim*perc/100.) ; width of the edge in y-dimension
       smooth_t=fix(tdim*perc/100.) ; width of the edge in t-dimension
     endif

     if ap ne 0 then begin
       av=mean(imeq)
       imeq=imeq-av
       xmask=fltarr(xdim)+1.
       ymask=fltarr(ydim)+1.
       tmask=fltarr(tdim)+1.
       print,'Now smoothing edges.'
;    smoothing with a cosine
       for i=0,smooth_x do xmask(i)=(1.-cos(!pi*float(i)/float(smooth_x)))/2.
       for i=0,smooth_y do ymask(i)=(1.-cos(!pi*float(i)/float(smooth_x)))/2.
       for i=0,smooth_t do tmask(i)=(1.-cos(!pi*float(i)/float(smooth_t)))/2.
       xmask(xdim-smooth_x:xdim-1)=reverse(xmask(1:smooth_x))
       ymask(ydim-smooth_y:ydim-1)=reverse(ymask(1:smooth_y))
       tmask(tdim-smooth_t:tdim-1)=reverse(tmask(1:smooth_t))
       for i=0,tdim-1 do imeq(i,*,*)=imeq(i,*,*)*tmask(i)
       for i=0,xdim-1 do imeq(*,i,*)=imeq(*,i,*)*xmask(i)
       for i=0,ydim-1 do imeq(*,*,i)=imeq(*,*,i)*ymask(i)    
       imeq=imeq+av
     endif
;
;  apply FFT
;
     print,'Now do FFT.'
     imeq=fft(imeq,-1,/overwrite)
;
;  prepares the filter
;
     nx=xdim/2+1
     ny=ydim/2+1
     nt=tdim/2+1
     filter_mask=fltarr(tdim,xdim,ydim)
;
;  calculate subsonic filter with abrupt cutoff;
;

   if cut eq 0 or cut eq 1 then begin
     print,'Now calculating filter.'
     for j=0,ny-1 do begin
       for i=0,nx-1 do begin
         k_by_v=sqrt((i*kx_step)^2+(j*ky_step)^2)*v_ph
         for k=1,nt-1 do if k*w_step le k_by_v then filter_mask(k,i,j)=1.
       endfor
     endfor
   endif

;optional smoothing of the cutoff in the filter:  WEDGE-MODE transition

   if cut eq 1 then begin      ;optional smoothing of the cutoff. Wedge-mode
     five=nint(50./perc)
     smoo=fltarr(smooth_t,smooth_t)
     for i=0,smooth_t-1 do begin
       for j=0,smooth_t-i-1 do begin
         smoo(i,j)=(1.+cos(!pi*float(j)/float(smooth_t-i)))/2.
       endfor
     endfor
     for i=0,nx-1 do begin
       for j=0,ny-1 do begin
         k=0
         repeat begin
           if filter_mask(k+1,i,j) ne 1. then begin 
           if k+1 ge five then trans=(k+1)/five else trans=1
           if k ge 1 then filter_mask(k-trans+1:k-trans+smooth_t,i,j)=$
                                     smoo(smooth_t-trans,*)
             k=nt-1
           endif else begin
             k=k+1
           endelse
         endrep until k eq nt-1
       endfor
     endfor
   endif

;  subsonic filter with optional smoothing of the cutoff in the filter: 
;  CONSTANT-WIDTH transition

   if cut eq 2 then begin 
     print,'Now calculating filter.'
     ntrans=smooth_t-1
     for j=0,ny-1 do begin
       for i=0,nx-1 do begin
         k_by_v=sqrt((i*kx_step)^2+(j*ky_step)^2)*v_ph
         for k=1,nt-1 do begin
           dif=k*w_step-k_by_v
           if dif le 0 then begin
	     filter_mask(k,i,j)=1.
	   endif else begin
	     dif=dif/w_step	     
	     in=0
	     for kk=k,k+ntrans do begin
	       if kk le nt-1 then begin
	         filter_mask(kk,i,j)=0.5+0.5*cos(!pi*(dif+in)/smooth_t)
	         in=in+1
	       endif else begin
	         goto,label
	       endelse
	     endfor
	     goto,label
	   endelse
         endfor
         label:
       endfor
     endfor
   endif
;
;
;  calculate the other 7 octants of the filter
;
     filter_mask(*,*,ny:ydim-1)=reverse(filter_mask(*,*,1:ny-2),3)
     filter_mask(*,nx:xdim-1,*)=reverse(filter_mask(*,1:nx-2,*),2)
     filter_mask(nt:tdim-1,*,*)=reverse(filter_mask(1:nt-2,*,*),1)
     filter_mask(0,*,*)=1.
;
;  multiply transformed image by the filter
;
     imeq=temporary(imeq)*filter_mask
     filter_mask=0	;to save space
;
;  apply the IFFT
;
     print,'Now do Inverse FFT:'
     imeq=FFT(imeq,1,/overwrite)
     imeq=nint(float(imeq))	;taking real part and converting to integer
;
;  saving filtered 3D array
;
     print,'Saving the result into '+name
     openw,1,name
     for i=0,tdim-1 do writeu,1,imeq(i,*,*)
     close,1
;
fin: 
     end