subroutine c2r_ffttp_1d( nx, a, b, prim, prim_fact, omega_fix, omegan_fix )
! 1 次元実 FFT 逆変換計算ルーチン
!
! 入力配列は複素数数で行い,実数配列を返す.
  use Math_Const
  implicit none
  integer, intent(in) :: nx  ! 1 次元データ要素数
  complex, intent(in), dimension(0:nx/2-1) :: a  ! 入力複素数データ配列
  real, intent(inout), dimension(0:nx-1) :: b  ! 出力実数データ配列
  character(1), optional, intent(in) :: prim  ! 素因数分解をするかどうか
  ! [o=分解する, x=分解しない] default=分解しない. その場合は, 通常の DFT.
  ! 素因数分解する場合, nx=2^a*3^b*5^c*7^d までしか分解しないようにする.
  ! なお, 実 FFT 変換は現在, データ要素数が偶数の変換しか行わないように
  ! 実装されていることに注意 (奇数データで変換しようとするとエラーを返す).
  ! また, in 属性の引数は, データ数が半分になっていることに注意.
  ! 仕様として, 入力引数の n=0 虚部に nx/2 番目の実数データが入っているように
  ! データを渡す.
  ! ただし, 'x' のときでも, prim_fact が設定されていれば, そのべきで FFT する.
  integer, intent(in), optional, dimension(5) :: prim_fact  ! prim = x のとき設定すると
  ! そのべきで FFT を行う. prim_fact=(/a,b,c,d/) : 2^a*3^b*5^c*7^d
  complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional :: omega_fix  ! 余素因数に対応する回転行列
  complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional :: omegan_fix  ! 回転行列
  complex, dimension(0:nx/2-1) :: c, d
  integer :: i
  integer, dimension(5) :: prim_num

  if(mod(nx,2)/=0)then
     write(*,*) "*** ERROR ***"
     write(*,*) "nx must be even number. Stop."
     stop
  end if

!-- nx/2 個の複素数データ（独立なフーリエ係数）を同数の複素数データに置き換える.
  c(0)=a(0)
  do i=1,nx/2-1
     c(i)=(a(i)+conjg(a(nx/2-i))) +(img*cos(2.0*pi*i/real(nx))-sin(2.0*pi*i/real(nx))) *(a(i)-conjg(a(nx/2-i)))
  end do

!-- FFT 開始

  if(present(prim))then
     if(present(prim_fact))then
        prim_num=prim_fact
        prim_num(1)=prim_num(1)-1    ! データの要素数が半分になっているため,
                                     ! 2 のべき数は 1 つ減っている.
        if(present(omegan_fix))then
           call ffttp_1d( nx/2, c, d, 'i', prim, prim_num, omega_fix(0:nx/2-1,0:nx/2-1), omegan_fix(0:nx/2-1,0:nx/2-1) )
        else
           call ffttp_1d( nx/2, c, d, 'i', prim, prim_num )
        end if
     else
        call ffttp_1d( nx/2, c, d, 'i', prim )
     end if
  else
     call ffttp_1d( nx/2, c, d, 'i' )
  end if

!-- 変換後の配列を整理
  do i=0,nx/2-1
     b(2*i)=real(d(i))
     b(2*i+1)=aimag(d(i))
  end do

end subroutine

subroutine ffttp_1d( nx, a, b, csign, prim, prim_fact, omega_fix, omegan_fix )
!  Temperton's FFT
!
!  1d の fft を csign をもとに, 正変換, 逆変換するルーチン
  use Math_Const
  implicit none
  integer, intent(in) :: nx  ! 入力配列の要素数
  complex, intent(in), dimension(0:nx-1) :: a  ! 入力配列
  complex, intent(inout), dimension(0:nx-1) :: b  ! 出力配列
  character(1), intent(in) :: csign  ! 正逆変換判定 [r=正変換, i=逆変換]
  character(1), optional, intent(in) :: prim  ! 素因数分解をするかどうか
  ! [o=分解する, x=分解しない] default=分解しない. その場合は, 通常の DFT.
  ! 素因数分解する場合, nx=2^a*3^b*5^c*7^d までしか分解しないようにする.
  ! ただし, 'x' のときでも, prim_fact が設定されていれば, そのべきで FFT する.
  integer, dimension(5), optional :: prim_fact  ! prim = x のとき設定すると
  ! そのべきで FFT を行う. prim_fact=(/a,b,c,d,e/) : 2^a*3^b*5^c*7^d*e
  complex, dimension(0:nx-1,0:nx-1), intent(in), optional :: omega_fix
                       ! 余素因数に対応する回転行列
  complex, dimension(0:nx-1,0:nx-1), intent(in), optional :: omegan_fix  ! 回転行列
  integer, allocatable, dimension(:) :: l, m, n  ! 要素等の作業用配列
  complex :: ctmp
  integer :: stat, counter, base
  integer :: i, j, k, id, jd, kd  ! do loop 用配列
  integer, parameter, dimension(4) :: prim_dim=(/2, 3, 5, 7/)  ! 素因数
  integer, dimension(4) :: prim_num  ! 各素因数のべき数
  complex, allocatable, dimension(:,:) :: omega
  complex, dimension(0:nx-1,0:nx-1) :: omegan
  complex, dimension(0:nx-1) :: c  ! tmp array
  complex, dimension(0:1,0:1) :: omega2  ! 波数 2 の回転行列
  complex, dimension(0:2,0:2) :: omega3  ! 波数 3 の回転行列
  complex, dimension(0:4,0:4) :: omega5  ! 波数 5 の回転行列
  complex, dimension(0:6,0:6) :: omega7  ! 波数 7 の回転行列


  base=nx
  prim_num=0
  counter=0

  do i=0,nx-1
     b(i)=a(i)
  end do

  if(real(romega2(1,1))/=-1.0.or.aimag(romega2(1,1))/=1.0.or. real(iomega2(1,1))/=-1.0.or.aimag(iomega2(1,1))/=1.0)then
     ! math_const の rotate_array ルーチンが計算されていない場合, これを作動させる.
     call rotate_array()
  end if

!-- 素因数分解する処理
  if(present(prim))then
     if(prim=='o')then
        if(present(prim_fact))then
           do i=1,4
!              if(prim_fact(i)==0)then
!                 prim_num(i)=1
!              else
                 prim_num(i)=prim_fact(i)
!              end if
              counter=counter+prim_fact(i)
           end do
           base=prim_fact(5)
        else
           call prim_calc( nx, prim_num, base )
           do i=1,4
              counter=counter+prim_num(i)
           end do
        end if
 
        if(base==1)then
           counter=counter-1
        end if
 
        if(counter/=0)then  ! prim=='o' であっても, 素因数分解できなければ DFT に送る.
           allocate(l(counter+1))
           allocate(m(counter+1))
           allocate(n(counter+1))
           stat=0
           l=0
           m=0
           n=0
   
           do i=1,4
              if(prim_num(i)/=0)then
                 select case(prim_dim(i))
                 case(2)
                    n(stat+1:stat+prim_num(i))=2
                 case(3)
                    n(stat+1:stat+prim_num(i))=3
                 case(5)
                    n(stat+1:stat+prim_num(i))=5
                 case(7)
                    n(stat+1:stat+prim_num(i))=7
                 end select
                 stat=stat+prim_num(i)
              end if
           end do
           if(base/=1)then
              n(counter+1)=base
           end if
        end if
     end if
  end if

!-- 回転行列を定義

  if(counter/=0)then
     select case (csign)
     case ('r')
        do j=0,1
           do i=0,1
              omega2(i,j)=romega2(i,j)
           end do
        end do
        do j=0,2
           do i=0,2
              omega3(i,j)=romega3(i,j)
           end do
        end do
        do j=0,4
           do i=0,4
              omega5(i,j)=romega5(i,j)
           end do
        end do
        do j=0,6
           do i=0,6
              omega7(i,j)=romega7(i,j)
           end do
        end do
     case ('i')
        do j=0,1
           do i=0,1
              omega2(i,j)=iomega2(i,j)
           end do
        end do
        do j=0,2
           do i=0,2
              omega3(i,j)=iomega3(i,j)
           end do
        end do
        do j=0,4
           do i=0,4
              omega5(i,j)=iomega5(i,j)
           end do
        end do
        do j=0,6
           do i=0,6
              omega7(i,j)=iomega7(i,j)
           end do
        end do
     end select
  end if

  if(present(omegan_fix))then
     allocate(omega(0:base-1,0:base-1))
     do j=0,base-1
        do i=0,base-1
           omega(i,j)=omega_fix(i,j)
        end do
     end do
     do j=0,nx-1
        do i=0,nx-1
           omegan(i,j)=omegan_fix(i,j)
        end do
     end do
  else
     allocate(omega(0:base-1,0:base-1))
     call rotate_calc( nx, csign, (/prim_num(1), prim_num(2), prim_num(3), prim_num(4), base/), omega, omegan )
  end if

!-- FFT 開始

  if(counter/=0)then
!-- 係数行列定義
     m(1)=1
     l(1)=nx/(n(1)*m(1))
     do i=2,counter+1
        m(i)=m(i-1)*n(i-1)
        l(i)=nx/(n(i)*m(i))
     end do

!-- 変換行列 W の定義

     do kd=1,counter+1
        do jd=0,l(kd)-1
           do id=0,n(kd)-1
              do k=0,m(kd)-1
                 ctmp=b(jd*m(kd)+k)
                 do j=1,n(kd)-1
                    select case(n(kd))
                    case(2)
                       ctmp=ctmp+omega2(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k)

                    case(3)
                       ctmp=ctmp+omega3(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k)

                    case(5)
                       ctmp=ctmp+omega5(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k)

                    case(7)
                       ctmp=ctmp+omega7(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k)

                    case default
                       ctmp=ctmp+omega(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k)

                    end select
                 end do
                 c(jd*n(kd)*m(kd)+id*m(kd)+k)=ctmp*omegan(m(kd),(id*jd))
              end do
           end do
        end do
        do id=0,nx-1
           b(id)=c(id)
        end do
     end do

  else

     do j=0,nx-1
        b(j)=a(0)
        do i=1,nx-1
           b(j)=b(j)+a(i)*omegan(i,j)
        end do
     end do

  end if

  if(csign=='r')then
     do j=0,nx-1
        b(j)=b(j)/real(nx)
     end do
  end if

end subroutine

subroutine ffttp_2d( nx, ny, a, b, csign, prim, prim_fact, omega_fix, omegan_fix )
!  Temperton's FFT (2d ver)
!
!  2d の fft を csign をもとに, 正変換, 逆変換するルーチン
  use Math_Const
  implicit none
  integer, intent(in) :: nx  ! 入力配列の要素数 1
  integer, intent(in) :: ny  ! 入力配列の要素数 2
  complex, intent(in), dimension(0:nx-1,0:ny-1) :: a  ! 入力配列
  complex, intent(inout), dimension(0:nx-1,0:ny-1) :: b  ! 出力配列
  character(1), intent(in) :: csign  ! 正逆変換判定 [r=正変換, i=逆変換]
  character(1), optional, intent(in) :: prim  ! 素因数分解をするかどうか
  ! [o=分解する, x=分解しない] default=分解しない. その場合は, 通常の DFT.
  ! 素因数分解する場合, nx=2^a*3^b*5^c*7^d までしか分解しないようにする.
  ! ただし, 'x' のときでも, prim_fact が設定されていれば, そのべきで FFT する.
  integer, dimension(5), optional :: prim_fact  ! prim = x のとき設定すると
  ! そのべきで FFT を行う. prim_fact=(/a,b,c,d/) : 2^a*3^b*5^c*7^d
  ! ここで, prim_fact で設定するべき数は nx 方向のべき数であることに注意.
  complex, dimension(0:nx-1,0:nx-1), intent(in), optional :: omega_fix  ! 余素因数に対応する回転行列
  complex, dimension(0:nx-1,0:nx-1), intent(in), optional :: omegan_fix  ! 回転行列
  integer :: i

  if(present(prim))then
     if(present(prim_fact))then
        if(present(omegan_fix))then
           do i=0,ny-1
              call ffttp_1d( nx, a(0:nx-1,i), b(0:nx-1,i), csign, prim, prim_fact, omega_fix(0:nx-1,0:nx-1), omegan_fix(0:nx-1,0:nx-1) )
           end do
        else
           do i=0,ny-1
              call ffttp_1d( nx, a(0:nx-1,i), b(0:nx-1,i), csign, prim, prim_fact )
           end do
        end if
     else
        do i=0,ny-1
           call ffttp_1d( nx, a(0:nx-1,i), b(0:nx-1,i), csign, prim )
        end do
     end if
  else
     do i=0,ny-1
        call ffttp_1d( nx, a(0:nx-1,i), b(0:nx-1,i), csign )
     end do
  end if

end subroutine

subroutine prim_calc( n, factor, resid )
  ! 2,3,5,7 についての素因数分解を行う.
  implicit none
  integer, intent(in) :: n  ! 分解する数値
  integer, intent(inout), dimension(4) :: factor  ! 分解したときの各素因数
  ! このとき, factor(1:4)=(/a,b,c,d/) であり, factor(5)=e とすると,
  ! n=2^a*3^b*5^c*7^d*e という式になっている.
  integer, intent(inout) :: resid  ! residual prim factor
  integer :: i, base
  integer, dimension(4) :: prim_dim

  prim_dim=(/2,3,5,7/)
  base=n
  factor=(/0,0,0,0/)
  do i=1,4
     do while(mod(base,prim_dim(i))==0)
        base=base/prim_dim(i)
        factor(i)=factor(i)+1
     end do
  end do

  resid=base

end subroutine

subroutine r2c_ffttp_1d( nx, a, b, prim, prim_fact, omega_fix, omegan_fix )
! 1 次元実 FFT 計算ルーチン
!
! 入力配列は実数で行い, 複素数配列を返す.
  use Math_Const
  implicit none
  integer, intent(in) :: nx  ! 1 次元データ要素数
  real, intent(in), dimension(0:nx-1) :: a  ! 入力実数データ配列
  complex, intent(inout), dimension(0:nx/2-1) :: b  ! 出力複素数データ配列
  character(1), optional, intent(in) :: prim  ! 素因数分解をするかどうか
  ! [o=分解する, x=分解しない] default=分解しない. その場合は, 通常の DFT.
  ! 素因数分解する場合, nx=2^a*3^b*5^c*7^d までしか分解しないようにする.
  ! なお, 実 FFT 変換は現在, データ要素数が偶数の変換しか行わないように
  ! 実装されていることに注意 (奇数データで変換しようとするとエラーを返す).
  ! また, inout 属性の引数は, データ数が半分になっていることに注意.
  ! 仕様として, nx/2 番目の実数データは, n=0 番目の虚部に格納されて
  ! 出力配列に渡されることに注意.
  ! ただし, 'x' のときでも, prim_fact が設定されていれば, そのべきで FFT する.
  integer, dimension(5), intent(in), optional :: prim_fact  ! prim = x のとき設定すると
  ! そのべきで FFT を行う. prim_fact=(/a,b,c,d/) : 2^a*3^b*5^c*7^d
  ! ここで入力するべき数は nx のべき数である.
  complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional :: omega_fix  ! 余素因数に対応する回転行列
  complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional :: omegan_fix  ! 回転行列
  complex, dimension(0:nx/2-1) :: c, d
  integer :: i
  integer, dimension(5) :: prim_num

  if(mod(nx,2)/=0)then
     write(*,*) "*** ERROR ***"
     write(*,*) "nx must be even number. Stop."
     stop
  end if

!-- nx 個の実数データを nx/2 個の複素数データに置き換える.
  do i=0,nx/2-1
     c(i)=a(2*i)+img*a(2*i+1)
  end do

!-- FFT 開始

  if(present(prim))then
     if(present(prim_fact))then
        prim_num=prim_fact
        prim_num(1)=prim_num(1)-1  ! データの要素数が半分になっているため,
                                   ! 2 のべき数は 1 つ減っている.
        if(present(omegan_fix))then
           call ffttp_1d( nx/2, c, d, 'r', prim, prim_num, omega_fix(0:nx/2-1,0:nx/2-1), omegan_fix(0:nx/2-1,0:nx/2-1) )
        else
           call ffttp_1d( nx/2, c, d, 'r', prim, prim_num )
        end if
     else
        call ffttp_1d( nx/2, c, d, 'r', prim )
     end if
  else
     call ffttp_1d( nx/2, c, d, 'r' )
  end if

!-- 変換後の配列を整理
!-- b(N) の実部は b(0) の虚部に組み込むことにする.
!-- b(k) の計算で 0.25 をかけるのは, 上の fft で 2/N で規格化しており,
!-- もとの計算では, 1/2N で規格化しなければならないので,
!-- 1/4 をかけることで, 2/N -> 1/2N で規格化したことになる.
!-- b(0) にかかっている係数もその類.
  b(0)=0.5*((real(d(0))+aimag(d(0)))+img*(real(d(0))-aimag(d(0))))
  do i=1,nx/2-1
     b(i)=0.25*((conjg(d(nx/2-i))+d(i)) -(sin(2.0*pi*real(i)/real(nx))+img*cos(2.0*pi*real(i)/real(nx))) *(d(i)-conjg(d(nx/2-i))))
  end do

end subroutine

subroutine rotate_calc( nx, csign, prim_fact, omega, omegan )
!  FFT に使用する回転行列を計算する.
  use Math_Const
  implicit none
  integer, intent(in) :: nx  ! 入力配列の要素数
  character(1), intent(in) :: csign  ! 正逆変換判定 [r=正変換, i=逆変換]
  integer, dimension(5) :: prim_fact  ! prim = x のとき設定すると
  ! そのべきで FFT を行う. prim_fact=(/a,b,c,d,e/) : 2^a*3^b*5^c*7^d*e
  complex, dimension(0:prim_fact(5)-1,0:prim_fact(5)-1), intent(inout) :: omega
  ! 余因数に対応する回転行列
  complex, dimension(0:nx-1,0:nx-1), intent(inout) :: omegan ! 回転行列
  integer :: counter, base
  integer :: i, j  ! do loop 用配列

  base=prim_fact(5)

  counter=0
  do i=1,4
     counter=counter+prim_fact(i)
  end do

!-- 回転行列を定義

  select case(csign)
  case('r')
     if(counter/=0)then
        if(base/=1)then
           do j=0,base-1
              do i=0,base-1
                 omega(i,j)=cos(2.0*pi*i*j/real(base))-img*sin(2.0*pi*i*j/real(base))
              end do
           end do
        end if
        do j=0,nx-1
           do i=0,nx-1
              omegan(i,j)=cos(2.0*pi*i*j/real(nx))-img*sin(2.0*pi*i*j/real(nx))
           end do
        end do

     else
        do j=0,nx-1
           do i=0,nx-1
              omega(i,j)=cos(2.0*pi*i*j/real(nx))-img*sin(2.0*pi*i*j/real(nx))
              omegan(i,j)=omega(i,j)
           end do
        end do
     end if


  case('i')
     if(counter/=0)then
        if(base/=1)then
           do j=0,base-1
              do i=0,base-1
                 omega(i,j)=cos(2.0*pi*i*j/real(base))+img*sin(2.0*pi*i*j/real(base))
              end do
           end do
        end if
        do j=0,nx-1
           do i=0,nx-1
              omegan(i,j)=cos(2.0*pi*i*j/real(nx))+img*sin(2.0*pi*i*j/real(nx))
           end do
        end do
     else
        do j=0,nx-1
           do i=0,nx-1
              omega(i,j)=cos(2.0*pi*i*j/real(nx))+img*sin(2.0*pi*i*j/real(nx))
              omegan(i,j)=omega(i,j)
           end do
        end do
     end if

  case default
     write(*,*) "******** ERROR : csign is bad. **********"
     write(*,*) "Stop!"
     stop
  end select

end subroutine