cffti1 Subroutine

subroutine cffti1(n, wa, ifac)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n
real(kind=dp), intent(out) :: wa(*)
integer, intent(out) :: ifac(*)

Contents

Source Code


Variables

Type Visibility Attributes Name Initial
real(kind=dp), public :: arg
real(kind=dp), public :: argh
real(kind=dp), public :: argld
real(kind=dp), public :: fi
integer, public :: i
integer, public :: i1
integer, public :: ib
integer, public :: ido
integer, public :: idot
integer, public :: ii
integer, public :: ip
integer, public :: ipm
integer, public :: j
integer, public :: k1
integer, public :: l1
integer, public :: l2
integer, public :: ld
integer, public :: nf
integer, public :: nl
integer, public :: nq
integer, public :: nr
integer, public :: ntry
integer, public, parameter, dimension(4) :: ntryh = [3, 4, 2, 5]
real(kind=dp), public, parameter :: tpi = 2.0_dp*acos(-1.0_dp)

Source Code

      subroutine cffti1(n, wa, ifac)
         use fftpack_kind, only: dp => rk
         implicit none
         integer, intent(in) :: n
         integer, intent(out) :: ifac(*)
         real(dp), intent(out) :: wa(*)
         real(dp) :: arg, argh, argld, fi
         integer :: i, i1, ib, ido, idot, ii, ip, ipm, j, k1, &
                    l1, l2, ld, nf, nl, nq, nr, ntry
         integer, dimension(4), parameter :: ntryh = [3, 4, 2, 5]
         real(dp), parameter :: tpi = 2.0_dp*acos(-1.0_dp) ! 2 * pi
         nl = n
         nf = 0
         j = 0
100      j = j + 1
         if (j <= 4) then
            ntry = ntryh(j)
         else
            ntry = ntry + 2
         end if
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)
               end do
               ifac(3) = 2
            end if
         end if
         if (nl /= 1) goto 200
         ifac(1) = n
         ifac(2) = nf
         argh = tpi/real(n, kind=dp)
         i = 2
         l1 = 1
         do k1 = 1, nf
            ip = ifac(k1 + 2)
            ld = 0
            l2 = l1*ip
            ido = n/l2
            idot = ido + ido + 2
            ipm = ip - 1
            do j = 1, ipm
               i1 = i
               wa(i - 1) = 1.0_dp
               wa(i) = 0.0_dp
               ld = ld + l1
               fi = 0.0_dp
               argld = real(ld, kind=dp)*argh
               do ii = 4, idot, 2
                  i = i + 2
                  fi = fi + 1.0_dp
                  arg = fi*argld
                  wa(i - 1) = cos(arg)
                  wa(i) = sin(arg)
               end do
               if (ip > 5) then
                  wa(i1 - 1) = wa(i - 1)
                  wa(i1) = wa(i)
               end if
            end do
            l1 = l2
         end do
      end subroutine cffti1