
!! Copyright (C) 2002-2007 BigDFT group 
!! This file is distributed under the terms of the
!! GNU General Public License, see ~/COPYING file
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the list of contributors, see ~/AUTHORS 



! --------------------------------------------------------------
!   3-dimensional complex-complex FFT routine: 
!   When compared to the best vendor implementations on RISC architectures 
!   it gives close to optimal performance (perhaps loosing 20 percent in speed)
!   and it is significanly faster than many not so good vendor implementations 
!   as well as other portable FFT's. 
!   On all vector machines tested so far (Cray, NEC, Fujitsu) is 
!   was significantly faster than the vendor routines
! The theoretical background is described in :
! 1) S. Goedecker: Rotating a three-dimensional array in optimal
! positions for vector processing: Case study for a three-dimensional Fast
! Fourier Transform, Comp. Phys. Commun. \underline{76}, 294 (1993)
! Citing of this reference is greatly appreciated if the routines are used 
! for scientific work.
! Is parallelized with OpenMP

!-----------------------------------------------------------

!  HEADER PART: ntime=1: Checks accuracy 
!               ntime<>1: Measures speed

        Implicit real*8 (a-h,o-z)
        integer count1,count2,count_rate,count_max

        REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: z  
        REAL*8, ALLOCATABLE, DIMENSION(:,:) :: zin

        write(*,*) 'n1,n2,n3=?'
        read(*,*) n1,n2,n3
        write(*,*) 'ntime=?'
        read(*,*) ntime
!        ntime=1
!        n2=6 ; n3=7
!        n1=3
!do 
!        call fourier_dim(n1,n_next)
!        n1=n_next

        if (mod(n1,2).eq.0) then 
          nd1=n1+1
        else
          nd1=n1
        endif
        if (mod(n2,2).eq.0) then 
          nd2=n2+1
        else
          nd2=n2
        endif
        if (mod(n3,2).eq.0) then 
          nd3=n3+1
        else
          nd3=n3
        endif

! general array
        allocate(zin(2,n1*n2*n3))
! arrays for FFT 
        allocate(z(2,nd1*nd2*nd3,2))


        Print*,'                    '
        Print*,'                    '
        write(6,'(a,3(i7),3x,i3)') 'FFT TEST for n1 ,n2 ,n3 =',n1,n2,n3
        write(6,*) 'nd1 ,nd2 ,nd3=', nd1 ,nd2 ,nd3
        Print*,'                    '
        isign=1
        do 398,i=1,nd1*nd2*nd3
        z(1,i,1)=0.d0
        z(2,i,1)=0.d0
        z(1,i,2)=0.d0
        z(2,i,2)=0.d0
398        continue
1234        call init(n1,n2,n3,nd1,nd2,nd3,zin,z)

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

        call cpu_time(t1)
        call system_clock(count1,count_rate,count_max)      

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

        call cpu_time(t2)
        call system_clock(count2,count_rate,count_max)      
        time=(t2-t1)/ntime
        tela=(count2-count1)/(float(count_rate)*ntime)

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

        if (ntime.ne.1) write(6,'(a,2(x,e11.4),x,i4)')  & 
              'Time (CPU,ELA) per FFT call (sec):' ,time,tela,ntime
        flops=5*n1*n2*n3*log(1.d0*n1*n2*n3)/log(2.d0)
	write(6,*) 'Estimated floating point operations per FFT call',flops
        if (ntime.ne.1) print*, 'CPU:  Mflops:' ,1.d-6*flops/time
        if (ntime.ne.1) print*, 'ELAP: Mflops:' ,1.d-6*flops/tela

        if (isign.eq.-1) then
        goto 1234
        endif
        write(6,*) '   '
        write(6,*) 'The following output should be a real number and ',&
           'not a NANq. If this happens overflow occured because ',& 
           'the FFT was called too many times (ntime too big).:' 
        write(6,*) z(1,1,1)


!        if (n1.eq.24000) exit
!        n1=n1+1
        deallocate(z,zin)
!enddo

        end



! auxiliary subroutines --------------------------

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

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


! FFT PART -----------------------------------------------------------------

        subroutine FFT(n1,n2,n3,nd1,nd2,nd3,z,isign,inzee)
!        CALCULATES THE DISCRETE FOURIERTRANSFORM F(I1,I2,I3)=
!        S_(j1,j2,j3) EXP(isign*i*2*pi*(j1*i1/n1+j2*i2/n2+j3*i3/n3)) R(j1,j2,j3)
!       with optimal performance on vector computer, workstations and 
!       multiprocessor shared memory computers using OpenMP compiler directives
!        INPUT:
!            n1,n2,n3:physical dimension of the transform. It must be a 
!                     product of the prime factors 2,3,5, but greater than 3. 
!                     One factor of 7 is also allowed.  
!                    If two ni's are equal it is recommended to place them 
!                    behind each other.
!            nd1,nd2,nd3:memory dimension of Z. ndi must always be greater or 
!                        equal than ni. On a vector machine, it is recomended 
!                       to chose ndi=ni if ni is odd and ndi=ni+1 if ni is 
!                       even to obtain optimal execution speed. On RISC 
!                       machines ndi=ni is usually fine for odd ni, for even 
!                       ni one should try ndi=ni+1, ni+2, ni+4 to find the 
!                       optimal performance. 
!           inzee=1: first part of Z is data (input) array, 
!                    second part work array
!           inzee=2: first part of Z is work array, second part data array
!                Z(1,i1,i2,i3,inzee)=real(R(i1,i2,i3))
!                Z(2,i1,i2,i3,inzee)=imag(R(i1,i2,i3))
!        OUTPUT:
!           inzee=1: first part of Z is data (output) array, 
!                    second part work array
!           inzee=2: first part of Z is work array, second part data array
!                real(F(i1,i2,i3))=Z(1,i1,i2,i3,inzee)
!                imag(F(i1,i2,i3))=Z(2,i1,i2,i3,inzee)
!           inzee on output is in general different from inzee on input
!        The input data are always overwritten independently of the 
!       value of inzee.
! PERFORMANCE AND THE NCACHE
!       The most important feature for performance is the right choice of 
!       the parameter ncache. On a vector machine ncache has to be put to 0.
!       On a RISC machine with cache, it is very important to find the optimal 
!       value of NCACHE. NCACHE determines the size of the work array zw, that
!       has to fit into cache. It has therefore to be chosen to equal roughly 
!        half the size of the physical cache in units of real*8 numbers.
!       If the machine has 2 cache levels it can not be predicted which 
!       cache level will be the most relevant one for choosing ncache. 
!       The optimal value of ncache can easily be determined by numerical 
!       experimentation. A too large value of ncache leads to a dramatic 
!       and sudden decrease of performance, a too small value to a to a 
!       slow and less dramatic decrease of performance. If NCACHE is set 
!       to a value so small, that not even a single one dimensional transform 
!       can be done in the workarray zw, the program stops with an error 
!       message.
!  Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
!  Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1995, 1999
!  Copyright (C) Stefan Goedecker, CEA Grenoble, 2002
!  the part for radix 7 added by Alexey Neelov, Basel University,  2008
!  This file is distributed under the terms of the
!  GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .

        implicit real*8 (a-h,o-z)
!$      interface
!$        integer ( kind=4 ) function omp_get_num_threads ( )
!$        end function omp_get_num_threads
!$      end interface
!$      interface
!$        integer ( kind=4 ) function omp_get_thread_num ( )
!$        end function omp_get_thread_num
!$      end interface


        REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: zw  
        REAL*8, ALLOCATABLE, DIMENSION(:,:) :: trig
        INTEGER, ALLOCATABLE, DIMENSION(:) :: after,now,before
        dimension z(2,nd1*nd2*nd3,2)
        if (max(n1,n2,n3).gt.24000) stop ' dimension larger than 24000'

! some reasonable values of ncache: 
!   IBM/RS6000/590: 16*1024 ; IBM/RS6000/390: 3*1024 ; 
!   IBM/PwPC: 1*1024 ; SGI/MIPS/R8000: 16*1024 ; DEC/Alpha/EV5 and EV6 6*1024
!   But if you care about performance find the optimal value of ncache yourself!
!       On all vector machines: ncache=0

        ncache=6*1024

! check whether input values are reasonable
	if (inzee.le.0 .or. inzee.ge.3) stop 'wrong inzee'
	if (isign.ne.1 .and. isign.ne.-1) stop 'wrong isign'
	if (n1.gt.nd1) stop 'n1>nd1'
	if (n2.gt.nd2) stop 'n2>nd2'
	if (n3.gt.nd3) stop 'n3>nd3'
        ntrig=max(n1,n2,n3)
	

! vector computer with memory banks:
      if (ncache.eq.0) then
        allocate(trig(2,ntrig),after(20),now(20),before(20))

        call ctrig(n3,ntrig,trig,after,before,now,isign,ic)
        nfft=nd1*n2
        mm=nd1*nd2
        do 51093,i=1,ic-1
        call fftstp(mm,nfft,nd3,mm,nd3,z(1,1,inzee),z(1,1,3-inzee), &
                          ntrig,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), &
                          ntrig,trig,after(i),now(i),before(i),isign)
        inzee=3-inzee

        if (n2.ne.n3) call ctrig(n2,ntrig,trig,after,before,now,isign,ic)
        nfft=nd3*n1
        mm=nd3*nd1
        do 52093,i=1,ic-1
        call fftstp(mm,nfft,nd2,mm,nd2,z(1,1,inzee),z(1,1,3-inzee), &
                           ntrig,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), &
                       ntrig,trig,after(i),now(i),before(i),isign)
        inzee=3-inzee

        if (n1.ne.n2) call ctrig(n1,ntrig,trig,after,before,now,isign,ic)
        nfft=nd2*n3
        mm=nd2*nd3
        do 53093,i=1,ic-1
        call fftstp(mm,nfft,nd1,mm,nd1,z(1,1,inzee),z(1,1,3-inzee), &
                         ntrig,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), &
                         ntrig,trig,after(i),now(i),before(i),isign)
        inzee=3-inzee

! RISC machine with cache:
      else
! INtel IFC does not understand default(private)
!!$omp parallel  default(private) &
!$omp parallel & 
!$omp private(zw,trig,before,after,now,i,j,iam,npr,jj,ma,mb,mm,ic,n,m,jompa,jompb,lot,lotomp,inzeep,inzet,nn,nfft) &
!$omp shared(n1,n2,n3,ntrig,nd1,nd2,nd3,z,isign,inzee,ncache) 
        npr=1
!$       npr=omp_get_num_threads()
        iam=0
!$       iam=omp_get_thread_num()
!        write(6,*) 'npr,iam',npr,iam
! Critical section only necessary on Intel
!$omp critical
        allocate(zw(2,ncache/4,2),trig(2,ntrig),after(20),now(20),before(20))
!$omp end critical

        inzet=inzee
! TRANSFORM ALONG Z AXIS

        mm=nd1*nd2
        m=nd3
        lot=max(1,ncache/(4*n3))
        nn=lot
        n=n3
        if (2*n*lot*2.gt.ncache) stop 'ncache has to be larger than 4*n3'

        call ctrig(n3,ntrig,trig,after,before,now,isign,ic)

      if (ic.eq.1) then
        i=ic
        lotomp=(nd1*n2)/npr+1
        ma=iam*lotomp+1
        mb=min((iam+1)*lotomp,nd1*n2)
        nfft=mb-ma+1
        j=ma
        jj=j*nd3-nd3+1
        call fftrot(mm,nfft,m,mm,m,z(1,j,inzet),z(1,jj,3-inzet), &
                          ntrig,trig,after(i),now(i),before(i),isign)

      else

        lotomp=(nd1*n2)/npr+1
        jompa=iam*lotomp+1
        jompb=min((iam+1)*lotomp,nd1*n2)
        do 1000,j=jompa,jompb,lot
        ma=j
        mb=min(j+(lot-1),jompb)
        nfft=mb-ma+1
        jj=j*nd3-nd3+1

        i=1
        inzeep=2
        call fftstp(mm,nfft,m,nn,n,z(1,j,inzet),zw(1,1,3-inzeep), &
                         ntrig,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), &
                         ntrig,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-inzet), &
                         ntrig,trig,after(i),now(i),before(i),isign)
1000        continue
      endif

        inzet=3-inzet

!$omp barrier

! TRANSFORM ALONG Y AXIS
        mm=nd3*nd1
        m=nd2
        lot=max(1,ncache/(4*n2))
        nn=lot
        n=n2
        if (2*n*lot*2.gt.ncache) stop 'ncache has to be larger than 4*n2'

        if (n2.ne.n3) call ctrig(n2,ntrig,trig,after,before,now,isign,ic)

      if (ic.eq.1) then
        i=ic
        lotomp=(nd3*n1)/npr+1
        ma=iam*lotomp+1
        mb=min((iam+1)*lotomp,nd3*n1)
        nfft=mb-ma+1
        j=ma
        jj=j*nd2-nd2+1
        call fftrot(mm,nfft,m,mm,m,z(1,j,inzet),z(1,jj,3-inzet), &
                         ntrig,trig,after(i),now(i),before(i),isign)

      else

        lotomp=(nd3*n1)/npr+1
        jompa=iam*lotomp+1
        jompb=min((iam+1)*lotomp,nd3*n1)
        do 2000,j=jompa,jompb,lot
        ma=j
        mb=min(j+(lot-1),jompb)
        nfft=mb-ma+1
        jj=j*nd2-nd2+1

        i=1
        inzeep=2
        call fftstp(mm,nfft,m,nn,n,z(1,j,inzet),zw(1,1,3-inzeep), &
                         ntrig,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), &
                         ntrig,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-inzet), &
                         ntrig,trig,after(i),now(i),before(i),isign)
2000        continue
      endif
        inzet=3-inzet

!$omp barrier

! TRANSFORM ALONG X AXIS
        mm=nd2*nd3
        m=nd1
        lot=max(1,ncache/(4*n1))
        nn=lot
        n=n1
        if (2*n*lot*2.gt.ncache) stop 'ncache has to be larger than 4*n1'

        if (n1.ne.n2) call ctrig(n1,ntrig,trig,after,before,now,isign,ic)

      if (ic.eq.1) then
        i=ic
        lotomp=(nd2*n3)/npr+1
        ma=iam*lotomp+1
        mb=min((iam+1)*lotomp,nd2*n3)
        nfft=mb-ma+1
        j=ma
        jj=j*nd1-nd1+1
        call fftrot(mm,nfft,m,mm,m,z(1,j,inzet),z(1,jj,3-inzet), &
                         ntrig,trig,after(i),now(i),before(i),isign)

      else

        lotomp=(nd2*n3)/npr+1
        jompa=iam*lotomp+1
        jompb=min((iam+1)*lotomp,nd2*n3)
        do 3000,j=jompa,jompb,lot
        ma=j
        mb=min(j+(lot-1),jompb)
        nfft=mb-ma+1
        jj=j*nd1-nd1+1

        i=1
        inzeep=2
        call fftstp(mm,nfft,m,nn,n,z(1,j,inzet),zw(1,1,3-inzeep), &
                         ntrig,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), &
                         ntrig,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-inzet), &
                         ntrig,trig,after(i),now(i),before(i),isign)
3000        continue
      endif
        inzet=3-inzet
        
        deallocate(zw,trig,after,now,before)
        if (iam.eq.0) inzee=inzet
!$omp end parallel  


      endif
        return
        end


!!****f* BigDFT/fourier_dim
!! NAME
!!   fourier_dim
!!
!! FUNCTION
!!   Give a number n_next > n compatible for the FFT
!!
!! SOURCE
!!
!!***
!  Copyright (C) Stefan Goedecker, CEA Grenoble, 2002
!  This file is distributed under the terms of the
!  GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
subroutine fourier_dim(n,n_next)
  implicit none
  !Arguments
  integer  :: n
  integer :: n_next

  !Local variables
  integer, parameter :: ndata = 286
  !Multiple of 2,3,5
  integer, dimension(ndata), parameter :: idata = (/   &
       3,      4,      5,      6,      7,      8,      9,     12,  &
      14,     15,     16,     18,     20,     21,     24,     25,  &
      27,     28,     30,     32,     35,     36,     40,     42,  &
      45,     48,     50,     54,     56,     60,     63,     64,  &
      70,     72,     75,     80,     81,     84,     90,     96,  &
     100,    105,    108,    112,    120,    125,    126,    128,  &
     135,    140,    144,    150,    160,    162,    168,    175,  &
     180,    189,    192,    200,    210,    216,    224,    225,  &
     240,    243,    252,    256,    270,    280,    288,    300,  &
     315,    320,    324,    336,    360,    375,    378,    384,  &
     400,    405,    420,    432,    448,    450,    480,    486,  &
     500,    504,    512,    525,    540,    560,    567,    576,  &
     600,    625,    630,    640,    648,    672,    675,    700,  &
     720,    729,    750,    756,    768,    800,    810,    840,  &
     864,    875,    896,    900,    945,    960,    972,   1000,  &
    1008,   1024,   1050,   1080,   1120,   1125,   1134,   1152,  &
    1200,   1215,   1260,   1280,   1296,   1344,   1350,   1400,  &
    1440,   1458,   1500,   1512,   1536,   1575,   1600,   1620,  &
    1680,   1701,   1728,   1792,   1800,   1875,   1890,   1920,  &
    1944,   2000,   2016,   2025,   2048,   2100,   2160,   2240,  &
    2250,   2268,   2304,   2400,   2430,   2500,   2520,   2560,  &
    2592,   2625,   2688,   2700,   2800,   2835,   2880,   3000,  &
    3024,   3072,   3125,   3150,   3200,   3240,   3360,   3375,  &
            3456,   3500,   3584,   3600,   3750,   3780,   3840,  &
    3888,   4000,   4032,   4050,   4096,   4200,   4320,   4375,  &
    4480,   4500,   4536,   4608,   4725,   4800,   5000,   5040,  &
    5120,   5184,   5250,   5376,   5400,   5600,   5625,          &
    5760,   6000,   6048,   6144,   6300,   6400,   6480,   6720,  &
    6750,   6912,   7000,   7168,   7200,   7500,   7560,   7680,  &
    7875,   8000,   8064,   8192,   8400,   8640,   8960,   9000,  &
            9216,   9375,            600,  10000,  10080,  10240,  &
   10368,  10500,  10752,  10800,  11200,  11250,  11520,  12000,  &
   12096,  12288,  12500,  12600,  12800,  13125,  13440,  13824,  &
   14000,  14336,  14400,  15000,          15360,  15625,          &
   16000,  16128,  16384,  16800,  17280,  17500,  17920,  18000,  &
   18432,  18750,  19200,  20000,  20160,  20480,  21000,  21504,  &
   21875,  22400,  23040,  24000 /)
  integer :: i

  loop_data: do i=1,ndata
     if (n <= idata(i)) then
        n_next = idata(i)
        return
     end if
  end do loop_data
  write(unit=*,fmt=*) "fourier_dim: ",n," is bigger than ",idata(ndata)
  stop
end subroutine fourier_dim


        subroutine ctrig(n,ntrig,trig,after,before,now,isign,ic)
!  Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
!  Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
!  Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
!  This file is distributed under the terms of the
!  GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .

!     Different factorizations affect the performance
!     Factoring 64 as 4*4*4 might for example be faster on some machines than 8*8.

        implicit real*8 (a-h,o-z)
        integer after,before
        dimension now(7),after(7),before(7),trig(2,ntrig)
		integer,parameter::num_fac=286
        INTEGER, DIMENSION(7,num_fac) :: idata 
! The factor 6 is only allowed in the first place!
        data ((idata(i,j),i=1,7),j=1,num_fac) /                        &
       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,  &
       7,   7,   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,  &
      14,   7,   2,   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,  &
      21,   7,   3,   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,  &
      28,   7,   4,   1,   1,   1,   1,  &
      30,   6,   5,   1,   1,   1,   1,  &
      32,   8,   4,   1,   1,   1,   1,  &
      35,   7,   5,   1,   1,   1,   1,  &
      36,   4,   3,   3,   1,   1,   1,  &
      40,   8,   5,   1,   1,   1,   1,  &
      42,   7,   3,   2,   1,   1,   1,  &
      45,   5,   3,   3,   1,   1,   1,  &
      48,   4,   4,   3,   1,   1,   1,  &
      50,   5,   5,   2,   1,   1,   1,  &
      54,   6,   3,   3,   1,   1,   1,  &
      56,   7,   8,   1,   1,   1,   1,  &
      60,   5,   4,   3,   1,   1,   1,  &
      63,   7,   3,   3,   1,   1,   1,  &
      64,   8,   8,   1,   1,   1,   1,  &
      70,   7,   5,   2,   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,  &
      84,   7,   4,   3,   1,   1,   1,  &
      90,   6,   5,   3,   1,   1,   1,  &
      96,   8,   4,   3,   1,   1,   1,  &
     100,   5,   5,   4,   1,   1,   1,  &
     105,   7,   5,   3,   1,   1,   1,  &
     108,   4,   3,   3,   3,   1,   1,  &
     112,   7,   4,   4,   1,   1,   1,  &
     120,   8,   5,   3,   1,   1,   1,  &
     125,   5,   5,   5,   1,   1,   1,  &
     126,   7,   3,   3,   2,   1,   1,  &
     128,   8,   4,   4,   1,   1,   1,  &
     135,   5,   3,   3,   3,   1,   1,  &
     140,   7,   5,   4,   1,   1,   1,  &
     144,   6,   8,   3,   1,   1,   1,  &
     150,   6,   5,   5,   1,   1,   1,  &
     160,   8,   5,   4,   1,   1,   1,  &
     162,   6,   3,   3,   3,   1,   1,  &
     168,   7,   8,   3,   1,   1,   1,  &
     175,   7,   5,   5,   1,   1,   1,  &
     180,   5,   4,   3,   3,   1,   1,  &
     189,   7,   3,   3,   3,   1,   1,  &
     192,   6,   8,   4,   1,   1,   1,  &
     200,   8,   5,   5,   1,   1,   1,  &
     210,   7,   3,   5,   2,   1,   1,  &
     216,   8,   3,   3,   3,   1,   1,  &
     224,   7,   8,   4,   1,   1,   1,  &
     225,   5,   5,   3,   3,   1,   1,  &
     240,   6,   8,   5,   1,   1,   1,  &
     243,   3,   3,   3,   3,   3,   1,  &
     252,   7,   4,   3,   3,   1,   1,  &
     256,   8,   8,   4,   1,   1,   1,  &
     270,   6,   5,   3,   3,   1,   1,  &
     280,   7,   8,   5,   1,   1,   1,  &
     288,   8,   4,   3,   3,   1,   1,  &
     300,   5,   5,   4,   3,   1,   1,  &
     315,   7,   5,   3,   3,   1,   1,  &
     320,   5,   4,   4,   4,   1,   1,  &
     324,   4,   3,   3,   3,   3,   1,  &
     336,   7,   4,   4,   3,   1,   1,  &
     360,   8,   5,   3,   3,   1,   1,  &
     375,   5,   5,   5,   3,   1,   1,  &
     378,   7,   3,   3,   3,   2,   1,  &
     384,   8,   4,   4,   3,   1,   1,  &
     400,   5,   5,   4,   4,   1,   1,  &
     405,   5,   3,   3,   3,   3,   1,  &
     420,   7,   5,   4,   3,   1,   1,  &
     432,   4,   4,   3,   3,   3,   1,  &
     448,   7,   8,   8,   1,   1,   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,  &
     504,   7,   8,   3,   3,   1,   1,  &
     512,   8,   8,   8,   1,   1,   1,  &
     525,   7,   5,   5,   3,   1,   1,  &
     540,   5,   4,   3,   3,   3,   1,  &
     560,   7,   5,   4,   4,   1,   1,  &
     567,   7,   3,   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,  &
     630,   7,   3,   5,   3,   2,   1,  &
     640,   8,   5,   4,   4,   1,   1,  &
     648,   8,   3,   3,   3,   3,   1,  &
     672,   7,   8,   4,   3,   1,   1,  &
     675,   5,   5,   3,   3,   3,   1,  &
     700,   7,   5,   5,   4,   1,   1,  &
     720,   5,   4,   4,   3,   3,   1,  &
     729,   3,   3,   3,   3,   3,   3,  &
     750,   6,   5,   5,   5,   1,   1,  &
     756,   7,   4,   3,   3,   3,   1,  &
     768,   4,   4,   4,   4,   3,   1,  &
     800,   8,   5,   5,   4,   1,   1,  &
     810,   6,   5,   3,   3,   3,   1,  &
     840,   7,   8,   5,   3,   1,   1,  &
     864,   8,   4,   3,   3,   3,   1,  &
     875,   7,   5,   5,   5,   1,   1,  &
     896,   7,   8,   4,   4,   1,   1,  &
     900,   5,   5,   4,   3,   3,   1,  &
     945,   7,   5,   3,   3,   3,   1,  &
     960,   5,   4,   4,   4,   3,   1,  &
     972,   4,   3,   3,   3,   3,   3,  &
    1000,   8,   5,   5,   5,   1,   1,  &
    1008,   7,   8,   3,   3,   2,   1,  &
    1024,   4,   4,   4,   4,   4,   1,  &
    1050,   7,   3,   5,   5,   2,   1,  &
    1080,   6,   5,   4,   3,   3,   1,  &
    1120,   7,   8,   5,   4,   1,   1,  &
    1125,   5,   5,   5,   3,   3,   1,  &
    1134,   7,   3,   3,   3,   3,   2,  &
    1152,   6,   4,   4,   4,   3,   1,  &
    1200,   6,   8,   5,   5,   1,   1,  &
    1215,   5,   3,   3,   3,   3,   3,  &
    1260,   7,   5,   4,   3,   3,   1,  &
    1280,   8,   8,   5,   4,   1,   1,  &
    1296,   6,   8,   3,   3,   3,   1,  &
    1344,   7,   8,   3,   8,   1,   1,  &
    1350,   6,   5,   5,   3,   3,   1,  &
    1400,   7,   8,   5,   5,   1,   1,  &
    1440,   6,   5,   4,   4,   3,   1,  &
    1458,   6,   3,   3,   3,   3,   3,  &
    1500,   5,   5,   5,   4,   3,   1,  &
    1512,   7,   8,   3,   3,   3,   1,  &
    1536,   6,   8,   8,   4,   1,   1,  &
    1575,   7,   5,   5,   3,   3,   1,  &
    1600,   8,   8,   5,   5,   1,   1,  &
    1620,   5,   4,   3,   3,   3,   3,  &
    1680,   7,   8,   3,   5,   2,   1,  &
    1701,   7,   3,   3,   3,   3,   3,  &
    1728,   6,   8,   4,   3,   3,   1,  &
    1792,   7,   8,   8,   4,   1,   1,  &
    1800,   6,   5,   5,   4,   3,   1,  &
    1875,   5,   5,   5,   5,   3,   1,  &
    1890,   7,   3,   5,   3,   3,   2,  &
    1920,   6,   5,   4,   4,   4,   1,  &
    1944,   6,   4,   3,   3,   3,   3,  &
    2000,   5,   5,   5,   4,   4,   1,  &
    2016,   7,   8,   4,   3,   3,   1,  &
    2025,   5,   5,   3,   3,   3,   3,  &
    2048,   8,   4,   4,   4,   4,   1,  &
    2100,   7,   5,   5,   4,   3,   1,  &
    2160,   6,   8,   5,   3,   3,   1,  &
    2240,   7,   5,   4,   4,   4,   1,  &
    2250,   6,   5,   5,   5,   3,   1,  &
    2268,   7,   4,   3,   3,   3,   3,  &
    2304,   6,   8,   4,   4,   3,   1,  &
    2400,   6,   5,   5,   4,   4,   1,  &
    2430,   6,   5,   3,   3,   3,   3,  &
    2500,   5,   5,   5,   5,   4,   1,  &
    2520,   7,   8,   5,   3,   3,   1,  &
    2560,   8,   5,   4,   4,   4,   1,  &
    2592,   6,   4,   4,   3,   3,   3,  &
    2625,   7,   5,   5,   5,   3,   1,  &
    2688,   7,   8,   4,   4,   3,   1,  &
    2700,   5,   5,   4,   3,   3,   3,  &
    2800,   7,   5,   5,   4,   4,   1,  &
    2835,   7,   5,   3,   3,   3,   3,  &
    2880,   6,   8,   5,   4,   3,   1,  &
    3000,   6,   5,   5,   5,   4,   1,  &
    3024,   7,   4,   4,   3,   3,   3,  &
    3072,   6,   8,   4,   4,   4,   1,  &
    3125,   5,   5,   5,   5,   5,   1,  &
    3150,   7,   3,   5,   5,   3,   2,  &
    3200,   8,   5,   5,   4,   4,   1,  &
    3240,   6,   5,   4,   3,   3,   3,  &
    3360,   7,   8,   5,   4,   3,   1,  &
    3375,   5,   5,   5,   3,   3,   3,  &
    3456,   6,   4,   4,   4,   3,   3,  &
    3500,   7,   5,   5,   5,   4,   1,  &
    3584,   7,   8,   8,   8,   1,   1,  &
    3600,   6,   8,   5,   5,   3,   1,  &
    3750,   6,   5,   5,   5,   5,   1,  &
    3780,   7,   5,   4,   3,   3,   3,  &
    3840,   6,   8,   5,   4,   4,   1,  &
    3888,   6,   8,   3,   3,   3,   3,  &
    4000,   8,   5,   5,   5,   4,   1,  &
    4032,   7,   4,   4,   4,   3,   3,  &
    4050,   6,   5,   5,   3,   3,   3,  &
    4096,   8,   8,   4,   4,   4,   1,  &
    4200,   7,   8,   5,   5,   3,   1,  &
    4320,   6,   5,   4,   4,   3,   3,  &
    4375,   7,   5,   5,   5,   5,   1,  &
    4480,   7,   8,   5,   4,   4,   1,  &
    4500,   5,   5,   5,   4,   3,   3,  &
    4536,   7,   8,   3,   3,   3,   3,  &
    4608,   6,   8,   8,   4,   3,   1,  &
    4725,   7,   5,   5,   3,   3,   3,  &
    4800,   6,   8,   5,   5,   4,   1,  &
    5000,   8,   5,   5,   5,   5,   1,  &
    5040,   7,   5,   4,   4,   3,   3,  &
    5120,   8,   8,   5,   4,   4,   1,  &
    5184,   6,   8,   4,   3,   3,   3,  &
    5250,   7,   3,   5,   5,   5,   2,  &
    5376,   7,   4,   4,   4,   4,   3,  &
    5400,   6,   5,   5,   4,   3,   3,  &
    5600,   7,   8,   5,   5,   4,   1,  &
    5625,   5,   5,   5,   5,   3,   3,  &
    5760,   6,   8,   8,   5,   3,   1,  &
    6000,   6,   8,   5,   5,   5,   1,  &
    6048,   7,   8,   4,   3,   3,   3,  &
    6144,   6,   8,   8,   4,   4,   1,  &
    6300,   7,   5,   5,   4,   3,   3,  &
    6400,   8,   8,   5,   5,   4,   1,  &
    6480,   6,   8,   5,   3,   3,   3,  &
    6720,   7,   5,   4,   4,   4,   3,  &
    6750,   6,   5,   5,   5,   3,   3,  &
    6912,   6,   8,   4,   4,   3,   3,  &
    7000,   7,   8,   5,   5,   5,   1,  &
    7168,   7,   4,   4,   4,   4,   4,  &
    7200,   6,   5,   5,   4,   4,   3,  &
    7500,   5,   5,   5,   5,   4,   3,  &
    7560,   7,   3,   5,   8,   3,   3,  &
    7680,   6,   8,   8,   5,   4,   1,  &
    7875,   7,   5,   5,   5,   3,   3,  &
    8000,   8,   8,   5,   5,   5,   1,  &
    8064,   7,   3,   8,   4,   4,   3,  &
    8192,   8,   8,   8,   4,   4,   1,  &
    8400,   7,   8,   3,   5,   5,   2,  &
    8640,   8,   8,   5,   3,   3,   3,  &
    8960,   7,   8,   8,   5,   4,   1,  &
    9000,   8,   5,   5,   5,   3,   3,  &
    9216,   6,   8,   8,   8,   3,   1,  &
    9375,   5,   5,   5,   5,   5,   3,  &
    9600,   8,   5,   5,   4,   4,   3,  &
   10000,   5,   5,   5,   5,   4,   4,  &
   10080,   7,   3,   5,   8,   4,   3,  &
   10240,   8,   8,   8,   5,   4,   1,  &
   10368,   6,   8,   8,   3,   3,   3,  &
   10500,   7,   5,   5,   5,   4,   3,  &
   10752,   7,   8,   8,   3,   8,   1,  &
   10800,   6,   8,   5,   5,   3,   3,  &
   11200,   7,   8,   8,   5,   5,   1,  &
   11250,   6,   5,   5,   5,   5,   3,  &
   11520,   8,   8,   5,   4,   3,   3,  &
   12000,   8,   5,   5,   5,   4,   3,  &
   12096,   7,   8,   3,   8,   3,   3,  &
   12288,   8,   8,   8,   8,   3,   1,  &
   12500,   5,   5,   5,   5,   5,   4,  &
   12600,   7,   3,   5,   5,   8,   3,  &
   12800,   8,   8,   8,   5,   5,   1,  &
   13125,   7,   5,   5,   5,   5,   3,  &
   13440,   7,   3,   5,   8,   4,   4,  &
   13824,   8,   8,   8,   3,   3,   3,  &
   14000,   7,   5,   5,   5,   4,   4,  &
   14336,   7,   8,   4,   4,   4,   4,  &
   14400,   8,   8,   5,   5,   3,   3,  &
   15000,   8,   5,   5,   5,   5,   3,  &
   15360,   6,   8,   8,   8,   5,   1,  &
   15625,   5,   5,   5,   5,   5,   5,  &
   16000,   8,   5,   5,   5,   4,   4,  &
   16128,   7,   8,   3,   8,   4,   3,  &
   16384,   8,   8,   8,   8,   4,   1,  &
   16800,   7,   3,   5,   5,   8,   4,  &
   17280,   6,   8,   8,   5,   3,   3,  &
   17500,   7,   5,   5,   5,   5,   4,  &
   17920,   7,   8,   5,   4,   4,   4,  &
   18000,   6,   8,   5,   5,   5,   3,  &
   18432,   8,   8,   8,   4,   3,   3,  &
   18750,   6,   5,   5,   5,   5,   5,  &
   19200,   8,   8,   5,   5,   4,   3,  &
   20000,   8,   5,   5,   5,   5,   4,  &
   20160,   7,   8,   3,   5,   8,   3,  &
   20480,   8,   8,   8,   8,   5,   1,  &
   21000,   7,   3,   5,   5,   5,   8,  &
   21504,   7,   8,   3,   8,   4,   4,  &
   21875,   7,   5,   5,   5,   5,   5,  &
   22400,   7,   8,   5,   5,   4,   4,  &
   23040,   8,   8,   8,   5,   3,   3,  &
   24000,   8,   8,   5,   5,   5,   3   /  
        do 111,i=1,num_fac
        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,num_fac)
        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))
!        write(6,12) (after(i),i=1,ic)
!        write(6,12) (now(i),i=1,ic)
!        write(6,12) (before(i),i=1,ic)

        twopi=6.283185307179586d0
        angle=isign*twopi/n
        if (mod(n,2).eq.0) then
        nh=n/2
        trig(1,1)=1.d0
        trig(2,1)=0.d0
        trig(1,nh+1)=-1.d0
        trig(2,nh+1)=0.d0
        do 40,i=1,nh-1
        trigc=cos(i*angle)
        trigs=sin(i*angle)
        trig(1,i+1)=trigc
        trig(2,i+1)=trigs
        trig(1,n-i+1)=trigc
        trig(2,n-i+1)=-trigs
40      continue
        else
        nh=(n-1)/2
        trig(1,1)=1.d0
        trig(2,1)=0.d0
        do 20,i=1,nh
        trigc=cos(i*angle)
        trigs=sin(i*angle)
        trig(1,i+1)=trigc
        trig(2,i+1)=trigs
        trig(1,n-i+1)=trigc
        trig(2,n-i+1)=-trigs
20      continue
        endif


        return
        end


!ccccccccccccccccccccccccccccccccccccccccccccccc

         subroutine fftstp(mm,nfft,m,nn,n,zin,zout,ntrig,trig,after,now,before,isign)
!  Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
!  Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1995, 1999
!  This file is distributed under the terms of the
!  GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .

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

!         sqrt(.5d0)
        rt2i=0.7071067811865475d0
        if (now.eq.2) then
        ia=1
        nin1=ia-after
        nout1=ia-atn
        do 2001,ib=1,before
        nin1=nin1+after
        nin2=nin1+atb
        nout1=nout1+atn
        nout2=nout1+after
        do 2001,j=1,nfft
        r1=zin(1,j,nin1)
        s1=zin(2,j,nin1)
        r2=zin(1,j,nin2)
        s2=zin(2,j,nin2)
        zout(1,j,nout1)= r2 + r1
        zout(2,j,nout1)= s2 + s1
        zout(1,j,nout2)= r1 - r2
        zout(2,j,nout2)= s1 - s2
2001        continue
        do 2000,ia=2,after
        ias=ia-1
        if (2*ias.eq.after) then
                if (isign.eq.1) then
                        nin1=ia-after
                        nout1=ia-atn
                        do 2010,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2010,j=1,nfft
                        r1=zin(1,j,nin1)
                        s1=zin(2,j,nin1)
                        r2=zin(2,j,nin2)
                        s2=zin(1,j,nin2)
                        zout(1,j,nout1)= r1 - r2
                        zout(2,j,nout1)= s2 + s1
                        zout(1,j,nout2)= r2 + r1
                        zout(2,j,nout2)= s1 - s2
2010                        continue
                else
                        nin1=ia-after
                        nout1=ia-atn
                        do 2020,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2020,j=1,nfft
                        r1=zin(1,j,nin1)
                        s1=zin(2,j,nin1)
                        r2=zin(2,j,nin2)
                        s2=zin(1,j,nin2)
                        zout(1,j,nout1)= r2 + r1
                        zout(2,j,nout1)= s1 - s2
                        zout(1,j,nout2)= r1 - r2
                        zout(2,j,nout2)= s2 + s1
2020                        continue
                endif
        else if (4*ias.eq.after) then
                if (isign.eq.1) then
                        nin1=ia-after
                        nout1=ia-atn
                        do 2030,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2030,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
                        zout(1,j,nout1)= r2 + r1
                        zout(2,j,nout1)= s2 + s1
                        zout(1,j,nout2)= r1 - r2
                        zout(2,j,nout2)= s1 - s2
2030                        continue
                else
                        nin1=ia-after
                        nout1=ia-atn
                        do 2040,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2040,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
                        zout(1,j,nout1)= r2 + r1
                        zout(2,j,nout1)= s2 + s1
                        zout(1,j,nout2)= r1 - r2
                        zout(2,j,nout2)= s1 - s2
2040                        continue
                endif
        else if (4*ias.eq.3*after) then
                if (isign.eq.1) then
                        nin1=ia-after
                        nout1=ia-atn
                        do 2050,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2050,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
                        zout(1,j,nout1)= r1 - r2
                        zout(2,j,nout1)= s2 + s1
                        zout(1,j,nout2)= r2 + r1
                        zout(2,j,nout2)= s1 - s2
2050                        continue
                else
                        nin1=ia-after
                        nout1=ia-atn
                        do 2060,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2060,j=1,nfft
                        r1=zin(1,j,nin1)
                        s1=zin(2,j,nin1)
                        r=zin(1,j,nin2)
                        s=zin(2,j,nin2)
                        r2=(s - r)*rt2i
                        s2=(r + s)*rt2i
                        zout(1,j,nout1)= r2 + r1
                        zout(2,j,nout1)= s1 - s2
                        zout(1,j,nout2)= r1 - r2
                        zout(2,j,nout2)= s2 + s1
2060                        continue
                endif
        else
                itrig=ias*before+1
                cr2=trig(1,itrig)
                ci2=trig(2,itrig)
                nin1=ia-after
                nout1=ia-atn
                do 2090,ib=1,before
                nin1=nin1+after
                nin2=nin1+atb
                nout1=nout1+atn
                nout2=nout1+after
                do 2090,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
                zout(1,j,nout1)= r2 + r1
                zout(2,j,nout1)= s2 + s1
                zout(1,j,nout2)= r1 - r2
                zout(2,j,nout2)= s1 - s2
2090                continue
        endif
2000        continue
        else 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 = ( dm - cp)*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
                do 8000,ia=2,after
                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)
                        itrig=itrig+itt
                        cr6=trig(1,itrig)
                        ci6=trig(2,itrig)
                        itrig=itrig+itt
                        cr7=trig(1,itrig)
                        ci7=trig(2,itrig)
                        itrig=itrig+itt
                        cr8=trig(1,itrig)
                        ci8=trig(2,itrig)
                        nin1=ia-after
                        nout1=ia-atn
                        do 8020,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 8020,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
                        r=zin(1,j,nin6)
                        s=zin(2,j,nin6)
                        r6=r*cr6 - s*ci6
                        s6=r*ci6 + s*cr6
                        r=zin(1,j,nin7)
                        s=zin(2,j,nin7)
                        r7=r*cr7 - s*ci7
                        s7=r*ci7 + s*cr7
                        r=zin(1,j,nin8)
                        s=zin(2,j,nin8)
                        r8=r*cr8 - s*ci8
                        s8=r*ci8 + s*cr8
                        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 = ( dm - cp)*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

8020                        continue
8000                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= ( dp - cm)*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

                do 8001,ia=2,after
                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)
                        itrig=itrig+itt
                        cr6=trig(1,itrig)
                        ci6=trig(2,itrig)
                        itrig=itrig+itt
                        cr7=trig(1,itrig)
                        ci7=trig(2,itrig)
                        itrig=itrig+itt
                        cr8=trig(1,itrig)
                        ci8=trig(2,itrig)
                        nin1=ia-after
                        nout1=ia-atn
                        do 8021,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 8021,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
                        r=zin(1,j,nin6)
                        s=zin(2,j,nin6)
                        r6=r*cr6 - s*ci6
                        s6=r*ci6 + s*cr6
                        r=zin(1,j,nin7)
                        s=zin(2,j,nin7)
                        r7=r*cr7 - s*ci7
                        s7=r*ci7 + s*cr7
                        r=zin(1,j,nin8)
                        s=zin(2,j,nin8)
                        r8=r*cr8 - s*ci8
                        s8=r*ci8 + s*cr8
                        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= ( dp - cm)*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
8021                        continue
8001                continue

        endif
        else if (now.eq.3) then 
!         .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=r3 + r2
                s=s2 - s3
                zout(1,j,nout1) = r1 - r
                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) = s1 - s
                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=(s - r)*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
!         cos(2.d0*pi/5.d0)
        cos2=0.3090169943749474d0
!         cos(4.d0*pi/5.d0)
        cos4=-0.8090169943749474d0
!        sin(2.d0*pi/5.d0)
        sin2=isign*0.9510565162951536d0
!         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 = r1 + cos2*r25 + cos4*r34
        s = sin2*s25 + sin4*s34
        zout(1,j,nout2) = r - s
        zout(1,j,nout5) = r + s
        r = r1 + cos4*r25 + cos2*r34
        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 = s1 + cos2*s25 + cos4*s34
        s = sin2*r25 + sin4*r34
        zout(2,j,nout2) = r + s
        zout(2,j,nout5) = r - s
        r = s1 + cos4*s25 + cos2*s34
        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 = r1 + cos2*r25 - cos4*r34 
                        s = sin2*s25 + sin4*s34
                        zout(1,j,nout2) = r - s
                        zout(1,j,nout5) = r + s
                        r = r1 + cos4*r25 - cos2*r34 
                        s = sin4*s25 - sin2*s34
                        zout(1,j,nout3) = r - s
                        zout(1,j,nout4) = r + s
                        r25 = r2 + r5
                        r34 = r4 - r3
                        s25 = s2 - s5
                        s34 = s3 + s4
                        zout(2,j,nout1) = s1 + s25 + s34
                        r = s1 + cos2*s25 + cos4*s34
                        s = sin2*r25 + sin4*r34
                        zout(2,j,nout2) = r + s
                        zout(2,j,nout5) = r - s
                        r = s1 + cos4*s25 + cos2*s34
                        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=(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
                        r5=zin(1,j,nin5)
                        s5=zin(2,j,nin5)
                        r25 = r2 - r5
                        r34 = r3 + r4
                        s25 = s2 + s5
                        s34 = s4 - s3
                        zout(1,j,nout1) = r1 + r25 + r34
                        r = r1 + cos2*r25 + cos4*r34
                        s = sin2*s25 + sin4*s34
                        zout(1,j,nout2) = r - s
                        zout(1,j,nout5) = r + s
                        r = r1 + cos4*r25 + cos2*r34
                        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 = s1 + cos2*s25 - cos4*s34
                        s = sin2*r25 + sin4*r34
                        zout(2,j,nout2) = r + s
                        zout(2,j,nout5) = r - s
                        r = s1 + cos4*s25 - cos2*s34
                        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 = r1 + cos2*r25 + cos4*r34
                s = sin2*s25 + sin4*s34
                zout(1,j,nout2) = r - s
                zout(1,j,nout5) = r + s
                r = r1 + cos4*r25 + cos2*r34
                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 = s1 + cos2*s25 + cos4*s34
                s = sin2*r25 + sin4*r34
                zout(2,j,nout2) = r + s
                zout(2,j,nout5) = r - s
                r = s1 + cos4*s25 + cos2*s34
                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
!         .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 if (now.eq.7) then

      	  ia=1
      	  nin1=ia-after
      	  nout1=ia-atn

      	  do 7000,ib=1,before
      	 	 nin1=nin1+after
      	 	 nin2=nin1+atb
      	 	 nin3=nin2+atb
      	 	 nin4=nin3+atb
      	 	 nin5=nin4+atb
      	 	 nin6=nin5+atb
      	 	 nin7=nin6+atb
      	 	 nout1=nout1+atn
      	 	 nout2=nout1+after
      	 	 nout3=nout2+after
      	 	 nout4=nout3+after
      	 	 nout5=nout4+after
      	 	 nout6=nout5+after
      	 	 nout7=nout6+after
      	 	 do 7000,j=1,nfft

				cos2= 0.62348980185873353053d0
				cos4=-0.22252093395631440429d0
				cos6=-0.90096886790241912624d0
				
				sin2= 0.78183148246802980871d0*isign
				sin4= 0.97492791218182360702d0*isign
				sin6= 0.43388373911755812048d0*isign
				
				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)
				
				rp2=r2+r7; rp3=r3+r6; rp4=r4+r5
				rm2=r2-r7; rm3=r3-r6; rm4=r4-r5
				sp2=s2+s7; sp3=s3+s6; sp4=s4+s5
				sm2=s2-s7; sm3=s3-s6; sm4=s4-s5 !3*4=12 flops
				
				ur2=r1+cos2*rp2+cos4*rp3+cos6*rp4
				vr2=   sin2*sm2+sin4*sm3+sin6*sm4
				
				ui2=s1+cos2*sp2+cos4*sp3+cos6*sp4
				vi2=   sin2*rm2+sin4*rm3+sin6*rm4
				          
				ur3=r1+cos4*rp2+cos6*rp3+cos2*rp4
				vr3=   sin4*sm2-sin6*sm3-sin2*sm4
				
				ui3=s1+cos4*sp2+cos6*sp3+cos2*sp4
				vi3=   sin4*rm2-sin6*rm3-sin2*rm4
				
				ur4=r1+cos6*rp2+cos2*rp3+cos4*rp4
				vr4=   sin6*sm2-sin2*sm3+sin4*sm4
				
				ui4=s1+cos6*sp2+cos2*sp3+cos4*sp4
				vi4=   sin6*rm2-sin2*rm3+sin4*rm4   ! 11*6=66 flops, of them
				!6*6=36 multiplications
				
				zout(1,j,nout1)=r1+rp2+rp3+rp4; zout(2,j,nout1)=s1+sp2+sp3+sp4
				zout(1,j,nout2)=ur2-vr2;        zout(2,j,nout2)=ui2+vi2
				zout(1,j,nout3)=ur3-vr3;        zout(2,j,nout3)=ui3+vi3
				zout(1,j,nout4)=ur4-vr4;        zout(2,j,nout4)=ui4+vi4
				zout(1,j,nout5)=ur4+vr4;        zout(2,j,nout5)=ui4-vi4
				zout(1,j,nout6)=ur3+vr3;        zout(2,j,nout6)=ui3-vi3
				zout(1,j,nout7)=ur2+vr2;        zout(2,j,nout7)=ui2-vi2 
				!(3+6)*2 =18 flops

				! In total: 12+66+18=96flops
		
7000  	  continue
        else 
        stop 'error fftstp'
        endif

        return
        end



         subroutine fftrot(mm,nfft,m,nn,n,zin,zout,ntrig,trig,after,now,before,isign)
!  Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
!  Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1995, 1999
!  This file is distributed under the terms of the
!  GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .

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

!         sqrt(.5d0)
        rt2i=0.7071067811865475d0
        if (now.eq.2) then
        ia=1
        nin1=ia-after
        nout1=ia-atn
        do 2001,ib=1,before
        nin1=nin1+after
        nin2=nin1+atb
        nout1=nout1+atn
        nout2=nout1+after
        do 2001,j=1,nfft
        r1=zin(1,j,nin1)
        s1=zin(2,j,nin1)
        r2=zin(1,j,nin2)
        s2=zin(2,j,nin2)
        zout(1,nout1,j)= r2 + r1
        zout(2,nout1,j)= s2 + s1
        zout(1,nout2,j)= r1 - r2
        zout(2,nout2,j)= s1 - s2
2001        continue
        do 2000,ia=2,after
        ias=ia-1
        if (2*ias.eq.after) then
                if (isign.eq.1) then
                        nin1=ia-after
                        nout1=ia-atn
                        do 2010,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2010,j=1,nfft
                        r1=zin(1,j,nin1)
                        s1=zin(2,j,nin1)
                        r2=zin(2,j,nin2)
                        s2=zin(1,j,nin2)
                        zout(1,nout1,j)= r1 - r2
                        zout(2,nout1,j)= s2 + s1
                        zout(1,nout2,j)= r2 + r1
                        zout(2,nout2,j)= s1 - s2
2010                        continue
                else
                        nin1=ia-after
                        nout1=ia-atn
                        do 2020,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2020,j=1,nfft
                        r1=zin(1,j,nin1)
                        s1=zin(2,j,nin1)
                        r2=zin(2,j,nin2)
                        s2=zin(1,j,nin2)
                        zout(1,nout1,j)= r2 + r1
                        zout(2,nout1,j)= s1 - s2
                        zout(1,nout2,j)= r1 - r2
                        zout(2,nout2,j)= s2 + s1
2020                        continue
                endif
        else if (4*ias.eq.after) then
                if (isign.eq.1) then
                        nin1=ia-after
                        nout1=ia-atn
                        do 2030,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2030,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
                        zout(1,nout1,j)= r2 + r1
                        zout(2,nout1,j)= s2 + s1
                        zout(1,nout2,j)= r1 - r2
                        zout(2,nout2,j)= s1 - s2
2030                        continue
                else
                        nin1=ia-after
                        nout1=ia-atn
                        do 2040,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2040,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
                        zout(1,nout1,j)= r2 + r1
                        zout(2,nout1,j)= s2 + s1
                        zout(1,nout2,j)= r1 - r2
                        zout(2,nout2,j)= s1 - s2
2040                        continue
                endif
        else if (4*ias.eq.3*after) then
                if (isign.eq.1) then
                        nin1=ia-after
                        nout1=ia-atn
                        do 2050,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2050,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
                        zout(1,nout1,j)= r1 - r2
                        zout(2,nout1,j)= s2 + s1
                        zout(1,nout2,j)= r2 + r1
                        zout(2,nout2,j)= s1 - s2
2050                        continue
                else
                        nin1=ia-after
                        nout1=ia-atn
                        do 2060,ib=1,before
                        nin1=nin1+after
                        nin2=nin1+atb
                        nout1=nout1+atn
                        nout2=nout1+after
                        do 2060,j=1,nfft
                        r1=zin(1,j,nin1)
                        s1=zin(2,j,nin1)
                        r=zin(1,j,nin2)
                        s=zin(2,j,nin2)
                        r2=(s - r)*rt2i
                        s2=(r + s)*rt2i
                        zout(1,nout1,j)= r2 + r1
                        zout(2,nout1,j)= s1 - s2
                        zout(1,nout2,j)= r1 - r2
                        zout(2,nout2,j)= s2 + s1
2060                        continue
                endif
        else
                itrig=ias*before+1
                cr2=trig(1,itrig)
                ci2=trig(2,itrig)
                nin1=ia-after
                nout1=ia-atn
                do 2090,ib=1,before
                nin1=nin1+after
                nin2=nin1+atb
                nout1=nout1+atn
                nout2=nout1+after
                do 2090,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
                zout(1,nout1,j)= r2 + r1
                zout(2,nout1,j)= s2 + s1
                zout(1,nout2,j)= r1 - r2
                zout(2,nout2,j)= s1 - s2
2090                continue
        endif
2000        continue
        else 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 = ( dm - cp)*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
                do 8000,ia=2,after
                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)
                        itrig=itrig+itt
                        cr6=trig(1,itrig)
                        ci6=trig(2,itrig)
                        itrig=itrig+itt
                        cr7=trig(1,itrig)
                        ci7=trig(2,itrig)
                        itrig=itrig+itt
                        cr8=trig(1,itrig)
                        ci8=trig(2,itrig)
                        nin1=ia-after
                        nout1=ia-atn
                        do 8020,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 8020,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
                        r=zin(1,j,nin6)
                        s=zin(2,j,nin6)
                        r6=r*cr6 - s*ci6
                        s6=r*ci6 + s*cr6
                        r=zin(1,j,nin7)
                        s=zin(2,j,nin7)
                        r7=r*cr7 - s*ci7
                        s7=r*ci7 + s*cr7
                        r=zin(1,j,nin8)
                        s=zin(2,j,nin8)
                        r8=r*cr8 - s*ci8
                        s8=r*ci8 + s*cr8
                        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 = ( dm - cp)*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

8020                        continue
8000                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= ( dp - cm)*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

                do 8001,ia=2,after
                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)
                        itrig=itrig+itt
                        cr6=trig(1,itrig)
                        ci6=trig(2,itrig)
                        itrig=itrig+itt
                        cr7=trig(1,itrig)
                        ci7=trig(2,itrig)
                        itrig=itrig+itt
                        cr8=trig(1,itrig)
                        ci8=trig(2,itrig)
                        nin1=ia-after
                        nout1=ia-atn
                        do 8021,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 8021,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
                        r=zin(1,j,nin6)
                        s=zin(2,j,nin6)
                        r6=r*cr6 - s*ci6
                        s6=r*ci6 + s*cr6
                        r=zin(1,j,nin7)
                        s=zin(2,j,nin7)
                        r7=r*cr7 - s*ci7
                        s7=r*ci7 + s*cr7
                        r=zin(1,j,nin8)
                        s=zin(2,j,nin8)
                        r8=r*cr8 - s*ci8
                        s8=r*ci8 + s*cr8
                        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= ( dp - cm)*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
8021                        continue
8001                continue

        endif
        else if (now.eq.3) then 
!         .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) = r1 - r
                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) = s1 - s
                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=(s - r)*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
!         cos(2.d0*pi/5.d0)
        cos2=0.3090169943749474d0
!         cos(4.d0*pi/5.d0)
        cos4=-0.8090169943749474d0
!        sin(2.d0*pi/5.d0)
        sin2=isign*0.9510565162951536d0
!         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 = r1 + cos2*r25 + cos4*r34
        s = sin2*s25 + sin4*s34
        zout(1,nout2,j) = r - s
        zout(1,nout5,j) = r + s
        r = r1 + cos4*r25 + cos2*r34
        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 = s1 + cos2*s25 + cos4*s34
        s = sin2*r25 + sin4*r34
        zout(2,nout2,j) = r + s
        zout(2,nout5,j) = r - s
        r = s1 + cos4*s25 + cos2*s34
        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 = r1 + cos2*r25 - cos4*r34
                        s = sin2*s25 + sin4*s34
                        zout(1,nout2,j) = r - s
                        zout(1,nout5,j) = r + s
                        r = r1 + cos4*r25 - cos2*r34
                        s = sin4*s25 - sin2*s34
                        zout(1,nout3,j) = r - s
                        zout(1,nout4,j) = r + s
                        r25 = r2 + r5
                        r34 = r4 - r3
                        s25 = s2 - s5
                        s34 = s3 + s4
                        zout(2,nout1,j) = s1 + s25 + s34
                        r = s1 + cos2*s25 + cos4*s34
                        s = sin2*r25 + sin4*r34
                        zout(2,nout2,j) = r + s
                        zout(2,nout5,j) = r - s
                        r = s1 + cos4*s25 + cos2*s34
                        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=(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
                        r5=zin(1,j,nin5)
                        s5=zin(2,j,nin5)
                        r25 = r2 - r5
                        r34 = r3 + r4
                        s25 = s2 + s5
                        s34 = s4 - s3
                        zout(1,nout1,j) = r1 + r25 + r34
                        r = r1 + cos2*r25 + cos4*r34
                        s = sin2*s25 + sin4*s34
                        zout(1,nout2,j) = r - s
                        zout(1,nout5,j) = r + s
                        r = r1 + cos4*r25 + cos2*r34
                        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 = s1 + cos2*s25 - cos4*s34
                        s = sin2*r25 + sin4*r34
                        zout(2,nout2,j) = r + s
                        zout(2,nout5,j) = r - s
                        r = s1 + cos4*s25 - cos2*s34
                        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 = r1 + cos2*r25 + cos4*r34
                s = sin2*s25 + sin4*s34
                zout(1,nout2,j) = r - s
                zout(1,nout5,j) = r + s
                r = r1 + cos4*r25 + cos2*r34
                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 = s1 + cos2*s25 + cos4*s34
                s = sin2*r25 + sin4*r34
                zout(2,nout2,j) = r + s
                zout(2,nout5,j) = r - s
                r = s1 + cos4*s25 + cos2*s34
                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
!         .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 if (now.eq.7) then

      	  ia=1
      	  nin1=ia-after
      	  nout1=ia-atn

      	  do 7000,ib=1,before
      	 	 nin1=nin1+after
      	 	 nin2=nin1+atb
      	 	 nin3=nin2+atb
      	 	 nin4=nin3+atb
      	 	 nin5=nin4+atb
      	 	 nin6=nin5+atb
      	 	 nin7=nin6+atb
      	 	 nout1=nout1+atn
      	 	 nout2=nout1+after
      	 	 nout3=nout2+after
      	 	 nout4=nout3+after
      	 	 nout5=nout4+after
      	 	 nout6=nout5+after
      	 	 nout7=nout6+after
      	 	 do 7000,j=1,nfft

				cos2= 0.62348980185873353053d0
				cos4=-0.22252093395631440429d0
				cos6=-0.90096886790241912624d0
				
				sin2= 0.78183148246802980871d0*isign
				sin4= 0.97492791218182360702d0*isign
				sin6= 0.43388373911755812048d0*isign
				
				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)
				
				rp2=r2+r7; rp3=r3+r6; rp4=r4+r5
				rm2=r2-r7; rm3=r3-r6; rm4=r4-r5
				sp2=s2+s7; sp3=s3+s6; sp4=s4+s5
				sm2=s2-s7; sm3=s3-s6; sm4=s4-s5
				
				ur2=r1+cos2*rp2+cos4*rp3+cos6*rp4
				vr2=   sin2*sm2+sin4*sm3+sin6*sm4
				
				ui2=s1+cos2*sp2+cos4*sp3+cos6*sp4
				vi2=   sin2*rm2+sin4*rm3+sin6*rm4
				          
				ur3=r1+cos4*rp2+cos6*rp3+cos2*rp4
				vr3=   sin4*sm2-sin6*sm3-sin2*sm4
				
				ui3=s1+cos4*sp2+cos6*sp3+cos2*sp4
				vi3=   sin4*rm2-sin6*rm3-sin2*rm4
				
				ur4=r1+cos6*rp2+cos2*rp3+cos4*rp4
				vr4=   sin6*sm2-sin2*sm3+sin4*sm4
				
				ui4=s1+cos6*sp2+cos2*sp3+cos4*sp4
				vi4=   sin6*rm2-sin2*rm3+sin4*rm4
				
				zout(1,nout1,j)=r1+rp2+rp3+rp4; zout(2,nout1,j)=s1+sp2+sp3+sp4
				zout(1,nout2,j)=ur2-vr2;        zout(2,nout2,j)=ui2+vi2
				zout(1,nout3,j)=ur3-vr3;        zout(2,nout3,j)=ui3+vi3
				zout(1,nout4,j)=ur4-vr4;        zout(2,nout4,j)=ui4+vi4
				zout(1,nout5,j)=ur4+vr4;        zout(2,nout5,j)=ui4-vi4
				zout(1,nout6,j)=ur3+vr3;        zout(2,nout6,j)=ui3-vi3
				zout(1,nout7,j)=ur2+vr2;        zout(2,nout7,j)=ui2-vi2
		
7000  	  continue
       else
        stop 'error fftrot'
       endif

        return
        end

