ezfft1 Subroutine

subroutine ezfft1(n, Wa, Ifac)

Uses

Arguments

Type IntentOptional Attributes Name
integer :: n
real(kind=rk) :: Wa
integer :: Ifac

Contents

Variables

Source Code

ezfft1

Variables

Type Visibility Attributes Name Initial
real(kind=rk), public :: arg1
real(kind=rk), public :: argh
real(kind=rk), public :: ch1
real(kind=rk), public :: ch1h
real(kind=rk), public :: dch1
real(kind=rk), public :: dsh1
integer, public :: i
integer, public :: ib
integer, public :: ido
integer, public :: ii
integer, public :: ip
integer, public :: ipm
integer, public :: is
integer, public :: j
integer, public :: k1
integer, public :: l1
integer, public :: l2
integer, public :: nf
integer, public :: nfm1
integer, public :: nl
integer, public :: nq
integer, public :: nr
integer, public :: ntry
integer, public, parameter, dimension(4) :: ntryh = [4, 2, 3, 5]
real(kind=rk), public :: sh1
real(kind=rk), public, parameter :: tpi = 2.0_rk*acos(-1.0_rk)

Source Code

      subroutine ezfft1(n,Wa,Ifac)
      use fftpack_kind
      implicit none
      real(rk) :: arg1 , argh , ch1 , ch1h , dch1 , dsh1 , sh1 , &
                  Wa
      integer :: i , ib , ido , Ifac , ii , ip , ipm , is , j , k1 , l1 , &
                 l2 , n , nf , nfm1 , nl , nq , nr , ntry
      dimension Wa(*) , Ifac(*)  
      integer,dimension(4),parameter :: ntryh = [4 , 2 , 3 , 5]
      real(rk),parameter :: tpi = 2.0_rk * acos(-1.0_rk) ! 2 * pi
      nl = n
      nf = 0
      j = 0
 100  j = j + 1
      if ( j<=4 ) then
         ntry = ntryh(j)
      else
         ntry = ntry + 2
      endif
 200  nq = nl/ntry
      nr = nl - ntry*nq
      if ( nr/=0 ) goto 100
      nf = nf + 1
      Ifac(nf+2) = ntry
      nl = nq
      if ( ntry==2 ) then
         if ( nf/=1 ) then
            do i = 2 , nf
               ib = nf - i + 2
               Ifac(ib+2) = Ifac(ib+1)
            enddo
            Ifac(3) = 2
         endif
      endif
      if ( nl/=1 ) goto 200
      Ifac(1) = n
      Ifac(2) = nf
      argh = tpi/real(n, rk)
      is = 0
      nfm1 = nf - 1
      l1 = 1
      if ( nfm1==0 ) return
      do k1 = 1 , nfm1
         ip = Ifac(k1+2)
         l2 = l1*ip
         ido = n/l2
         ipm = ip - 1
         arg1 = real(l1, rk)*argh
         ch1 = 1.0_rk
         sh1 = 0.0_rk
         dch1 = cos(arg1)
         dsh1 = sin(arg1)
         do j = 1 , ipm
            ch1h = dch1*ch1 - dsh1*sh1
            sh1 = dch1*sh1 + dsh1*ch1
            ch1 = ch1h
            i = is + 2
            Wa(i-1) = ch1
            Wa(i) = sh1
            if ( ido>=5 ) then
               do ii = 5 , ido , 2
                  i = i + 2
                  Wa(i-1) = ch1*Wa(i-3) - sh1*Wa(i-2)
                  Wa(i) = ch1*Wa(i-2) + sh1*Wa(i-3)
               enddo
            endif
            is = is + ido
         enddo
         l1 = l2
      enddo
      end subroutine ezfft1