* how to compile: f77 -dW -c -h ary90 cft_asl.f
*
***********************************************************************
*     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     *
*               modified for the new ASL  April 2, 1999 Y. Hanaoka    *
***********************************************************************
c
      subroutine cft512(a,b,isw)
      
c      integer isw,iw(2052),iflag
      integer isw,ifax(40)
c      real a(512,512),b(512,512),w1(9216),w2(524288)
      real a(512,512),b(512,512),aa(513,513),bb(513,513),
     -     trigs(2048),wk(526338)
      parameter(is=256)

      if(isw.eq.0) then
c        call vfctsi(512,512,iw,w1,ierr)
        call rfc2fb(512,512,aa,bb,513,513,isw,ifax,trigs,wk,ierr)
        if(ierr.ne.0) write(6,*)'initialize error'
        return
      end if


c      is=256

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

      aa(1:512,1:512)=a(1:512,1:512)
      bb(1:512,1:512)=b(1:512,1:512)

c      call vfctsf(512,512,a,b,512,isw,iw,w1,w2,ierr)
      call rfc2bf(512,512,aa,bb,513,513,isw,ifax,trigs,wk,ierr)
      if(ierr.ne.0) then
        write(6,*)'FFT error. ierr= ',ierr
      end if

      a(1:512,1:512)=aa(1:512,1:512)
      b(1:512,1:512)=bb(1:512,1:512)

      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     *
*               modified for the new ASL  April 2, 1999 Y. Hanaoka    *
***********************************************************************
      subroutine cft1024(a,b,isw)

****   1024x1024 array
c      integer isw,iw(4100)
      integer isw,ifax(40)
c      real w1(20480),w2(2097152)
c      real a(1024,1024),b(1024,1024)
      real a(1024,1024),b(1024,1024),aa(1025,1025),bb(1025,1025),
     -     trigs(4096),wk(2101250)
      parameter(is=512)

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


c      is=512

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

      aa(1:1024,1:1024)=a(1:1024,1:1024)
      bb(1:1024,1:1024)=b(1:1024,1:1024)

c      call vfctsf(1024,1024,a,b,1024,isw,iw,w1,w2,ierr)
      call rfc2bf(1024,1024,aa,bb,1025,1025,isw,ifax,trigs,wk,ierr)

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

      a(1:1024,1:1024)=aa(1:1024,1:1024)
      b(1:1024,1:1024)=bb(1:1024,1:1024)

      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     *
*               modified for the new ASL  April 2, 1999 Y. Hanaoka    *
***********************************************************************
      subroutine cft128(a,b,isw)
      
c      integer isw,iw(516)
      integer isw,ifax(40)
c      real a(128,128),b(128,128),w1(1792),w2(32768)
      real a(128,128),b(128,128),aa(129,129),bb(129,129),
     -     trigs(512),wk(33282)

      if(isw.eq.0) then
        write(6,*) 'Initialize....'
c        call vfctsi(128,128,iw,w1,ierr)
        call rfc2fb(128,128,aa,bb,129,129,isw,ifax,trigs,wk,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)

      aa(1:128,1:128)=a(1:128,1:128)
      bb(1:128,1:128)=b(1:128,1:128)

c      call vfctsf(128,128,a,b,128,isw,iw,w1,w2,ierr)
      call rfc2bf(128,128,aa,bb,129,129,isw,ifax,trigs,wk,ierr)
      if(ierr.ne.0) then
        write(6,*)'FFT error. ierr= ',ierr
      end if

      a(1:128,1:128)=aa(1:128,1:128)
      b(1:128,1:128)=bb(1:128,1:128)

      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