c  Copyright (C) Stefan Goedecker, CEA Grenoble, 2002
c  This file is distributed under the terms of the
c  GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .

c The theoretical background is described in :
c 1) S. Goedecker: Rotating a three-dimensional array in optimal
C positions for vector processing: Case study for a three-dimensional Fast
C Fourier Transform, Comp. Phys. Commun. \underline{76}, 294 (1993)
C 2) S. Goedecker: Fast radix 2,3,4 and 5 kernels for Fast Fourier
C Transformations on computers with overlapping multiply-add instructions,
C SIAM Journal on Scientific Computing \underline{18}, 1605 (1997)
C Citing these 2 references in any scientific publication using
C these FFT routines is greatly appreciated.

c  4-dim generic (not processor optimized) FFT:
  
c  HEADER PART: ntime=1: Checks accuracy against NAG routine
c               ntime<>1: Measures speed

	Implicit real*8 (a-h,o-z)

	parameter(ntime=1)
C some reasonable values of ncache: IBM/RS6000/590: 16*1024 ; IBM/RS6000/390: 3*1024 ; 
c       IBM/PwPC: 1*1024 ; SGI/MIPS/R8000: 16*1024 ; DEC/Alpha/EV5 8*1024
c       ncache=0 on all vector machines
c       But if you care about performance find the optimal value of ncache yourself!
	parameter(ncache=16*1024)
C   dimension parameters
	parameter(n1=3, n2=5, n3=6, n4=128)   
C   parameters for FFT
        parameter(nd1=n1+1, nd2=n2+1, nd3=n3+1, nd4=n4+1)
C   parameters for NAG
	parameter(lwork=3*(n1+n2+n3+n4))
c   timing array on DEC
	dimension tarray(2)

C general array
	dimension zin(2,n1*n2*n3*n4)
C arrays for FFT (add 1 to size of zw so no F90 standard is required if ncache=0)
	dimension z(2,nd1*nd2*nd3*nd4,2),zw(ncache+1)
C arrays for NAG
	dimension xr(n1*n2*n3*n4),xi(n1*n2*n3*n4),nd(4),work(lwork)
	nd(1)=n1
	nd(2)=n2
	nd(3)=n3
	nd(4)=n4


	Print*,'                    '
	Print*,'                    '
66	format(a,4(i4),3x,i3)
	write(6,66) 'FFT TEST for n1 ,n2 ,n3, ,n4, ncache=',
     1               n1,n2,n3,n4,ncache/1024
	Print*,'                    '
	isign=1
1234	call init(n1,n2,n3,n4,nd1,nd2,nd3,nd4,zin,z,xr,xi)

	inzee=1
	call fft(n1,n2,n3,n4,nd1,nd2,nd3,nd4,z,isign,inzee,zw,ncache)

c  accuracy test agains NAG library 
c  (comment the next 8 lines if you do not have NAG):
	if (ntime.eq.1) then
        if (isign.eq.1) call C06GCF(xi,n1*n2*n3*n4,ifail)
        call C06FJF(4,nd,n1*n2*n3*n4,xr,xi,work,lwork,ifail)
        if (isign.eq.1) call C06GCF(xi,n1*n2*n3*n4,ifail)
	call vglp(n1,n2,n3,n4,nd1,nd2,nd3,nd4,
     1       z(1,1,inzee),xr,xi,tta,ttm)
	if (ntime.eq.1) print*,'NAG<>FFT:ttm=,tta=',ttm,tta
	endif


c unix timer (IBM,NEC,SGI)
        t1=mclock()
c Dec timer
C        t1=etime(tarray)
c CRAY timer
C	call second(t1)

	isign=-isign
	do i=1,ntime
	call fft(n1,n2,n3,n4,nd1,nd2,nd3,nd4,z,isign,inzee,zw,ncache)
	enddo

c unix timer (IBM,NEC,SGI)
        t2=mclock()
	time=1.d+1*(t2-t1)/ntime
c Dec timer
C        t2=etime(tarray)
C	time=1.d+3*(t2-t1)/ntime
c CRAY timer
C	call second(t2)
C	time=1.d+3*(t2-t1)/ntime

	call vgl(n1,n2,n3,n4,nd1,nd2,nd3,nd4,z(1,1,inzee),
     1	         n1,n2,n3,n4,zin,1.d0/(n1*n2*n3*n4),tta,ttm)
	if (ntime.eq.1) print*,'Backw<>Forw:ttm=,tta=',ttm,tta
        if (ntime.eq.1 .and. ttm.gt.1.d-8) print*, 'WARNING'

	if (ntime.ne.1) print*, 'Time fft(msec):' ,time,ntime

	if (isign.eq.-1) then
	goto 1234
	endif
	write(6,*) z(1,1,1)

	end





	subroutine init(n1,n2,n3,n4,nd1,nd2,nd3,nd4,zin,z,xr,xi)
	implicit real*8 (a-h,o-z)
	dimension zin(2,n1,n2,n3,n4),z(2,nd1,nd2,nd3,nd4),
     1	          xr(n1,n2,n3,n4),xi(n1,n2,n3,n4)
	do 9763,i4=1,n4
	do 9763,i3=1,n3
	do 9763,i2=1,n2
	do 9763,i1=1,n1
	zin(1,i1,i2,i3,i4) = cos(1.23*float(i1*111+ i2*11 + i3 +i4*4))
	zin(2,i1,i2,i3,i4) = sin(3.21*float(i4*4 + i3*111 + i2*11 + i1))
	z(1,i1,i2,i3,i4) = zin(1,i1,i2,i3,i4) 
	z(2,i1,i2,i3,i4) = zin(2,i1,i2,i3,i4) 
	xr(i1,i2,i3,i4) = zin(1,i1,i2,i3,i4) 
	xi(i1,i2,i3,i4) = zin(2,i1,i2,i3,i4) 
9763	continue
	return
	end

	subroutine vgl(n1,n2,n3,n4,nd1,nd2,nd3,nd4,x,md1,md2,md3,md4,y,
     1                 scale,tta,ttm)
	implicit real*8 (a-h,o-z)
	dimension x(2,nd1,nd2,nd3,nd4),y(2,md1,md2,md3,md4)
	ttm=0.d0
	tta=0.d0
	do 976,i4=1,n4
	do 976,i3=1,n3
	do 976,i2=1,n2
	do 976,i1=1,n1
	ttr=abs(x(1,i1,i2,i3,i4)*
     1       scale-y(1,i1,i2,i3,i4))/abs(y(1,i1,i2,i3,i4))
	tti=abs(x(2,i1,i2,i3,i4)*
     1       scale-y(2,i1,i2,i3,i4))/abs(y(2,i1,i2,i3,i4))
	ttm=max(ttr,tti,ttm)
	tta=tta+ttr+tti
976	continue
	tta=tta/(n1*n2*n3*n4)
	return
	end



	subroutine vglp(n1,n2,n3,n4,nd1,nd2,nd3,nd4,z,xr,xi,tta,ttm)
	implicit real*8 (a-h,o-z)
	dimension z(2,nd1,nd2,nd3,nd4),xr(n1,n2,n3,n4),xi(n1,n2,n3,n4)
	scale=sqrt(1.d0/(n1*n2*n3*n4))
	ttm=0.d0
	tta=0.d0
	do 976,i4=1,n4
	do 976,i3=1,n3
	do 976,i2=1,n2
	do 976,i1=1,n1
	ttr=abs(scale*z(1,i1,i2,i3,i4)-xr(i1,i2,i3,i4))/
     1      abs(xr(i1,i2,i3,i4))
	tti=abs(scale*z(2,i1,i2,i3,i4)-xi(i1,i2,i3,i4))/
     1      abs(xi(i1,i2,i3,i4))
C	write(6,*) scale*z(1,i1,i2,i3,i4),xr(i1,i2,i3,i4)
C	write(6,*) scale*z(2,i1,i2,i3,i4),xi(i1,i2,i3,i4)
	ttm=max(ttr,tti,ttm)
	tta=tta+ttr+tti
976	continue
	tta=tta/(n1*n2*n3*n4)
	return
	end


C FFT PART -----------------------------------------------------------------

	subroutine FFT(n1,n2,n3,n4,nd1,nd2,nd3,nd4,z,isign,inzee,
     1             zw,ncache)
C	CALCULATES THE DISCRETE FOURIERTRANSFORM F(I1,I2,I3)=
C	S_(j1,j2,j3,j4) EXP(isign*i*2*pi*(j1*i1/n1+j2*i2/n2+j3*i3/n3+j4*i4/n4)) 
C                       R(j1,j2,j3,j4)
C	INPUT:
C	    n1,n2,n3,n4:physical dimension of the transform. It must be a 
C		     product of the prime factors 2,3,5, but greater than 3. 
C                    If two ni's are equal it is recommended to place them 
C                    behind each other.
C	    nd1,nd2,nd3,nd4:memory dimension of Z. ndi must always be greater or 
C		        equal than ni. It is recomended to chose ndi=ni 
C		        if ni is odd and ndi=ni+1 if ni is even to obtain 
C                       optimal execution speed. For small ni it can however 
C                       sometimes be advantageous to set ndi=ni
C           inzee=1: first part of Z is data array, second part work array
C           inzee=2: first part of Z is work array, second part data array
C		Z(1,i1,i2,i3,i4,inzee)=real(R(i1,i2,i3,i4))
C		Z(2,i1,i2,i3,i4,inzee)=imag(R(i1,i2,i3,i4))
C	OUTPUT:
C           inzee=1: first part of Z is data array, second part work array
C           inzee=2: first part of Z is work array, second part data array
C		real(F(i1,i2,i3,i4))=Z(1,i1,i2,i3,i4,inzee)
C		imag(F(i1,i2,i3,i4))=Z(2,i1,i2,i3,i4,inzee)
C            inzee on output is in general different from inzee on input
C	THE ARRAYELEMENTS Z( , , , ,3-inzee) ARE USED AS WORKING SPACE
C       On a RISC machine with cache, it is very important to find the optimal 
C       value of NCACHE. NCACHE determines the size of the work array zw, that
C       has to fit into cache. It has therefore to be chosen to equal roughly 
C	half the size of the physical cache in units of real*8 numbers.
C       The optimal value of ncache can easily be determined by numerical 
C       experimentation. A too large value of ncache leads to a dramatic 
C       and sudden decrease of performance, a too small value to a to a 
C       slow and less dramatic decrease of performance. If NCACHE is set 
C       to a value so small, that not even a single one dimensional transform 
C       can be done in the workarray zw, the program stops with an error message.
C	On a vector machine ncache has to be put to 0 
C 	Copyright by Stefan Goedecker, Cornell, Ithaca, USA, March 25, 1994
C 	modified by Stefan Goedecker, Stuttgart, Germany, October 15, 1995
C       Commercial use is prohibited without the explicit permission of the author.


	implicit real*8 (a-h,o-z)
	integer after,before,now
	dimension trig(2,1024),zw(2,ncache/4,2),
     1	z(2,nd1*nd2*nd3*nd4,2),after(20),now(20),before(20)
	if (max(n1,n2,n3,n4).gt.1024) stop '1024'

c vector computer with memory banks:
      if (ncache.eq.0) then

	call ctrig(n4,trig,after,before,now,isign,ic)
	nfft=nd1*nd2*n3
	mm=nd1*nd2*nd3
	do 54093,i=1,ic-1
	call fftstp(mm,nfft,nd4,mm,nd4,z(1,1,inzee),z(1,1,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
54093	inzee=3-inzee
	i=ic
	call fftrot(mm,nfft,nd4,mm,nd4,z(1,1,inzee),z(1,1,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
	inzee=3-inzee

	if (n3.ne.n4) call ctrig(n3,trig,after,before,now,isign,ic)
	nfft=nd4*nd1*n2
	mm=nd4*nd1*nd2
	do 51093,i=1,ic-1
	call fftstp(mm,nfft,nd3,mm,nd3,z(1,1,inzee),z(1,1,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
51093	inzee=3-inzee
	i=ic
	call fftrot(mm,nfft,nd3,mm,nd3,z(1,1,inzee),z(1,1,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
	inzee=3-inzee

	if (n2.ne.n3) call ctrig(n2,trig,after,before,now,isign,ic)
	nfft=nd3*nd4*n1
	mm=nd3*nd4*nd1
	do 52093,i=1,ic-1
	call fftstp(mm,nfft,nd2,mm,nd2,z(1,1,inzee),z(1,1,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
52093	inzee=3-inzee
	i=ic
	call fftrot(mm,nfft,nd2,mm,nd2,z(1,1,inzee),z(1,1,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
	inzee=3-inzee

	if (n1.ne.n2) call ctrig(n1,trig,after,before,now,isign,ic)
	nfft=nd2*nd3*n4
	mm=nd2*nd3*nd4
	do 53093,i=1,ic-1
	call fftstp(mm,nfft,nd1,mm,nd1,z(1,1,inzee),z(1,1,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
53093	inzee=3-inzee
	i=ic
	call fftrot(mm,nfft,nd1,mm,nd1,z(1,1,inzee),z(1,1,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
	inzee=3-inzee

c RISC machine with cache:
      else

C TRANSFORM ALONG T AXIS
	mm=nd1*nd2*nd3
	m=nd4
        lot=max(1,ncache/(4*n4))
C	print*,'lot',lot
	nn=lot
	n=n4
	if (2*n*lot*2.gt.ncache) stop 'ncache1'
	call ctrig(n4,trig,after,before,now,isign,ic)

	ja=1
	jb=nd1*nd2*n3

      if (ic.eq.1) then
	i=ic
	jj=ja*nd4-nd4+1
	nfft=jb-ja+1
	call fftrot(mm,nfft,m,mm,m,z(1,ja,inzee),z(1,jj,3-inzee),
     1	trig,after(i),now(i),before(i),isign)

      else

        do 5000,j=ja,jb,lot
        ma=j
        mb=min(j+(lot-1),jb)
	nfft=mb-ma+1
	jj=j*nd4-nd4+1

	i=1
	inzeep=2
	call fftstp(mm,nfft,m,nn,n,z(1,j,inzee),zw(1,1,3-inzeep),
     1	trig,after(i),now(i),before(i),isign)
	inzeep=1

	do 5093,i=2,ic-1
	call fftstp(nn,nfft,n,nn,n,zw(1,1,inzeep),zw(1,1,3-inzeep),
     1	trig,after(i),now(i),before(i),isign)
5093	inzeep=3-inzeep
	i=ic
	call fftrot(nn,nfft,n,mm,m,zw(1,1,inzeep),z(1,jj,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
5000	continue
      endif

	inzee=3-inzee


C TRANSFORM ALONG Z AXIS
	mm=nd4*nd1*nd2
	m=nd3
        lot=max(1,ncache/(4*n3))
C	print*,'lot',lot
	nn=lot
	n=n3
	if (2*n*lot*2.gt.ncache) stop 'ncache1'
	if (n3.ne.n4) call ctrig(n3,trig,after,before,now,isign,ic)

	ja=1
	jb=nd4*nd1*n2

      if (ic.eq.1) then
	i=ic
	jj=ja*nd3-nd3+1
	nfft=jb-ja+1
	call fftrot(mm,nfft,m,mm,m,z(1,ja,inzee),z(1,jj,3-inzee),
     1	trig,after(i),now(i),before(i),isign)

      else

        do 1000,j=ja,jb,lot
        ma=j
        mb=min(j+(lot-1),jb)
	nfft=mb-ma+1
	jj=j*nd3-nd3+1

	i=1
	inzeep=2
	call fftstp(mm,nfft,m,nn,n,z(1,j,inzee),zw(1,1,3-inzeep),
     1	trig,after(i),now(i),before(i),isign)
	inzeep=1

	do 1093,i=2,ic-1
	call fftstp(nn,nfft,n,nn,n,zw(1,1,inzeep),zw(1,1,3-inzeep),
     1	trig,after(i),now(i),before(i),isign)
1093	inzeep=3-inzeep
	i=ic
	call fftrot(nn,nfft,n,mm,m,zw(1,1,inzeep),z(1,jj,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
1000	continue
      endif

	inzee=3-inzee


C TRANSFORM ALONG Y AXIS
	mm=nd3*nd4*nd1
	m=nd2
        lot=max(1,ncache/(4*n2))
C	print*,'lot',lot
	nn=lot
	n=n2
	if (2*n*lot*2.gt.ncache) stop 'ncache2'
	if (n2.ne.n3) call ctrig(n2,trig,after,before,now,isign,ic)

	ja=1
	jb=nd3*nd4*n1

      if (ic.eq.1) then
	i=ic
	jj=ja*nd2-nd2+1
	nfft=jb-ja+1
	call fftrot(mm,nfft,m,mm,m,z(1,ja,inzee),z(1,jj,3-inzee),
     1	trig,after(i),now(i),before(i),isign)

      else

        do 2000,j=ja,jb,lot
        ma=j
        mb=min(j+(lot-1),jb)
	nfft=mb-ma+1
	jj=j*nd2-nd2+1

	i=1
	inzeep=2
	call fftstp(mm,nfft,m,nn,n,z(1,j,inzee),zw(1,1,3-inzeep),
     1	trig,after(i),now(i),before(i),isign)
	inzeep=1

	do 2093,i=2,ic-1
	call fftstp(nn,nfft,n,nn,n,zw(1,1,inzeep),zw(1,1,3-inzeep),
     1	trig,after(i),now(i),before(i),isign)
2093	inzeep=3-inzeep

	i=ic
	call fftrot(nn,nfft,n,mm,m,zw(1,1,inzeep),z(1,jj,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
2000	continue
      endif
	inzee=3-inzee


C TRANSFORM ALONG X AXIS
	mm=nd2*nd3*nd4
	m=nd1
        lot=max(1,ncache/(4*n1))
C	print*,'lot',lot
	nn=lot
	n=n1
	if (2*n*lot*2.gt.ncache) stop 'ncache3'
	if (n1.ne.n2) call ctrig(n1,trig,after,before,now,isign,ic)

      if (ic.eq.1) then
	i=ic
	j=1
	jp=nd2*nd3*n4
	jj=j*nd1-nd1+1
	nfft=jp-j+1
	call fftrot(mm,nfft,m,mm,m,z(1,j,inzee),z(1,jj,3-inzee),
     1	trig,after(i),now(i),before(i),isign)

      else

        do 3000,j=1,nd2*nd3*n4,lot
        ma=j
        mb=min(j+(lot-1),nd2*nd3*n4)
	nfft=mb-ma+1
	jj=j*nd1-nd1+1

	i=1
	inzeep=2
	call fftstp(mm,nfft,m,nn,n,z(1,j,inzee),zw(1,1,3-inzeep),
     1	trig,after(i),now(i),before(i),isign)
	inzeep=1

	do 3093,i=2,ic-1
	call fftstp(nn,nfft,n,nn,n,zw(1,1,inzeep),zw(1,1,3-inzeep),
     1	trig,after(i),now(i),before(i),isign)
3093	inzeep=3-inzeep
	i=ic
	call fftrot(nn,nfft,n,mm,m,zw(1,1,inzeep),z(1,jj,3-inzee),
     1	trig,after(i),now(i),before(i),isign)
3000	continue
      endif
	inzee=3-inzee


      endif
	return
	end



	subroutine ctrig(n,trig,after,before,now,isign,ic)
C       Copyright by Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
C       modified by Stefan Goedecker, Cornell, Ithaca, USA, March 25, 1994
C       modified by Stefan Goedecker, Stuttgart, Germany, October 6, 1995
C       Commercial use is prohibited without the explicit permission of the author.

	implicit real*8 (a-h,o-z)
	integer after,before
	dimension idata(7,82),now(7),after(7),before(7),
     1  trig(2,1024)

	data idata /
     &      3,   3, 1, 1, 1, 1, 1,       4,   4, 1, 1, 1, 1, 1,
     &      5,   5, 1, 1, 1, 1, 1,       6,   6, 1, 1, 1, 1, 1,
     &      8,   8, 1, 1, 1, 1, 1,       9,   3, 3, 1, 1, 1, 1,
     &     12,   4, 3, 1, 1, 1, 1,      15,   5, 3, 1, 1, 1, 1,
     &     16,   4, 4, 1, 1, 1, 1,      18,   6, 3, 1, 1, 1, 1,
     &     20,   5, 4, 1, 1, 1, 1,      24,   8, 3, 1, 1, 1, 1,
     &     25,   5, 5, 1, 1, 1, 1,      27,   3, 3, 3, 1, 1, 1,
     &     30,   6, 5, 1, 1, 1, 1,      32,   8, 4, 1, 1, 1, 1,
     &     36,   4, 3, 3, 1, 1, 1,      40,   8, 5, 1, 1, 1, 1,
     &     45,   5, 3, 3, 1, 1, 1,      48,   4, 4, 3, 1, 1, 1,
     &     54,   6, 3, 3, 1, 1, 1,      60,   5, 4, 3, 1, 1, 1,
     &     64,   4, 4, 4, 1, 1, 1,      72,   8, 3, 3, 1, 1, 1,
     &     75,   5, 5, 3, 1, 1, 1,      80,   5, 4, 4, 1, 1, 1,
     &     81,   3, 3, 3, 3, 1, 1,      90,   6, 5, 3, 1, 1, 1,
     &     96,   8, 4, 3, 1, 1, 1,     100,   5, 5, 4, 1, 1, 1,
     &    108,   4, 3, 3, 3, 1, 1,     120,   8, 5, 3, 1, 1, 1,
     &    125,   5, 5, 5, 1, 1, 1,     128,   8, 4, 4, 1, 1, 1,
     &    135,   5, 3, 3, 3, 1, 1,     144,   4, 4, 3, 3, 1, 1,
     &    150,   6, 5, 5, 1, 1, 1,     160,   8, 5, 4, 1, 1, 1,
     &    162,   6, 3, 3, 3, 1, 1,     180,   5, 4, 3, 3, 1, 1,
     &    192,   4, 4, 4, 3, 1, 1,     200,   8, 5, 5, 1, 1, 1,
     &    216,   8, 3, 3, 3, 1, 1,     225,   5, 5, 3, 3, 1, 1,
     &    240,   5, 4, 4, 3, 1, 1,     243,   3, 3, 3, 3, 3, 1,
     &    256,   4, 4, 4, 4, 1, 1,     270,   6, 5, 3, 3, 1, 1,
     &    288,   8, 4, 3, 3, 1, 1,     300,   5, 5, 4, 3, 1, 1,
     &    320,   5, 4, 4, 4, 1, 1,     324,   4, 3, 3, 3, 3, 1,
     &    360,   8, 5, 3, 3, 1, 1,     375,   5, 5, 5, 3, 1, 1,
     &    384,   8, 4, 4, 3, 1, 1,     400,   5, 5, 4, 4, 1, 1,
     &    405,   5, 3, 3, 3, 3, 1,     432,   4, 4, 3, 3, 3, 1,
     &    450,   6, 5, 5, 3, 1, 1,     480,   8, 5, 4, 3, 1, 1,
     &    486,   6, 3, 3, 3, 3, 1,     500,   5, 5, 5, 4, 1, 1,
     &    512,   8, 4, 4, 4, 1, 1,     540,   5, 4, 3, 3, 3, 1,
     &    576,   4, 4, 4, 3, 3, 1,     600,   8, 5, 5, 3, 1, 1,
     &    625,   5, 5, 5, 5, 1, 1,     640,   8, 5, 4, 4, 1, 1,
     &    648,   8, 3, 3, 3, 3, 1,     675,   5, 5, 3, 3, 3, 1,
     &    720,   5, 4, 4, 3, 3, 1,     729,   3, 3, 3, 3, 3, 3,
     &    750,   6, 5, 5, 5, 1, 1,     768,   4, 4, 4, 4, 3, 1,
     &    800,   8, 5, 5, 4, 1, 1,     810,   6, 5, 3, 3, 3, 1,
     &    864,   8, 4, 3, 3, 3, 1,     900,   5, 5, 4, 3, 3, 1,
     &    960,   5, 4, 4, 4, 3, 1,     972,   4, 3, 3, 3, 3, 3,
     &   1000,   8, 5, 5, 5, 1, 1,    1024,   4, 4, 4, 4, 4, 1  /
	do 111,i=1,82
	if (n.eq.idata(1,i)) then
	ic=0
	do 11,j=1,6
	itt=idata(1+j,i)
	if (itt.gt.1) then
	ic=ic+1
	now(j)=idata(1+j,i)
	else
	goto 1000
	endif
11	continue
	goto 1000
	endif
111	continue
	print*,'VALUE OF',n,'NOT ALLOWED FOR FFT, ALLOWED VALUES ARE:'
37	format(15(i5))
	write(6,37) (idata(1,i),i=1,82)
	stop
1000	continue

	after(1)=1
	before(ic)=1
	do 22,i=2,ic
	after(i)=after(i-1)*now(i-1)
22	before(ic-i+1)=before(ic-i+2)*now(ic-i+2)

12	format(6(i3))
C	write(6,12) (after(i),i=1,ic)
C	write(6,12) (now(i),i=1,ic)
C	write(6,12) (before(i),i=1,ic)

	twopi=8.d0*datan(1.d0)
	angle=isign*twopi/n
	trig(1,1)=1.d0
	trig(2,1)=0.d0
	do 40,i=1,n-1
	trig(1,i+1)=cos(i*angle)
	trig(2,i+1)=sin(i*angle)
40	continue
check
C	if (mod(n,2).eq.0) then
C	trig(2,n/2+1)=-100000.d0
C	endif

	return
	end


cccccccccccccccccccccccccccccccccccccccccccccccc

 	subroutine fftstp(mm,nfft,m,nn,n,zin,zout,
     1	                 trig,after,now,before,isign)
C 	Copyright by Stefan Goedecker, Cornell, Ithaca, USA, March 25, 1994
C 	modified by Stefan Goedecker, Stuttgart, Germany, October 15, 1995
C       Commercial use is prohibited without the explicit permission of the author.

	implicit real*8 (a-h,o-z)
	integer after,before,atn,atb
	dimension trig(2,1024),zin(2,mm,m),zout(2,nn,n)
	atn=after*now
	atb=after*before

c	 sqrt(.5d0)
	rt2i=0.7071067811865475d0
	if (now.eq.4) then
	if (isign.eq.1) then 
		ia=1
		nin1=ia-after
		nout1=ia-atn
		do 4001,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nin4=nin3+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		nout4=nout3+after
		do 4001,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r2=zin(1,j,nin2)
		s2=zin(2,j,nin2)
		r3=zin(1,j,nin3)
		s3=zin(2,j,nin3)
		r4=zin(1,j,nin4)
		s4=zin(2,j,nin4)
		r=r1 + r3
		s=r2 + r4
		zout(1,j,nout1) = r + s
		zout(1,j,nout3) = r - s
		r=r1 - r3
		s=s2 - s4
		zout(1,j,nout2) = r - s 
		zout(1,j,nout4) = r + s
		r=s1 + s3
		s=s2 + s4
		zout(2,j,nout1) = r + s 
		zout(2,j,nout3) = r - s
		r=s1 - s3
		s=r2 - r4
		zout(2,j,nout2) = r + s 
		zout(2,j,nout4) = r - s
4001		continue
		do 4000,ia=2,after
		ias=ia-1
		if (2*ias.eq.after) then
			nin1=ia-after
			nout1=ia-atn
			do 4010,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			do 4010,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=(r-s)*rt2i
			s2=(r+s)*rt2i
			r3=-zin(2,j,nin3)
			s3=zin(1,j,nin3)
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=-(r + s)*rt2i
			s4=(r - s)*rt2i
			r=r1 + r3
			s=r2 + r4
			zout(1,j,nout1) = r + s
			zout(1,j,nout3) = r - s
			r=r1 - r3
			s=s2 - s4
			zout(1,j,nout2) = r - s 
			zout(1,j,nout4) = r + s
			r=s1 + s3
			s=s2 + s4
			zout(2,j,nout1) = r + s 
			zout(2,j,nout3) = r - s
			r=s1 - s3
			s=r2 - r4
			zout(2,j,nout2) = r + s 
			zout(2,j,nout4) = r - s
4010			continue
		else
			itt=ias*before
			itrig=itt+1
			cr2=trig(1,itrig)
			ci2=trig(2,itrig)
			itrig=itrig+itt
			cr3=trig(1,itrig)
			ci3=trig(2,itrig)
			itrig=itrig+itt
			cr4=trig(1,itrig)
			ci4=trig(2,itrig)
			nin1=ia-after
			nout1=ia-atn
			do 4020,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			do 4020,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=r*cr2 - s*ci2
			s2=r*ci2 + s*cr2
			r=zin(1,j,nin3)
			s=zin(2,j,nin3)
			r3=r*cr3 - s*ci3
			s3=r*ci3 + s*cr3
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=r*cr4 - s*ci4
			s4=r*ci4 + s*cr4
			r=r1 + r3
			s=r2 + r4
			zout(1,j,nout1) = r + s
			zout(1,j,nout3) = r - s
			r=r1 - r3
			s=s2 - s4
			zout(1,j,nout2) = r - s 
			zout(1,j,nout4) = r + s
			r=s1 + s3
			s=s2 + s4
			zout(2,j,nout1) = r + s 
			zout(2,j,nout3) = r - s
			r=s1 - s3
			s=r2 - r4
			zout(2,j,nout2) = r + s 
			zout(2,j,nout4) = r - s
4020			continue
		endif
4000		continue
	else
		ia=1
		nin1=ia-after
		nout1=ia-atn
		do 4101,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nin4=nin3+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		nout4=nout3+after
		do 4101,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r2=zin(1,j,nin2)
		s2=zin(2,j,nin2)
		r3=zin(1,j,nin3)
		s3=zin(2,j,nin3)
		r4=zin(1,j,nin4)
		s4=zin(2,j,nin4)
		r=r1 + r3
		s=r2 + r4
		zout(1,j,nout1) = r + s
		zout(1,j,nout3) = r - s
		r=r1 - r3
		s=s2 - s4
		zout(1,j,nout2) = r + s
		zout(1,j,nout4) = r - s
		r=s1 + s3
		s=s2 + s4
		zout(2,j,nout1) = r + s
		zout(2,j,nout3) = r - s
		r=s1 - s3
		s=r2 - r4
		zout(2,j,nout2) = r - s
		zout(2,j,nout4) = r + s
4101		continue
		do 4100,ia=2,after
		ias=ia-1
		if (2*ias.eq.after) then
			nin1=ia-after
			nout1=ia-atn
			do 4110,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			do 4110,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=(r + s)*rt2i
			s2=(s - r)*rt2i
			r3=zin(2,j,nin3)
			s3=-zin(1,j,nin3)
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=(s - r)*rt2i
			s4=-(r + s)*rt2i
			r=r1 + r3
			s=r2 + r4
			zout(1,j,nout1) = r + s
			zout(1,j,nout3) = r - s
			r=r1 - r3
			s=s2 - s4
			zout(1,j,nout2) = r + s
			zout(1,j,nout4) = r - s
			r=s1 + s3
			s=s2 + s4
			zout(2,j,nout1) = r + s
			zout(2,j,nout3) = r - s
			r=s1 - s3
			s=r2 - r4
			zout(2,j,nout2) = r - s
			zout(2,j,nout4) = r + s
4110			continue
		else
			itt=ias*before
			itrig=itt+1
			cr2=trig(1,itrig)
			ci2=trig(2,itrig)
			itrig=itrig+itt
			cr3=trig(1,itrig)
			ci3=trig(2,itrig)
			itrig=itrig+itt
			cr4=trig(1,itrig)
			ci4=trig(2,itrig)
			nin1=ia-after
			nout1=ia-atn
			do 4120,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			do 4120,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=r*cr2 - s*ci2
			s2=r*ci2 + s*cr2
			r=zin(1,j,nin3)
			s=zin(2,j,nin3)
			r3=r*cr3 - s*ci3
			s3=r*ci3 + s*cr3
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=r*cr4 - s*ci4
			s4=r*ci4 + s*cr4
			r=r1 + r3
			s=r2 + r4
			zout(1,j,nout1) = r + s
			zout(1,j,nout3) = r - s
			r=r1 - r3
			s=s2 - s4
			zout(1,j,nout2) = r + s
			zout(1,j,nout4) = r - s
			r=s1 + s3
			s=s2 + s4
			zout(2,j,nout1) = r + s
			zout(2,j,nout3) = r - s
			r=s1 - s3
			s=r2 - r4
			zout(2,j,nout2) = r - s
			zout(2,j,nout4) = r + s
4120			continue
		endif
4100		continue
	endif
	else if (now.eq.8) then
	if (isign.eq.-1) then 
		ia=1
			nin1=ia-after
			nout1=ia-atn
			do 8120,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nin5=nin4+atb
			nin6=nin5+atb
			nin7=nin6+atb
			nin8=nin7+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			nout5=nout4+after
			nout6=nout5+after
			nout7=nout6+after
			nout8=nout7+after
			do 8120,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r2=zin(1,j,nin2)
			s2=zin(2,j,nin2)
			r3=zin(1,j,nin3)
			s3=zin(2,j,nin3)
			r4=zin(1,j,nin4)
			s4=zin(2,j,nin4)
			r5=zin(1,j,nin5)
			s5=zin(2,j,nin5)
			r6=zin(1,j,nin6)
			s6=zin(2,j,nin6)
			r7=zin(1,j,nin7)
			s7=zin(2,j,nin7)
			r8=zin(1,j,nin8)
			s8=zin(2,j,nin8)
			r=r1 + r5
			s=r3 + r7
			ap=r + s
			am=r - s
			r=r2 + r6
			s=r4 + r8
			bp=r + s
			bm=r - s
			r=s1 + s5
			s=s3 + s7
			cp=r + s
			cm=r - s
			r=s2 + s6
			s=s4 + s8
			dp=r + s
			dm=r - s
			zout(1,j,nout1) = ap + bp
			zout(2,j,nout1) = cp + dp
			zout(1,j,nout5) = ap - bp
			zout(2,j,nout5) = cp - dp
			zout(1,j,nout3) = am + dm
			zout(2,j,nout3) = cm - bm
			zout(1,j,nout7) = am - dm
			zout(2,j,nout7) = cm + bm
			r=r1 - r5
			s=s3 - s7
			ap=r + s
			am=r - s
			r=s1 - s5
			s=r3 - r7
			bp=r + s
			bm=r - s
			r=s4 - s8
			s=r2 - r6
			cp=r + s
			cm=r - s
			r=s2 - s6
			s=r4 - r8
			dp=r + s
			dm=r - s
			r = ( cp + dm)*rt2i
			s = (-cp + dm)*rt2i
			cp= ( cm + dp)*rt2i
			dp = ( cm - dp)*rt2i
			zout(1,j,nout2) = ap + r
			zout(2,j,nout2) = bm + s
			zout(1,j,nout6) = ap - r
			zout(2,j,nout6) = bm - s
			zout(1,j,nout4) = am + cp
			zout(2,j,nout4) = bp + dp
			zout(1,j,nout8) = am - cp
			zout(2,j,nout8) = bp - dp
8120			continue
	else
		ia=1
			nin1=ia-after
			nout1=ia-atn
			do 8121,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nin5=nin4+atb
			nin6=nin5+atb
			nin7=nin6+atb
			nin8=nin7+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			nout5=nout4+after
			nout6=nout5+after
			nout7=nout6+after
			nout8=nout7+after
			do 8121,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r2=zin(1,j,nin2)
			s2=zin(2,j,nin2)
			r3=zin(1,j,nin3)
			s3=zin(2,j,nin3)
			r4=zin(1,j,nin4)
			s4=zin(2,j,nin4)
			r5=zin(1,j,nin5)
			s5=zin(2,j,nin5)
			r6=zin(1,j,nin6)
			s6=zin(2,j,nin6)
			r7=zin(1,j,nin7)
			s7=zin(2,j,nin7)
			r8=zin(1,j,nin8)
			s8=zin(2,j,nin8)
			r=r1 + r5
			s=r3 + r7
			ap=r + s
			am=r - s
			r=r2 + r6
			s=r4 + r8
			bp=r + s
			bm=r - s
			r=s1 + s5
			s=s3 + s7
			cp=r + s
			cm=r - s
			r=s2 + s6
			s=s4 + s8
			dp=r + s
			dm=r - s
			zout(1,j,nout1) = ap + bp
			zout(2,j,nout1) = cp + dp
			zout(1,j,nout5) = ap - bp
			zout(2,j,nout5) = cp - dp
			zout(1,j,nout3) = am - dm
			zout(2,j,nout3) = cm + bm
			zout(1,j,nout7) = am + dm
			zout(2,j,nout7) = cm - bm
			r= r1 - r5
			s=-s3 + s7
			ap=r + s
			am=r - s
			r=s1 - s5
			s=r7 - r3
			bp=r + s
			bm=r - s
			r=-s4 + s8
			s= r2 - r6
			cp=r + s
			cm=r - s
			r=-s2 + s6
			s= r4 - r8
			dp=r + s
			dm=r - s
			r = ( cp + dm)*rt2i
			s = ( cp - dm)*rt2i
			cp= ( cm + dp)*rt2i
			dp= (-cm + dp)*rt2i
			zout(1,j,nout2) = ap + r
			zout(2,j,nout2) = bm + s
			zout(1,j,nout6) = ap - r
			zout(2,j,nout6) = bm - s
			zout(1,j,nout4) = am + cp
			zout(2,j,nout4) = bp + dp
			zout(1,j,nout8) = am - cp
			zout(2,j,nout8) = bp - dp
8121			continue
	endif
	else if (now.eq.3) then 
c	 .5d0*sqrt(3.d0)
	bb=isign*0.8660254037844387d0
	ia=1
	nin1=ia-after
	nout1=ia-atn
	do 3001,ib=1,before
	nin1=nin1+after
	nin2=nin1+atb
	nin3=nin2+atb
	nout1=nout1+atn
	nout2=nout1+after
	nout3=nout2+after
	do 3001,j=1,nfft
	r1=zin(1,j,nin1)
	s1=zin(2,j,nin1)
	r2=zin(1,j,nin2)
	s2=zin(2,j,nin2)
	r3=zin(1,j,nin3)
	s3=zin(2,j,nin3)
	r=r2 + r3
	s=s2 + s3
	zout(1,j,nout1) = r + r1
	zout(2,j,nout1) = s + s1
	r1=r1 - .5d0*r
	s1=s1 - .5d0*s
	r2=bb*(r2-r3)
	s2=bb*(s2-s3)
	zout(1,j,nout2) = r1 - s2 
	zout(2,j,nout2) = s1 + r2
	zout(1,j,nout3) = r1 + s2 
	zout(2,j,nout3) = s1 - r2
3001	continue
	do 3000,ia=2,after
	ias=ia-1
	if (4*ias.eq.3*after) then
	if (isign.eq.1) then
		nin1=ia-after
		nout1=ia-atn
		do 3010,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		do 3010,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r2=-zin(2,j,nin2)
		s2=zin(1,j,nin2)
		r3=-zin(1,j,nin3)
		s3=-zin(2,j,nin3)
		r=r2 + r3
		s=s2 + s3
		zout(1,j,nout1) = r + r1
		zout(2,j,nout1) = s + s1
		r1=r1 - .5d0*r
		s1=s1 - .5d0*s	
		r2=bb*(r2-r3)	
		s2=bb*(s2-s3)
		zout(1,j,nout2) = r1 - s2 
		zout(2,j,nout2) = s1 + r2
		zout(1,j,nout3) = r1 + s2 
		zout(2,j,nout3) = s1 - r2
3010		continue
	else
		nin1=ia-after
		nout1=ia-atn
		do 3020,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		do 3020,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r2=zin(2,j,nin2)
		s2=-zin(1,j,nin2)
		r3=-zin(1,j,nin3)
		s3=-zin(2,j,nin3)
		r=r2 + r3
		s=s2 + s3
		zout(1,j,nout1) = r + r1
		zout(2,j,nout1) = s + s1
		r1=r1 - .5d0*r
		s1=s1 - .5d0*s	
		r2=bb*(r2-r3)	
		s2=bb*(s2-s3)
		zout(1,j,nout2) = r1 - s2 
		zout(2,j,nout2) = s1 + r2
		zout(1,j,nout3) = r1 + s2 
		zout(2,j,nout3) = s1 - r2
3020		continue
	endif
	else if (8*ias.eq.3*after) then
	if (isign.eq.1) then
		nin1=ia-after
		nout1=ia-atn
		do 3030,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		do 3030,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r=zin(1,j,nin2)
		s=zin(2,j,nin2)
		r2=(r - s)*rt2i
		s2=(r + s)*rt2i
		r3=-zin(2,j,nin3)
		s3=zin(1,j,nin3) 
		r=r2 + r3
		s=s2 + s3
		zout(1,j,nout1) = r + r1
		zout(2,j,nout1) = s + s1
		r1=r1 - .5d0*r
		s1=s1 - .5d0*s	
		r2=bb*(r2-r3)	
		s2=bb*(s2-s3)
		zout(1,j,nout2) = r1 - s2 
		zout(2,j,nout2) = s1 + r2
		zout(1,j,nout3) = r1 + s2 
		zout(2,j,nout3) = s1 - r2
3030		continue
	else
		nin1=ia-after
		nout1=ia-atn
		do 3040,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		do 3040,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r=zin(1,j,nin2)
		s=zin(2,j,nin2)
		r2=(r + s)*rt2i
		s2=(-r + s)*rt2i
		r3=zin(2,j,nin3)
		s3=-zin(1,j,nin3)
		r=r2 + r3
		s=s2 + s3
		zout(1,j,nout1) = r + r1
		zout(2,j,nout1) = s + s1
		r1=r1 - .5d0*r
		s1=s1 - .5d0*s	
		r2=bb*(r2-r3)	
		s2=bb*(s2-s3)
		zout(1,j,nout2) = r1 - s2 
		zout(2,j,nout2) = s1 + r2
		zout(1,j,nout3) = r1 + s2 
		zout(2,j,nout3) = s1 - r2
3040		continue
	endif
	else
	itt=ias*before
	itrig=itt+1
	cr2=trig(1,itrig)
	ci2=trig(2,itrig)
	itrig=itrig+itt
	cr3=trig(1,itrig)
	ci3=trig(2,itrig)
	nin1=ia-after
	nout1=ia-atn
	do 3090,ib=1,before
	nin1=nin1+after
	nin2=nin1+atb
	nin3=nin2+atb
	nout1=nout1+atn
	nout2=nout1+after
	nout3=nout2+after
	do 3090,j=1,nfft
	r1=zin(1,j,nin1)
	s1=zin(2,j,nin1)
	r=zin(1,j,nin2)
	s=zin(2,j,nin2)
	r2=r*cr2 - s*ci2
	s2=r*ci2 + s*cr2
	r=zin(1,j,nin3)
	s=zin(2,j,nin3)
	r3=r*cr3 - s*ci3
	s3=r*ci3 + s*cr3
	r=r2 + r3
	s=s2 + s3
	zout(1,j,nout1) = r + r1
	zout(2,j,nout1) = s + s1
	r1=r1 - .5d0*r
	s1=s1 - .5d0*s
	r2=bb*(r2-r3)
	s2=bb*(s2-s3)
	zout(1,j,nout2) = r1 - s2 
	zout(2,j,nout2) = s1 + r2
	zout(1,j,nout3) = r1 + s2 
	zout(2,j,nout3) = s1 - r2
3090	continue
	endif
3000	continue
	else if (now.eq.5) then
c	 cos(2.d0*pi/5.d0)
	cos2=0.3090169943749474d0
c	 cos(4.d0*pi/5.d0)
	cos4=-0.8090169943749474d0
c	sin(2.d0*pi/5.d0)
	sin2=isign*0.9510565162951536d0
c	 sin(4.d0*pi/5.d0)
	sin4=isign*0.5877852522924731d0
	ia=1
	nin1=ia-after
	nout1=ia-atn
	do 5001,ib=1,before
	nin1=nin1+after
	nin2=nin1+atb
	nin3=nin2+atb
	nin4=nin3+atb
	nin5=nin4+atb
	nout1=nout1+atn
	nout2=nout1+after
	nout3=nout2+after
	nout4=nout3+after
	nout5=nout4+after
	do 5001,j=1,nfft
	r1=zin(1,j,nin1)
	s1=zin(2,j,nin1)
	r2=zin(1,j,nin2)
	s2=zin(2,j,nin2)
	r3=zin(1,j,nin3)
	s3=zin(2,j,nin3)
	r4=zin(1,j,nin4)
	s4=zin(2,j,nin4)
	r5=zin(1,j,nin5)
	s5=zin(2,j,nin5)
	r25 = r2 + r5
	r34 = r3 + r4
	s25 = s2 - s5
	s34 = s3 - s4
	zout(1,j,nout1) = r1 + r25 + r34
	r = cos2*r25 + cos4*r34 + r1
	s = sin2*s25 + sin4*s34
	zout(1,j,nout2) = r - s
	zout(1,j,nout5) = r + s
	r = cos4*r25 + cos2*r34 + r1
	s = sin4*s25 - sin2*s34
	zout(1,j,nout3) = r - s
	zout(1,j,nout4) = r + s
	r25 = r2 - r5
	r34 = r3 - r4
	s25 = s2 + s5
	s34 = s3 + s4
	zout(2,j,nout1) = s1 + s25 + s34
	r = cos2*s25 + cos4*s34 + s1
	s = sin2*r25 + sin4*r34
	zout(2,j,nout2) = r + s
	zout(2,j,nout5) = r - s
	r = cos4*s25 + cos2*s34 + s1
	s = sin4*r25 - sin2*r34
	zout(2,j,nout3) = r + s
	zout(2,j,nout4) = r - s
5001	continue
	do 5000,ia=2,after
	ias=ia-1
	if (8*ias.eq.5*after) then
		if (isign.eq.1) then
			nin1=ia-after
			nout1=ia-atn
			do 5010,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nin5=nin4+atb	
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			nout5=nout4+after
			do 5010,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=(r - s)*rt2i
			s2=(r + s)*rt2i
			r3=-zin(2,j,nin3)
			s3=zin(1,j,nin3) 
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=-(r + s)*rt2i
			s4=(r - s)*rt2i
			r5=-zin(1,j,nin5)
			s5=-zin(2,j,nin5)
			r25 = r2 + r5
			r34 = r3 + r4
			s25 = s2 - s5
			s34 = s3 - s4
			zout(1,j,nout1) = r1 + r25 + r34
			r = cos2*r25 + cos4*r34 + r1
			s = sin2*s25 + sin4*s34
			zout(1,j,nout2) = r - s
			zout(1,j,nout5) = r + s
			r = cos4*r25 + cos2*r34 + r1
			s = sin4*s25 - sin2*s34
			zout(1,j,nout3) = r - s
			zout(1,j,nout4) = r + s
			r25 = r2 - r5
			r34 = r3 - r4
			s25 = s2 + s5
			s34 = s3 + s4
			zout(2,j,nout1) = s1 + s25 + s34
			r = cos2*s25 + cos4*s34 + s1
			s = sin2*r25 + sin4*r34
			zout(2,j,nout2) = r + s
			zout(2,j,nout5) = r - s
			r = cos4*s25 + cos2*s34 + s1
			s = sin4*r25 - sin2*r34
			zout(2,j,nout3) = r + s
			zout(2,j,nout4) = r - s
5010			continue
		else
			nin1=ia-after
			nout1=ia-atn
			do 5020,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nin5=nin4+atb	
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			nout5=nout4+after
			do 5020,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=(r + s)*rt2i
			s2=(-r + s)*rt2i
			r3=zin(2,j,nin3)
			s3=-zin(1,j,nin3)
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=(s - r)*rt2i
			s4=-(r + s)*rt2i
			r5=-zin(1,j,nin5)
			s5=-zin(2,j,nin5)
			r25 = r2 + r5
			r34 = r3 + r4
			s25 = s2 - s5
			s34 = s3 - s4
			zout(1,j,nout1) = r1 + r25 + r34
			r = cos2*r25 + cos4*r34 + r1
			s = sin2*s25 + sin4*s34
			zout(1,j,nout2) = r - s
			zout(1,j,nout5) = r + s
			r = cos4*r25 + cos2*r34 + r1
			s = sin4*s25 - sin2*s34
			zout(1,j,nout3) = r - s
			zout(1,j,nout4) = r + s
			r25 = r2 - r5
			r34 = r3 - r4
			s25 = s2 + s5
			s34 = s3 + s4
			zout(2,j,nout1) = s1 + s25 + s34
			r = cos2*s25 + cos4*s34 + s1
			s = sin2*r25 + sin4*r34
			zout(2,j,nout2) = r + s
			zout(2,j,nout5) = r - s
			r = cos4*s25 + cos2*s34 + s1
			s = sin4*r25 - sin2*r34
			zout(2,j,nout3) = r + s
			zout(2,j,nout4) = r - s
5020			continue
		endif
	else
		ias=ia-1
		itt=ias*before
		itrig=itt+1
		cr2=trig(1,itrig)
		ci2=trig(2,itrig)
		itrig=itrig+itt
		cr3=trig(1,itrig)
		ci3=trig(2,itrig)
		itrig=itrig+itt
		cr4=trig(1,itrig)
		ci4=trig(2,itrig)
		itrig=itrig+itt
		cr5=trig(1,itrig)
		ci5=trig(2,itrig)
		nin1=ia-after
		nout1=ia-atn
		do 5100,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nin4=nin3+atb
		nin5=nin4+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		nout4=nout3+after
		nout5=nout4+after
		do 5100,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r=zin(1,j,nin2)
		s=zin(2,j,nin2)
		r2=r*cr2 - s*ci2
		s2=r*ci2 + s*cr2
		r=zin(1,j,nin3)
		s=zin(2,j,nin3)
		r3=r*cr3 - s*ci3
		s3=r*ci3 + s*cr3
		r=zin(1,j,nin4)
		s=zin(2,j,nin4)
		r4=r*cr4 - s*ci4
		s4=r*ci4 + s*cr4
		r=zin(1,j,nin5)
		s=zin(2,j,nin5)
		r5=r*cr5 - s*ci5
		s5=r*ci5 + s*cr5
		r25 = r2 + r5
		r34 = r3 + r4
		s25 = s2 - s5
		s34 = s3 - s4
		zout(1,j,nout1) = r1 + r25 + r34
		r = cos2*r25 + cos4*r34 + r1
		s = sin2*s25 + sin4*s34
		zout(1,j,nout2) = r - s
		zout(1,j,nout5) = r + s
		r = cos4*r25 + cos2*r34 + r1
		s = sin4*s25 - sin2*s34
		zout(1,j,nout3) = r - s
		zout(1,j,nout4) = r + s
		r25 = r2 - r5
		r34 = r3 - r4
		s25 = s2 + s5
		s34 = s3 + s4
		zout(2,j,nout1) = s1 + s25 + s34
		r = cos2*s25 + cos4*s34 + s1
		s = sin2*r25 + sin4*r34
		zout(2,j,nout2) = r + s
		zout(2,j,nout5) = r - s
		r = cos4*s25 + cos2*s34 + s1
		s = sin4*r25 - sin2*r34
		zout(2,j,nout3) = r + s
		zout(2,j,nout4) = r - s
5100		continue
	endif
5000	continue
       else if (now.eq.6) then
c	 .5d0*sqrt(3.d0)
	bb=isign*0.8660254037844387d0

	ia=1
	nin1=ia-after
	nout1=ia-atn
	do 6120,ib=1,before
	nin1=nin1+after
	nin2=nin1+atb
	nin3=nin2+atb
	nin4=nin3+atb
	nin5=nin4+atb
	nin6=nin5+atb
	nout1=nout1+atn
	nout2=nout1+after
	nout3=nout2+after
	nout4=nout3+after
	nout5=nout4+after
	nout6=nout5+after
	do 6120,j=1,nfft
	r2=zin(1,j,nin3)
	s2=zin(2,j,nin3)
	r3=zin(1,j,nin5)
	s3=zin(2,j,nin5)
	r=r2 + r3
	s=s2 + s3
	r1=zin(1,j,nin1)
	s1=zin(2,j,nin1)
	ur1 = r + r1
	ui1 = s + s1
	r1=r1 - .5d0*r
	s1=s1 - .5d0*s
	r=r2-r3
	s=s2-s3
	ur2 = r1 - s*bb
	ui2 = s1 + r*bb
	ur3 = r1 + s*bb
	ui3 = s1 - r*bb

	r2=zin(1,j,nin6)
	s2=zin(2,j,nin6)
	r3=zin(1,j,nin2)
	s3=zin(2,j,nin2)
	r=r2 + r3
	s=s2 + s3
	r1=zin(1,j,nin4)
	s1=zin(2,j,nin4)
	vr1 = r + r1
	vi1 = s + s1
	r1=r1 - .5d0*r
	s1=s1 - .5d0*s
	r=r2-r3
	s=s2-s3
	vr2 = r1 - s*bb
	vi2 = s1 + r*bb
	vr3 = r1 + s*bb
	vi3 = s1 - r*bb

	zout(1,j,nout1)=ur1+vr1
	zout(2,j,nout1)=ui1+vi1
	zout(1,j,nout5)=ur2+vr2
	zout(2,j,nout5)=ui2+vi2
	zout(1,j,nout3)=ur3+vr3
	zout(2,j,nout3)=ui3+vi3
	zout(1,j,nout4)=ur1-vr1
	zout(2,j,nout4)=ui1-vi1
	zout(1,j,nout2)=ur2-vr2
	zout(2,j,nout2)=ui2-vi2
	zout(1,j,nout6)=ur3-vr3
	zout(2,j,nout6)=ui3-vi3
6120	continue

	else 
	stop 'error fftstp'
	endif

	return
	end



 	subroutine fftrot(mm,nfft,m,nn,n,zin,zout,
     1	                 trig,after,now,before,isign)
C 	Copyright by Stefan Goedecker, Cornell, Ithaca, USA, March 25, 1994
C 	modified by Stefan Goedecker, Stuttgart, Germany, October 15, 1995
C       Commercial use is prohibited without the explicit permission of the author.

	implicit real*8 (a-h,o-z)
	integer after,before,atn,atb
	dimension trig(2,1024),zin(2,mm,m),zout(2,n,nn)
	atn=after*now
	atb=after*before

c	 sqrt(.5d0)
	rt2i=0.7071067811865475d0
	if (now.eq.4) then
	if (isign.eq.1) then 
		ia=1
		nin1=ia-after
		nout1=ia-atn
		do 4001,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nin4=nin3+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		nout4=nout3+after
		do 4001,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r2=zin(1,j,nin2)
		s2=zin(2,j,nin2)
		r3=zin(1,j,nin3)
		s3=zin(2,j,nin3)
		r4=zin(1,j,nin4)
		s4=zin(2,j,nin4)
		r=r1 + r3
		s=r2 + r4
		zout(1,nout1,j) = r + s
		zout(1,nout3,j) = r - s
		r=r1 - r3
		s=s2 - s4
		zout(1,nout2,j) = r - s 
		zout(1,nout4,j) = r + s
		r=s1 + s3
		s=s2 + s4
		zout(2,nout1,j) = r + s 
		zout(2,nout3,j) = r - s
		r=s1 - s3
		s=r2 - r4
		zout(2,nout2,j) = r + s 
		zout(2,nout4,j) = r - s
4001		continue
		do 4000,ia=2,after
		ias=ia-1
		if (2*ias.eq.after) then
			nin1=ia-after
			nout1=ia-atn
			do 4010,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			do 4010,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=(r-s)*rt2i
			s2=(r+s)*rt2i
			r3=-zin(2,j,nin3)
			s3=zin(1,j,nin3)
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=-(r + s)*rt2i
			s4=(r - s)*rt2i
			r=r1 + r3
			s=r2 + r4
			zout(1,nout1,j) = r + s
			zout(1,nout3,j) = r - s
			r=r1 - r3
			s=s2 - s4
			zout(1,nout2,j) = r - s 
			zout(1,nout4,j) = r + s
			r=s1 + s3
			s=s2 + s4
			zout(2,nout1,j) = r + s 
			zout(2,nout3,j) = r - s
			r=s1 - s3
			s=r2 - r4
			zout(2,nout2,j) = r + s 
			zout(2,nout4,j) = r - s
4010			continue
		else
			itt=ias*before
			itrig=itt+1
			cr2=trig(1,itrig)
			ci2=trig(2,itrig)
			itrig=itrig+itt
			cr3=trig(1,itrig)
			ci3=trig(2,itrig)
			itrig=itrig+itt
			cr4=trig(1,itrig)
			ci4=trig(2,itrig)
			nin1=ia-after
			nout1=ia-atn
			do 4020,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			do 4020,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=r*cr2 - s*ci2
			s2=r*ci2 + s*cr2
			r=zin(1,j,nin3)
			s=zin(2,j,nin3)
			r3=r*cr3 - s*ci3
			s3=r*ci3 + s*cr3
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=r*cr4 - s*ci4
			s4=r*ci4 + s*cr4
			r=r1 + r3
			s=r2 + r4
			zout(1,nout1,j) = r + s
			zout(1,nout3,j) = r - s
			r=r1 - r3
			s=s2 - s4
			zout(1,nout2,j) = r - s 
			zout(1,nout4,j) = r + s
			r=s1 + s3
			s=s2 + s4
			zout(2,nout1,j) = r + s 
			zout(2,nout3,j) = r - s
			r=s1 - s3
			s=r2 - r4
			zout(2,nout2,j) = r + s 
			zout(2,nout4,j) = r - s
4020			continue
		endif
4000		continue
	else
		ia=1
		nin1=ia-after
		nout1=ia-atn
		do 4101,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nin4=nin3+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		nout4=nout3+after
		do 4101,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r2=zin(1,j,nin2)
		s2=zin(2,j,nin2)
		r3=zin(1,j,nin3)
		s3=zin(2,j,nin3)
		r4=zin(1,j,nin4)
		s4=zin(2,j,nin4)
		r=r1 + r3
		s=r2 + r4
		zout(1,nout1,j) = r + s
		zout(1,nout3,j) = r - s
		r=r1 - r3
		s=s2 - s4
		zout(1,nout2,j) = r + s
		zout(1,nout4,j) = r - s
		r=s1 + s3
		s=s2 + s4
		zout(2,nout1,j) = r + s
		zout(2,nout3,j) = r - s
		r=s1 - s3
		s=r2 - r4
		zout(2,nout2,j) = r - s
		zout(2,nout4,j) = r + s
4101		continue
		do 4100,ia=2,after
		ias=ia-1
		if (2*ias.eq.after) then
			nin1=ia-after
			nout1=ia-atn
			do 4110,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			do 4110,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=(r + s)*rt2i
			s2=(s - r)*rt2i
			r3=zin(2,j,nin3)
			s3=-zin(1,j,nin3)
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=(s - r)*rt2i
			s4=-(r + s)*rt2i
			r=r1 + r3
			s=r2 + r4
			zout(1,nout1,j) = r + s
			zout(1,nout3,j) = r - s
			r=r1 - r3
			s=s2 - s4
			zout(1,nout2,j) = r + s
			zout(1,nout4,j) = r - s
			r=s1 + s3
			s=s2 + s4
			zout(2,nout1,j) = r + s
			zout(2,nout3,j) = r - s
			r=s1 - s3
			s=r2 - r4
			zout(2,nout2,j) = r - s
			zout(2,nout4,j) = r + s
4110			continue
		else
			itt=ias*before
			itrig=itt+1
			cr2=trig(1,itrig)
			ci2=trig(2,itrig)
			itrig=itrig+itt
			cr3=trig(1,itrig)
			ci3=trig(2,itrig)
			itrig=itrig+itt
			cr4=trig(1,itrig)
			ci4=trig(2,itrig)
			nin1=ia-after
			nout1=ia-atn
			do 4120,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			do 4120,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=r*cr2 - s*ci2
			s2=r*ci2 + s*cr2
			r=zin(1,j,nin3)
			s=zin(2,j,nin3)
			r3=r*cr3 - s*ci3
			s3=r*ci3 + s*cr3
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=r*cr4 - s*ci4
			s4=r*ci4 + s*cr4
			r=r1 + r3
			s=r2 + r4
			zout(1,nout1,j) = r + s
			zout(1,nout3,j) = r - s
			r=r1 - r3
			s=s2 - s4
			zout(1,nout2,j) = r + s
			zout(1,nout4,j) = r - s
			r=s1 + s3
			s=s2 + s4
			zout(2,nout1,j) = r + s
			zout(2,nout3,j) = r - s
			r=s1 - s3
			s=r2 - r4
			zout(2,nout2,j) = r - s
			zout(2,nout4,j) = r + s
4120			continue
		endif
4100		continue
	endif

	else if (now.eq.8) then
	if (isign.eq.-1) then 
		ia=1
			nin1=ia-after
			nout1=ia-atn
			do 8120,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nin5=nin4+atb
			nin6=nin5+atb
			nin7=nin6+atb
			nin8=nin7+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			nout5=nout4+after
			nout6=nout5+after
			nout7=nout6+after
			nout8=nout7+after
			do 8120,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r2=zin(1,j,nin2)
			s2=zin(2,j,nin2)
			r3=zin(1,j,nin3)
			s3=zin(2,j,nin3)
			r4=zin(1,j,nin4)
			s4=zin(2,j,nin4)
			r5=zin(1,j,nin5)
			s5=zin(2,j,nin5)
			r6=zin(1,j,nin6)
			s6=zin(2,j,nin6)
			r7=zin(1,j,nin7)
			s7=zin(2,j,nin7)
			r8=zin(1,j,nin8)
			s8=zin(2,j,nin8)
			r=r1 + r5
			s=r3 + r7
			ap=r + s
			am=r - s
			r=r2 + r6
			s=r4 + r8
			bp=r + s
			bm=r - s
			r=s1 + s5
			s=s3 + s7
			cp=r + s
			cm=r - s
			r=s2 + s6
			s=s4 + s8
			dp=r + s
			dm=r - s
			zout(1,nout1,j) = ap + bp
			zout(2,nout1,j) = cp + dp
			zout(1,nout5,j) = ap - bp
			zout(2,nout5,j) = cp - dp
			zout(1,nout3,j) = am + dm
			zout(2,nout3,j) = cm - bm
			zout(1,nout7,j) = am - dm
			zout(2,nout7,j) = cm + bm
			r=r1 - r5
			s=s3 - s7
			ap=r + s
			am=r - s
			r=s1 - s5
			s=r3 - r7
			bp=r + s
			bm=r - s
			r=s4 - s8
			s=r2 - r6
			cp=r + s
			cm=r - s
			r=s2 - s6
			s=r4 - r8
			dp=r + s
			dm=r - s
			r = ( cp + dm)*rt2i
			s = (-cp + dm)*rt2i
			cp= ( cm + dp)*rt2i
			dp = ( cm - dp)*rt2i
			zout(1,nout2,j) = ap + r
			zout(2,nout2,j) = bm + s
			zout(1,nout6,j) = ap - r
			zout(2,nout6,j) = bm - s
			zout(1,nout4,j) = am + cp
			zout(2,nout4,j) = bp + dp
			zout(1,nout8,j) = am - cp
			zout(2,nout8,j) = bp - dp
8120			continue
	else
		ia=1
			nin1=ia-after
			nout1=ia-atn
			do 8121,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nin5=nin4+atb
			nin6=nin5+atb
			nin7=nin6+atb
			nin8=nin7+atb
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			nout5=nout4+after
			nout6=nout5+after
			nout7=nout6+after
			nout8=nout7+after
			do 8121,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r2=zin(1,j,nin2)
			s2=zin(2,j,nin2)
			r3=zin(1,j,nin3)
			s3=zin(2,j,nin3)
			r4=zin(1,j,nin4)
			s4=zin(2,j,nin4)
			r5=zin(1,j,nin5)
			s5=zin(2,j,nin5)
			r6=zin(1,j,nin6)
			s6=zin(2,j,nin6)
			r7=zin(1,j,nin7)
			s7=zin(2,j,nin7)
			r8=zin(1,j,nin8)
			s8=zin(2,j,nin8)
			r=r1 + r5
			s=r3 + r7
			ap=r + s
			am=r - s
			r=r2 + r6
			s=r4 + r8
			bp=r + s
			bm=r - s
			r=s1 + s5
			s=s3 + s7
			cp=r + s
			cm=r - s
			r=s2 + s6
			s=s4 + s8
			dp=r + s
			dm=r - s
			zout(1,nout1,j) = ap + bp
			zout(2,nout1,j) = cp + dp
			zout(1,nout5,j) = ap - bp
			zout(2,nout5,j) = cp - dp
			zout(1,nout3,j) = am - dm
			zout(2,nout3,j) = cm + bm
			zout(1,nout7,j) = am + dm
			zout(2,nout7,j) = cm - bm
			r= r1 - r5
			s=-s3 + s7
			ap=r + s
			am=r - s
			r=s1 - s5
			s=r7 - r3
			bp=r + s
			bm=r - s
			r=-s4 + s8
			s= r2 - r6
			cp=r + s
			cm=r - s
			r=-s2 + s6
			s= r4 - r8
			dp=r + s
			dm=r - s
			r = ( cp + dm)*rt2i
			s = ( cp - dm)*rt2i
			cp= ( cm + dp)*rt2i
			dp= (-cm + dp)*rt2i
			zout(1,nout2,j) = ap + r
			zout(2,nout2,j) = bm + s
			zout(1,nout6,j) = ap - r
			zout(2,nout6,j) = bm - s
			zout(1,nout4,j) = am + cp
			zout(2,nout4,j) = bp + dp
			zout(1,nout8,j) = am - cp
			zout(2,nout8,j) = bp - dp
8121			continue
	endif
	else if (now.eq.3) then 
c	 .5d0*sqrt(3.d0)
	bb=isign*0.8660254037844387d0
	ia=1
	nin1=ia-after
	nout1=ia-atn
	do 3001,ib=1,before
	nin1=nin1+after
	nin2=nin1+atb
	nin3=nin2+atb
	nout1=nout1+atn
	nout2=nout1+after
	nout3=nout2+after
	do 3001,j=1,nfft
	r1=zin(1,j,nin1)
	s1=zin(2,j,nin1)
	r2=zin(1,j,nin2)
	s2=zin(2,j,nin2)
	r3=zin(1,j,nin3)
	s3=zin(2,j,nin3)
	r=r2 + r3
	s=s2 + s3
	zout(1,nout1,j) = r + r1
	zout(2,nout1,j) = s + s1
	r1=r1 - .5d0*r
	s1=s1 - .5d0*s
	r2=bb*(r2-r3)
	s2=bb*(s2-s3)
	zout(1,nout2,j) = r1 - s2 
	zout(2,nout2,j) = s1 + r2
	zout(1,nout3,j) = r1 + s2 
	zout(2,nout3,j) = s1 - r2
3001	continue
	do 3000,ia=2,after
	ias=ia-1
	if (4*ias.eq.3*after) then
	if (isign.eq.1) then
		nin1=ia-after
		nout1=ia-atn
		do 3010,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		do 3010,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r2=-zin(2,j,nin2)
		s2=zin(1,j,nin2)
		r3=-zin(1,j,nin3)
		s3=-zin(2,j,nin3)
		r=r2 + r3
		s=s2 + s3
		zout(1,nout1,j) = r + r1
		zout(2,nout1,j) = s + s1
		r1=r1 - .5d0*r
		s1=s1 - .5d0*s	
		r2=bb*(r2-r3)	
		s2=bb*(s2-s3)
		zout(1,nout2,j) = r1 - s2 
		zout(2,nout2,j) = s1 + r2
		zout(1,nout3,j) = r1 + s2 
		zout(2,nout3,j) = s1 - r2
3010		continue
	else
		nin1=ia-after
		nout1=ia-atn
		do 3020,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		do 3020,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r2=zin(2,j,nin2)
		s2=-zin(1,j,nin2)
		r3=-zin(1,j,nin3)
		s3=-zin(2,j,nin3)
		r=r2 + r3
		s=s2 + s3
		zout(1,nout1,j) = r + r1
		zout(2,nout1,j) = s + s1
		r1=r1 - .5d0*r
		s1=s1 - .5d0*s	
		r2=bb*(r2-r3)	
		s2=bb*(s2-s3)
		zout(1,nout2,j) = r1 - s2 
		zout(2,nout2,j) = s1 + r2
		zout(1,nout3,j) = r1 + s2 
		zout(2,nout3,j) = s1 - r2
3020		continue
	endif
	else if (8*ias.eq.3*after) then
	if (isign.eq.1) then
		nin1=ia-after
		nout1=ia-atn
		do 3030,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		do 3030,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r=zin(1,j,nin2)
		s=zin(2,j,nin2)
		r2=(r - s)*rt2i
		s2=(r + s)*rt2i
		r3=-zin(2,j,nin3)
		s3=zin(1,j,nin3) 
		r=r2 + r3
		s=s2 + s3
		zout(1,nout1,j) = r + r1
		zout(2,nout1,j) = s + s1
		r1=r1 - .5d0*r
		s1=s1 - .5d0*s	
		r2=bb*(r2-r3)	
		s2=bb*(s2-s3)
		zout(1,nout2,j) = r1 - s2 
		zout(2,nout2,j) = s1 + r2
		zout(1,nout3,j) = r1 + s2 
		zout(2,nout3,j) = s1 - r2
3030		continue
	else
		nin1=ia-after
		nout1=ia-atn
		do 3040,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		do 3040,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r=zin(1,j,nin2)
		s=zin(2,j,nin2)
		r2=(r + s)*rt2i
		s2=(-r + s)*rt2i
		r3=zin(2,j,nin3)
		s3=-zin(1,j,nin3)
		r=r2 + r3
		s=s2 + s3
		zout(1,nout1,j) = r + r1
		zout(2,nout1,j) = s + s1
		r1=r1 - .5d0*r
		s1=s1 - .5d0*s	
		r2=bb*(r2-r3)	
		s2=bb*(s2-s3)
		zout(1,nout2,j) = r1 - s2 
		zout(2,nout2,j) = s1 + r2
		zout(1,nout3,j) = r1 + s2 
		zout(2,nout3,j) = s1 - r2
3040		continue
	endif
	else
	itt=ias*before
	itrig=itt+1
	cr2=trig(1,itrig)
	ci2=trig(2,itrig)
	itrig=itrig+itt
	cr3=trig(1,itrig)
	ci3=trig(2,itrig)
	nin1=ia-after
	nout1=ia-atn
	do 3090,ib=1,before
	nin1=nin1+after
	nin2=nin1+atb
	nin3=nin2+atb
	nout1=nout1+atn
	nout2=nout1+after
	nout3=nout2+after
	do 3090,j=1,nfft
	r1=zin(1,j,nin1)
	s1=zin(2,j,nin1)
	r=zin(1,j,nin2)
	s=zin(2,j,nin2)
	r2=r*cr2 - s*ci2
	s2=r*ci2 + s*cr2
	r=zin(1,j,nin3)
	s=zin(2,j,nin3)
	r3=r*cr3 - s*ci3
	s3=r*ci3 + s*cr3
	r=r2 + r3
	s=s2 + s3
	zout(1,nout1,j) = r + r1
	zout(2,nout1,j) = s + s1
	r1=r1 - .5d0*r
	s1=s1 - .5d0*s
	r2=bb*(r2-r3)
	s2=bb*(s2-s3)
	zout(1,nout2,j) = r1 - s2 
	zout(2,nout2,j) = s1 + r2
	zout(1,nout3,j) = r1 + s2 
	zout(2,nout3,j) = s1 - r2
3090	continue
	endif
3000	continue
	else if (now.eq.5) then
c	 cos(2.d0*pi/5.d0)
	cos2=0.3090169943749474d0
c	 cos(4.d0*pi/5.d0)
	cos4=-0.8090169943749474d0
c	sin(2.d0*pi/5.d0)
	sin2=isign*0.9510565162951536d0
c	 sin(4.d0*pi/5.d0)
	sin4=isign*0.5877852522924731d0
	ia=1
	nin1=ia-after
	nout1=ia-atn
	do 5001,ib=1,before
	nin1=nin1+after
	nin2=nin1+atb
	nin3=nin2+atb
	nin4=nin3+atb
	nin5=nin4+atb
	nout1=nout1+atn
	nout2=nout1+after
	nout3=nout2+after
	nout4=nout3+after
	nout5=nout4+after
	do 5001,j=1,nfft
	r1=zin(1,j,nin1)
	s1=zin(2,j,nin1)
	r2=zin(1,j,nin2)
	s2=zin(2,j,nin2)
	r3=zin(1,j,nin3)
	s3=zin(2,j,nin3)
	r4=zin(1,j,nin4)
	s4=zin(2,j,nin4)
	r5=zin(1,j,nin5)
	s5=zin(2,j,nin5)
	r25 = r2 + r5
	r34 = r3 + r4
	s25 = s2 - s5
	s34 = s3 - s4
	zout(1,nout1,j) = r1 + r25 + r34
	r = cos2*r25 + cos4*r34 + r1
	s = sin2*s25 + sin4*s34
	zout(1,nout2,j) = r - s
	zout(1,nout5,j) = r + s
	r = cos4*r25 + cos2*r34 + r1
	s = sin4*s25 - sin2*s34
	zout(1,nout3,j) = r - s
	zout(1,nout4,j) = r + s
	r25 = r2 - r5
	r34 = r3 - r4
	s25 = s2 + s5
	s34 = s3 + s4
	zout(2,nout1,j) = s1 + s25 + s34
	r = cos2*s25 + cos4*s34 + s1
	s = sin2*r25 + sin4*r34
	zout(2,nout2,j) = r + s
	zout(2,nout5,j) = r - s
	r = cos4*s25 + cos2*s34 + s1
	s = sin4*r25 - sin2*r34
	zout(2,nout3,j) = r + s
	zout(2,nout4,j) = r - s
5001	continue
	do 5000,ia=2,after
	ias=ia-1
	if (8*ias.eq.5*after) then
		if (isign.eq.1) then
			nin1=ia-after
			nout1=ia-atn
			do 5010,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nin5=nin4+atb	
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			nout5=nout4+after
			do 5010,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=(r - s)*rt2i
			s2=(r + s)*rt2i
			r3=-zin(2,j,nin3)
			s3=zin(1,j,nin3) 
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=-(r + s)*rt2i
			s4=(r - s)*rt2i
			r5=-zin(1,j,nin5)
			s5=-zin(2,j,nin5)
			r25 = r2 + r5
			r34 = r3 + r4
			s25 = s2 - s5
			s34 = s3 - s4
			zout(1,nout1,j) = r1 + r25 + r34
			r = cos2*r25 + cos4*r34 + r1
			s = sin2*s25 + sin4*s34
			zout(1,nout2,j) = r - s
			zout(1,nout5,j) = r + s
			r = cos4*r25 + cos2*r34 + r1
			s = sin4*s25 - sin2*s34
			zout(1,nout3,j) = r - s
			zout(1,nout4,j) = r + s
			r25 = r2 - r5
			r34 = r3 - r4
			s25 = s2 + s5
			s34 = s3 + s4
			zout(2,nout1,j) = s1 + s25 + s34
			r = cos2*s25 + cos4*s34 + s1
			s = sin2*r25 + sin4*r34
			zout(2,nout2,j) = r + s
			zout(2,nout5,j) = r - s
			r = cos4*s25 + cos2*s34 + s1
			s = sin4*r25 - sin2*r34
			zout(2,nout3,j) = r + s
			zout(2,nout4,j) = r - s
5010			continue
		else
			nin1=ia-after
			nout1=ia-atn
			do 5020,ib=1,before
			nin1=nin1+after
			nin2=nin1+atb
			nin3=nin2+atb
			nin4=nin3+atb
			nin5=nin4+atb	
			nout1=nout1+atn
			nout2=nout1+after
			nout3=nout2+after
			nout4=nout3+after
			nout5=nout4+after
			do 5020,j=1,nfft
			r1=zin(1,j,nin1)
			s1=zin(2,j,nin1)
			r=zin(1,j,nin2)
			s=zin(2,j,nin2)
			r2=(r + s)*rt2i
			s2=(-r + s)*rt2i
			r3=zin(2,j,nin3)
			s3=-zin(1,j,nin3)
			r=zin(1,j,nin4)
			s=zin(2,j,nin4)
			r4=(s - r)*rt2i
			s4=-(r + s)*rt2i
			r5=-zin(1,j,nin5)
			s5=-zin(2,j,nin5)
			r25 = r2 + r5
			r34 = r3 + r4
			s25 = s2 - s5
			s34 = s3 - s4
			zout(1,nout1,j) = r1 + r25 + r34
			r = cos2*r25 + cos4*r34 + r1
			s = sin2*s25 + sin4*s34
			zout(1,nout2,j) = r - s
			zout(1,nout5,j) = r + s
			r = cos4*r25 + cos2*r34 + r1
			s = sin4*s25 - sin2*s34
			zout(1,nout3,j) = r - s
			zout(1,nout4,j) = r + s
			r25 = r2 - r5
			r34 = r3 - r4
			s25 = s2 + s5
			s34 = s3 + s4
			zout(2,nout1,j) = s1 + s25 + s34
			r = cos2*s25 + cos4*s34 + s1
			s = sin2*r25 + sin4*r34
			zout(2,nout2,j) = r + s
			zout(2,nout5,j) = r - s
			r = cos4*s25 + cos2*s34 + s1
			s = sin4*r25 - sin2*r34
			zout(2,nout3,j) = r + s
			zout(2,nout4,j) = r - s
5020			continue
		endif
	else
		ias=ia-1
		itt=ias*before
		itrig=itt+1
		cr2=trig(1,itrig)
		ci2=trig(2,itrig)
		itrig=itrig+itt
		cr3=trig(1,itrig)
		ci3=trig(2,itrig)
		itrig=itrig+itt
		cr4=trig(1,itrig)
		ci4=trig(2,itrig)
		itrig=itrig+itt
		cr5=trig(1,itrig)
		ci5=trig(2,itrig)
		nin1=ia-after
		nout1=ia-atn
		do 5100,ib=1,before
		nin1=nin1+after
		nin2=nin1+atb
		nin3=nin2+atb
		nin4=nin3+atb
		nin5=nin4+atb
		nout1=nout1+atn
		nout2=nout1+after
		nout3=nout2+after
		nout4=nout3+after
		nout5=nout4+after
		do 5100,j=1,nfft
		r1=zin(1,j,nin1)
		s1=zin(2,j,nin1)
		r=zin(1,j,nin2)
		s=zin(2,j,nin2)
		r2=r*cr2 - s*ci2
		s2=r*ci2 + s*cr2
		r=zin(1,j,nin3)
		s=zin(2,j,nin3)
		r3=r*cr3 - s*ci3
		s3=r*ci3 + s*cr3
		r=zin(1,j,nin4)
		s=zin(2,j,nin4)
		r4=r*cr4 - s*ci4
		s4=r*ci4 + s*cr4
		r=zin(1,j,nin5)
		s=zin(2,j,nin5)
		r5=r*cr5 - s*ci5
		s5=r*ci5 + s*cr5
		r25 = r2 + r5
		r34 = r3 + r4
		s25 = s2 - s5
		s34 = s3 - s4
		zout(1,nout1,j) = r1 + r25 + r34
		r = cos2*r25 + cos4*r34 + r1
		s = sin2*s25 + sin4*s34
		zout(1,nout2,j) = r - s
		zout(1,nout5,j) = r + s
		r = cos4*r25 + cos2*r34 + r1
		s = sin4*s25 - sin2*s34
		zout(1,nout3,j) = r - s
		zout(1,nout4,j) = r + s
		r25 = r2 - r5
		r34 = r3 - r4
		s25 = s2 + s5
		s34 = s3 + s4
		zout(2,nout1,j) = s1 + s25 + s34
		r = cos2*s25 + cos4*s34 + s1
		s = sin2*r25 + sin4*r34
		zout(2,nout2,j) = r + s
		zout(2,nout5,j) = r - s
		r = cos4*s25 + cos2*s34 + s1
		s = sin4*r25 - sin2*r34
		zout(2,nout3,j) = r + s
		zout(2,nout4,j) = r - s
5100		continue
	endif
5000	continue
       else if (now.eq.6) then
c	 .5d0*sqrt(3.d0)
	bb=isign*0.8660254037844387d0

	ia=1
	nin1=ia-after
	nout1=ia-atn
	do 6120,ib=1,before
	nin1=nin1+after
	nin2=nin1+atb
	nin3=nin2+atb
	nin4=nin3+atb
	nin5=nin4+atb
	nin6=nin5+atb
	nout1=nout1+atn
	nout2=nout1+after
	nout3=nout2+after
	nout4=nout3+after
	nout5=nout4+after
	nout6=nout5+after
	do 6120,j=1,nfft
	r2=zin(1,j,nin3)
	s2=zin(2,j,nin3)
	r3=zin(1,j,nin5)
	s3=zin(2,j,nin5)
	r=r2 + r3
	s=s2 + s3
	r1=zin(1,j,nin1)
	s1=zin(2,j,nin1)
	ur1 = r + r1
	ui1 = s + s1
	r1=r1 - .5d0*r
	s1=s1 - .5d0*s
	r=r2-r3
	s=s2-s3
	ur2 = r1 - s*bb
	ui2 = s1 + r*bb
	ur3 = r1 + s*bb
	ui3 = s1 - r*bb

	r2=zin(1,j,nin6)
	s2=zin(2,j,nin6)
	r3=zin(1,j,nin2)
	s3=zin(2,j,nin2)
	r=r2 + r3
	s=s2 + s3
	r1=zin(1,j,nin4)
	s1=zin(2,j,nin4)
	vr1 = r + r1
	vi1 = s + s1
	r1=r1 - .5d0*r
	s1=s1 - .5d0*s
	r=r2-r3
	s=s2-s3
	vr2 = r1 - s*bb
	vi2 = s1 + r*bb
	vr3 = r1 + s*bb
	vi3 = s1 - r*bb

	zout(1,nout1,j)=ur1+vr1
	zout(2,nout1,j)=ui1+vi1
	zout(1,nout5,j)=ur2+vr2
	zout(2,nout5,j)=ui2+vi2
	zout(1,nout3,j)=ur3+vr3
	zout(2,nout3,j)=ui3+vi3
	zout(1,nout4,j)=ur1-vr1
	zout(2,nout4,j)=ui1-vi1
	zout(1,nout2,j)=ur2-vr2
	zout(2,nout2,j)=ui2-vi2
	zout(1,nout6,j)=ur3-vr3
	zout(2,nout6,j)=ui3-vi3
6120	continue

       else
	stop 'error fftrot'
       endif

	return
	end