* PACKAGE TTFFT99    !" テスト : フーリエ変換
*
*"  [HIS] 95/08/09 (takepiro) 
*
***********************************************************************
      PROGRAM TFFT99
*
* [DIMENSION]
      INTEGER    NX                        !" X方向の格子点数
      INTEGER    NZ                        !" Z方向の格子点数
	PARAMETER( NX = 8, NZ = 6 )
*
*
* [INPUT]
      REAL       ZETA ( 0:NX+1 , 0:NZ+1 )  !" 渦度
*
* [WORK]
      REAL Z    ( NX , NZ )	           !" 作業配列
      REAL WTBL ( 2, NX )		   !" 作業配列
*
* [INTERNAL WORK]
      INTEGER IX
      INTEGER IZ
      REAL    PI 		           !" 円周率
*
      PI =  ATAN2( 0.0, -1.0 )
*
*"   < 1. 渦度分布の設定 >
*
      DO 100 IZ = 0, NZ+1
        DO 100 IX = 0, NX+1
           ZETA( IX, IZ ) = COS( 2 * PI * (IZ-1) * (IX-1) / REAL(NX) )
  100 CONTINUE
*
      WRITE ( 6, * ) ' Before FFT'
      DO 110 IZ=1,NZ
         WRITE( 6, '(8F9.3)' ) (ZETA(IX,IZ), IX=1,NX)
  110 CONTINUE
*
*"     < 2. 三角関数表作成 >
*
      CALL RFFTIM     !" FFT 3角関数表作成
     I         ( NX    ,
     O           WTBL  , IFAC  )
*
*"     < 3. フーリエ正変換 >
*
      CALL FFT99Y     !" FFT を呼ぶ
     M         ( ZETA  ,
     M           Z     ,
     C           WTBL  , 
     C           2     , NX  , NX+2  , 
     C           2     , NZ  , NZ+2  , 
     C           0                      )
*
*"   < 4. 出力 >
*
      WRITE ( 6, * ) ' After FFT'
      DO 300 IZ=1,NZ
         WRITE( 6, '(8F9.3)' ) (Z(IX,IZ), IX=1,NX)
  300 CONTINUE
*
*"     < 5. フーリエ逆変換 >
*
      CALL FFT99Y     !" FFT を呼ぶ
     M         ( ZETA  ,
     M           Z     ,
     C           WTBL  , 
     C           2     , NX  , NX+2  , 
     C           2     , NZ  , NZ+2  , 
     C           1                      )
*
*"     < 6. 出力 >
*
      WRITE ( 6, * ) ' After Inverse FFT'
      DO 600 IZ=1,NZ
         WRITE( 6, '(8F9.3)' ) (ZETA(IX,IZ), IX=1,NX) 
  600 CONTINUE
*
      STOP
      END
