#<
# Cfft
#
# DESCRIPTION
# Complex fast fourier transform
#
# USAGE
# Cfft < in.h > out.h
#
# INPUT PARAMETERS
# sign      - integer     forward transform
#                         -1 : reverse transform
#                         ( FDGP sign convention )
#  n1,n2,n3  - integer   input cube dimensions ( complex data ).
#                        (n1 < 2048)
#                        (n1 = 2*(ntime1/2+1) if in frequency domain)
#  esize     - integer   [8]
#
# CATEGORY
# Seis:Filter
#
# COMPILE LEVEL
# DISTR
#>
#%
# Author - Peter Mora
#  modified - Ray Abma
#             increased buffer sizes to match old version.
#
# ----------
#  Keyword:  Fourier-transform
# ----------
   integer input,output,fetch,infd,outfd,ier,esize
   integer sreed, srite,maxnt
	
   complex data(16384)

		maxnt=2048
   call initpar()
   call doc(source)

   iop=1			;ier=fetch('op sign','i',iop)
   if(fetch('n1 nt','i',nt)==0)	call erexit('n1 is required')
   iop=-iop
   if(nt>maxnt)			call erexit('n1 too big')

   nx=1				;ier=fetch('n2 nx','i',nx)
   np=1				;ier=fetch('n3 np','i',np)
   esize=8	                ;ier=fetch('esize','i',esize)

   call putch('n1','i',nt)	;call putch('n2','i',nx) ;call putch('n3','i',np)

   call hclose()

#  loop over planes

   do ip=1,np {
      do ix=1,nx {
         ier=sreed('in',data,esize*nt)
         call sepcfft(data,nt,iop)
         ier=srite('out',data,esize*nt)
      }
   }

 end


subroutine sepcfft(x,n,iop)

#	complex one dimensional fast fourier transform
integer n
complex x(n),y(16384),w,t
#
save y
#
integer i,j,k,k0,l,m,n2,inc,iop,last
real pi,scale
pi=3.1415965

last=0

#	compute trigonometric values

if(last!=n) { last=n ; call cexps(y,n) }

#	compute fast Fourier transform

n2=n/2
j=1
do i=1,n {
   if(i<=j) { t=x(j) ; x(j)=x(i) ; x(i)=t }
   m=n2
   while (j>m&&m>=1) { j=j-m ; m=m/2 }
   j=j+m }
l=1
while(l<n) {
   inc=l+l ; k0=n/inc ; k=1
   do m=1,l {
      w=cmplx(real(y(k)),-iop*aimag(y(k))) ; k=k+k0
      do i=m,n,inc { t=w*x(i+l) ; x(i+l)=x(i)-t ; x(i)=x(i)+t }}
   l=inc }
if(iop<0) { scale=1./n ; do i=1,n { x(i)=x(i)*scale } }
return
end


subroutine cexps(y,n)

#	compute complex exponentials

integer last,i,j,n2,n
complex y(n),z(16384)
last=0
if(n==last) {
   do i=1,n { y(i)=z(i) }}
else {
   n2=n/2 ; last=n
   do i=1,n/4+1 { y(i)=cmplx(cos((i-1)*3.14159265*2/n),0.)	; z(i)=y(i) }
   do i=1,n/4+1 { y(i)=cmplx(real(y(i)),real(y(n/4+2-i)))	; z(i)=y(i) }
   do i=2,n/4 { j=n2-i+2 ; y(j)=cmplx(-real(y(i)),aimag(y(i)))	; z(j)=y(j) }}
return
end