fft2.f90

Go to the documentation of this file.

subroutine fft2 (datar, datai, n, inc)
!---------------------------------------------------------------------------
!
!     ******************************************************************
!     * fft2 fortran version clair nielson may 75.                     *
!     ******************************************************************
!
!---------------------------------------------------------------------------

  use itm_types

  implicit none

  integer(itm_i4), intent(in) :: n, inc 
  real(r8), dimension(iabs(inc) * n), intent(inout) :: datar, datai

  real(r8) :: tempr, tempi, theta, sinth, wstpr, wstpi, wr, wi
  integer(itm_i4) :: ktran, ks, ip0, ip1, ip2, ip3, irev, ibit, i
  integer(itm_i4) :: i1, i3, j0, j1

  ktran = sign(-1, inc)
  ks = abs(inc)
  ip0 = ks
  ip3 = ip0 * n
  irev = 1

  do i = 1, ip3, ip0
    if (i < irev) then
      tempr = datar(i)
      tempi = datai(i)
      datar(i) = datar(irev)
      datai(i) = datai(irev)
      datar(irev) = tempr
      datai(irev) = tempi
    end if
    ibit = ip3 / 2
    do
      if (irev <= ibit) exit
      irev = irev - ibit
      ibit = ibit / 2
      if (ibit < ip0) exit    
    end do
    irev = irev + ibit
  end do

  ip1 = ip0
  theta = dble(ktran) * 3.1415926535898_r8
  do
    if (ip1 >= ip3) exit    
    ip2 = ip1 + ip1
    sinth = sin(.5_r8 * theta)
    wstpr = -2._r8 * sinth * sinth
    wstpi = sin(theta)
    wr = 1._r8
    wi = 0._r8
    do i1 = 1, ip1, ip0
      do i3 = i1, ip3, ip2
        j0 = i3
        j1 = j0 + ip1
        tempr = wr * datar(j1) - wi * datai(j1)
        tempi = wr * datai(j1) + wi * datar(j1)
        datar(j1) = datar(j0) - tempr
        datai(j1) = datai(j0) - tempi
        datar(j0) = datar(j0) + tempr
        datai(j0) = datai(j0) + tempi
      end do
      tempr = wr
      wr = wr * wstpr - wi * wstpi + wr
      wi= wi * wstpr + tempr * wstpi + wi
    end do
    ip1 = ip2
    theta = .5_r8 * theta
  end do

  return
end subroutine fft2