PRO big3dfft,pm,first,last,dimx,dimy,bxdim,bydim,xdim,ydim,x_inf,y_inf
;
;PURPOSE:
;        Computes by parts the fft of a large 3D-array I(x,y,t) which cannot
;        be handled in one run by the CPU/RAM.
;
;METHOD: 
;       Let us assume that we have a series of n=first-last+1 images of 
;	dimensions dimx*dimy, and we are interested in the fft of a sub-array 
;	of dimensions xdim*ydim*n. The first step consist in computing the
;       fft of separated sub-images(xdim*ydim), thus creating in the hard
;       disk a series of n complex Fourier-transformed-sub-images.
;       From this series, 3D-complex-sub-arrays of dimensions bxdim*bydim*n,
;       able to enter in the cpu in one run, are created, and then for each
;	spatial position is extracted a temporary array and computed its 1D-fft.
;	The procedure is based on the splitting of the 3D Fourier integral
;	as follows:
;	S{exp(-2*pi*j*nt*t)[SS{f(x,y,t)exp(-2*pi*j*(nx*x+ny*y))dxdx}]dt}
;
;	WARNING!  bxdim & bydim should be as larger as possible to minimize
;	the intermediate read-write processes in the execution of the program.
;
;HISTORY:
;	Edited by J.A.Bonet on Jul,7,1995. Modified on Nov,17,1995,
;	as large3d_fft.pro.
;	I/O modifications by S.Stangl and M.Sobotka,
;	tuned Apr,18,2000.
;  
;
;INPUTS:
;      data files with filenames 0, 1, 2, ...
;      pm : direction of the FFT: direct (=-1), inverse (=1)
;      first: number of the first image of the series
;      last : number of the last image of the series
;      dimx : x-dimension of the original images of the series
;      dimy : y-dimension of the original images of the series
;      bxdim: x-dimension of the sub-box in the Fourier domain partition
;      bydim: y-dimension of the sub-box in the Fourier domain partition
;
;    OPTIONAL (if not set, the default is
;	       xdim=dimx, ydim=dimy, x_inf=0, y_inf=0):
;      xdim : x-dimension of the box extracted in the orig.images (even number)
;      ydim : y-dimension of the box extracted in the orig.images (even number)
;      x_inf: lower left corner of the box (x-dimension)
;      y_inf: lower left corner of the box (y-dimension)
;
;OUTPUTS:
;      The resulting Fourier Transform consist on a series of last-first+1
;      separated complex images, directly stored on the hard disk, with
;      names fft.n where n=0,...,first-last.
;      Subsequent filtering processes requires sequencial manipulation (not
;      in one go) of these complex images.
;      The inverse Fourier Transform will also produce a complex
;      formatted array, stored in files F.n
;____________________________________________________________________________
;
;  defining paths         TO BE CHANGED BY EDITOR
;
      str1=string('/scratch/work/')  ; path
;
;  it is assumed that the filenames are 0, 1, 2, 3 etc.
;____________________________________________________________________________

if n_params() lt 8 then begin
	xdim=dimx
	ydim=dimy
	x_inf=0
	y_inf=0
endif
if xdim mod 2 ne 0 then xdim=xdim-1
if ydim mod 2 ne 0 then ydim=ydim-1

selection:
case pm of
    -1: begin
	  strR=str1
	  strW=str1+'fft.'
	  print,'direct FFT'
	end
     1: begin
	  strR=str1+'fft.'
	  strW=str1+'F.'
	  print,'inverse FFT'
	end
  else: begin
	   read, 'valid pm values: -1 (forward), 1 (inverse) : ',pm
	   print, pm 
	   goto, selection
	end
endcase
;
;
;..............................................................................
;
;  reads separate images and computes their the 2D-Fourier Transform.
;

     for i=first,last do begin
       dcn=strtrim(i,1)
       if pm eq -1 then im=fltarr(dimx,dimy) else im=complexarr(xdim,ydim)
       print, 'reading image : '+ strR + dcn
       openr,1,strR + dcn
         readu,1,im
       close,1

       fim=fft(im(x_inf:x_inf+xdim-1,y_inf:y_inf+ydim-1),pm,/overwrite)
       print, '.... and writing its 2D-Fourier Transform ' + strW + dcn
       openw,1,strW + dcn
	 writeu,1,fim
       close,1
     endfor
     fim=0
     im=0
;
;creates 3d_subarrays (nu_x,nu_y,t) with dimensions in spatial freq. small
;enough to enter the 3d_subarray in the cpu, and computes the fft of 
;1D temporal arrays correponding to separated x,y pixels
;
     a=xdim mod bxdim
     b=ydim mod bydim
     n=fix(xdim/bxdim)
     m=fix(ydim/bydim)
     if a ne 0 then begin
       dix=intarr(n+1)+bxdim
       dix(n)=a
     endif else begin
       dix=intarr(n)+bxdim
     endelse
     if b ne 0 then begin
       diy=intarr(m+1)+bydim
       diy(m)=b
     endif else begin
       diy=intarr(m)+bydim
     endelse
     nelx=fix(n_elements(dix))
     nely=fix(n_elements(diy))
     print, 'Number of sub-arrays = ',strtrim(nelx*nely,1)

     fim=complexarr(xdim, ydim)		; loop over sub-boxes
     for jbox=0,nely-1 do begin
       for ibox=0,nelx-1 do begin
         box3d=complexarr(dix(ibox),diy(jbox),last-first+1)
         print, 'Now reading in the Fourier domain the sub-array: ',ibox,jbox

         for i=first,last do begin
	   dcn=strtrim(i,1)
           openr,1,strW + dcn
             readu,1,fim
           close,1
           box3d(*,*,i-first)=fim(ibox*bxdim:ibox*bxdim+dix(ibox)-1,$
                            	  jbox*bydim:jbox*bydim+diy(jbox)-1)
         endfor

         for y=0,diy(jbox)-1 do begin	; fft of time-columns
	   for x=0,dix(ibox)-1 do begin
	     tt=box3d(x,y,*)
             box3d(x,y,*)=fft(tt,pm,/overwrite)
	   endfor
	 endfor 

         print,'....and writing its t-Fourier transform ' + strW
         for i=first,last do begin
	   dcn=strtrim(i,1)
	   openr,1,strW + dcn
             readu,1,fim
           close,1
           fim(ibox*bxdim:ibox*bxdim+dix(ibox)-1,$
	       jbox*bydim:jbox*bydim+diy(jbox)-1)=box3d(*,*,i-first)
	   openw,1,strW + dcn
             writeu,1,fim
           close,1
         endfor

     endfor
     endfor

     print, ''
     print, '*** END OF FFT ***'

END