Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | n | |||
real(kind=dp), | intent(out) | :: | wa(*) | |||
integer, | intent(out) | :: | ifac(*) |
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) |
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