***********************************************************************
*     cft512 - 2-dimensional FFT (512x512)                            *
*           coordinate of the DC component of the output : (257,257)  *
*       isw   : input,integer direction of FFT                        *
*       isw   = 1                                                     *
*         a : input,real  real part of the original map               *
*             output,real cos of Fourier map                          *
*         b : input,real  imaginary part of the original map          *
*             output,real sin of Fourier map                          *
*       isw = -1                                                      *
*         a : input,real(512,512)  cos of Fourier map                 *
*             output,real(512,512) real part of the output map        *
*         b : input,real(512,512)  sin of Fourier map                 *
*             output,real(512,512) imaginary part of the output map   *
*                                         April 8, 1996 K. Fujiki     *
***********************************************************************
      subroutine cft512(a,b,isw)
      
      integer isw,iw(2052),iflag
      real a(512,512),b(512,512),w1(9216),w2(524288)

      if(isw.eq.0) then
        call vfctsi(512,512,iw,w1,ierr)
        if(ierr.ne.0) write(6,*)'initialize error'
        return
      end if


      is=256

      a=cshift(a,is,1)
      a=cshift(a,is,2)
      b=cshift(b,is,1)
      b=cshift(b,is,2)

      call vfctsf(512,512,a,b,512,isw,iw,w1,w2,ierr)
      if(ierr.ne.0) then
        write(6,*)'FFT error. ierr= ',ierr
      end if

      a=cshift(a,-is,1)
      a=cshift(a,-is,2)
      b=cshift(b,-is,1)
      b=cshift(b,-is,2)

      if(isw.eq.1) then
        a=a/512./512.
        b=b/512./512.
      end if

      return
      end


***********************************************************************
*     cft1024 - 2-dimensional FFT (1024x1024)                         *
*           coordinate of the DC component of the output : (513,513)  *
*       isw   : input,integer direction of FFT                        *
*       isw   = 1                                                     *
*         a : input,real  real part of the original map               *
*             output,real cos of Fourier map                          *
*         b : input,real  imaginary part of the original map          *
*             output,real sin of Fourier map                          *
*       isw = -1                                                      *
*         a : input,real(512,512)  cos of Fourier map                 *
*             output,real(512,512) real part of the output map        *
*         b : input,real(512,512)  sin of Fourier map                 *
*             output,real(512,512) imaginary part of the output map   *
*                                         April 8, 1996 K. Fujiki     *
***********************************************************************
      subroutine cft1024(a,b,isw)

****   1024x1024 array
      integer isw,iw(4100)
      real w1(20480),w2(2097152)
      real a(1024,1024),b(1024,1024)

      if(isw.eq.0) then
        write(6,*) 'Initialize routine for 1024x1024 matrix....'
        call vfctsi(1024,1024,iw,w1,ierr)
        if(ierr.ne.0) then
          write(6,*)'cft1024 initialize error ',ierr
        else
          write(6,*)'cft1024 initialize success ',ierr
        end if
c        write(8,*)'Done...'
        return
      end if


      is=512

      a=cshift(a,is,1)
      a=cshift(a,is,2)
      b=cshift(b,is,1)
      b=cshift(b,is,2)

      call vfctsf(1024,1024,a,b,1024,isw,iw,w1,w2,ierr)

      if(ierr.ne.0) write(6,*)'FFT error. ierr= ',ierr

      a=cshift(a,-is,1)
      a=cshift(a,-is,2)
      b=cshift(b,-is,1)
      b=cshift(b,-is,2)

      if(isw.eq.1) then
        a=a/1024./1024.
        b=b/1024./1024.
      end if

      return
      end
c
***********************************************************************
*     cft128 - 2-dimensional FFT (128x128)                            *
*           coordinate of the DC component of the output : (65,65)    *
*       istat : input,integer direction of FFT                        *
*       istat = 1                                                     *
*         a : input,real  real part of the original map               *
*             output,real cos of Fourier map                          *
*         b : input,real  imaginary part of the original map          *
*             output,real sin of Fourier map                          *
*       istat = -1                                                    *
*         a : input,real(128,128)  cos of Fourier map                 *
*             output,real(128,128) real part of the output map        *
*         b : input,real(128,128)  sin of Fourier map                 *
*             output,real(128,128) imaginary part of the output map   *
*                                         April 8, 1996 K. Fujiki     *
***********************************************************************
      subroutine cft128(a,b,isw)
      
      integer isw,iw(516)
      real a(128,128),b(128,128),w1(1792),w2(32768)

      if(isw.eq.0) then
        write(6,*) 'Initialize....'
        call vfctsi(128,128,iw,w1,ierr)
        if(ierr.ne.0) then
          write(6,*)'initialize error ',ierr
        else
          write(6,*)'initialize success ',ierr
        end if
        return
      end if


      is=64

      a=cshift(a,is,1)
      a=cshift(a,is,2)
      b=cshift(b,is,1)
      b=cshift(b,is,2)

      call vfctsf(128,128,a,b,128,isw,iw,w1,w2,ierr)
      if(ierr.ne.0) then
        write(6,*)'FFT error. ierr= ',ierr
      end if

      a=cshift(a,-is,1)
      a=cshift(a,-is,2)
      b=cshift(b,-is,1)
      b=cshift(b,-is,2)

      if(isw.eq.1) then
        a=a/128./128.
        b=b/128./128.
      end if

      return
      end