pro son3view,xdim,ydim,tdim,filter_mask
;
;Calculates the subsonic filter in 3 dimensions (k_x,k_y,w)
; for tuning and visualization
;
;
;input parameters: dimensions of the array, and...
;
;  time_step: mean time distance of the images
;      v_ph : maximum phase velocity
;
;    output : filtered images
;
;  History: 
;	   Based on the J.Hirzsberger's procedure "wave.pro"
;	   Adaptation made on 6,Nov,1995  by M.Sobotka and J.A.Bonet
;


;PARAMETERS TO BE EDITED-------------------------------------------------
;
      scale=0.167	;0.062  ;arcsec/pixel

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

;
;  reading input parameters
;
      if xdim mod 2 ne 0 then xdim=xdim-1
      if ydim mod 2 ne 0 then ydim=ydim-1
      print,''
      read,'Mean time distance : ', time_step
      read,'Max. Phase Velocity [km/s] : ', v_ph
      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
      if cut ne 0 then begin
	read,'length of the damped edge (% of total length)=',perc
      endif

;      
;  number of images should be even
;
      if tdim mod 2 ne 0 then tdim=tdim-1
;
;  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)
     if cut ne 0 then smooth_t=fix(tdim*perc/100.)
 				; width of the edge in t-dimension

;  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.

fin: 
     end