*********************************************************************
* 以下では, 変換法を用いた非線形移流方程式のスペクトル計算を行なう
* 変換法で計算する手順は, 以下のとおりである.
* 1) 初期値をフーリエ変換し, スペクトルにする
* 2) スペクトルを用いて実空間での値とその空間微分の格子点値を求める
* 3) 各格子点値から, 非線形項 F_m を求める
* 4) F_m から常微分方程式の具体的な形が得られたので, 常微分方程式を
*    計算する
* 5) 得られたスペクトルから再び格子点値を計算し,... の繰り返し
*********************************************************************

      program spectrum
      integer pmax, mmax, nmax, lmax
      parameter( mmax=42 )    ! 波数空間の格子点数
      parameter( nmax=2*mmax )   ! 実空間の格子点数
      parameter( pmax=3*mmax+1 ) ! 変換法で使用する格子点数
      parameter( lmax=300 )   ! 時間ステップの回数
      real x(0:nmax-1), y(0:nmax-1)
      real xmin, xmax, dx
      real dt, s
      real pi
      complex i
      complex a(lmax,-mmax:mmax), W, V
*      complex fs(lmax,-mmax:mmax)
      complex gfs(-mmax:mmax)
      complex gts(-mmax:mmax)
      complex dwr(0:pmax-1)
      complex fwr(0:pmax-1)
      integer iws
*-- 出力ファイルをダンプ (保存) するための命令 ---
      CALL SWpSET( 'LDUMP', .TRUE. )
*      CALL SWpSET( 'LWAIT', .FALSE. )
      CALL SWpSET( 'LWAIT1', .FALSE. )
      i=(0.0,1.0)
      pi=3.14159
*-- 描画領域の決定 ---
      xmin=0.0
      xmax=2.0*pi
      dx=(xmax-xmin)/(nmax-1)
*-- 実空間での初期値 ---
      do 10 j=0,nmax-1
         x(j)=xmin+dx*j
         y(j)=exp(-(x(j)-pi)**2)
  10  continue
*-- 時間ステップの決定 ---
*-- 実空間での初期値の最大値を求める手続き ---
      s=y(0)
      do 11 j=1,nmax-1
         if(y(j).gt.s)then
            s=y(j)
         end if
  11  continue
      dt=1.0/(abs(2.0*s*mmax))
*-- フーリエ変換 ---
      CALL DFT(nmax-1,2*mmax,y(0:nmax-1),gfs(-mmax:mmax))
      do 25 k=-mmax,mmax
         a(1,k)=gfs(k)
  25  continue
*-- 格子点値の計算 ---
*-- 1 回目は, 空間の値は初期値を用いればよいので. ---
*-- 空間微分の格子点のみ求める ---
      CALL CIDFT(pmax-1,2*mmax,gfs(-mmax:mmax),fwr(0:pmax-1))
      do m=-mmax,mmax
         gfs(m)=real(m)*gfs(m)
      end do
      CALL CIDFT(pmax-1,2*mmax,gfs(-mmax:mmax),dwr(0:pmax-1))
      do j=-mmax,mmax
         gfs(j)=(0.0,0.0)
      end do
      do j=0,pmax-1
         fwr(j)=dwr(j)*fwr(j)
         dwr(j)=(0.0,0.0)
      end do
*-- F_m を求める ---
      CALL CDFT(pmax-1,2*mmax,fwr(0:pmax-1),gfs(-mmax:mmax))
      do j=0,pmax-1
         fwr(j)=(0.0,0.0)
      end do
*-- 常微分方程式 --
*-- 第 1 ステップは時間半分のオイラー積分 ---
      do 31 k=-mmax,mmax
         a(2,k)=a(1,k)-i*dt*gfs(k)
         gfs(k)=a(2,k)
  31  continue
*-- 第 2 ステップからリープフロッグ法 ---
      do 30 j=2,lmax-1
*-- 変換法ここから ---
         CALL CIDFT(pmax-1,2*mmax,gfs(-mmax:mmax),fwr(0:pmax-1))
         do m=-mmax,mmax
            gfs(m)=real(m)*a(j,m)
         end do
         CALL CIDFT(pmax-1,2*mmax,gfs(-mmax:mmax),dwr(0:pmax-1))
         do k=0,pmax-1
            fwr(k)=dwr(k)*fwr(k)
            dwr(k)=(0.0,0.0)
         end do
         do k=-mmax,mmax
            gfs(k)=(0.0,0.0)
         end do
         CALL CDFT(pmax-1,2*mmax,fwr(0:pmax-1),gfs(-mmax:mmax))
         do k=0,pmax-1
            fwr(k)=(0.0,0.0)
         end do
*         do m=-mmax,mmax
*            fs(j,m)=fts(m)
*         end do
         do 35 k=-mmax,mmax
            a(j+1,k)=a(j-1,k)-2.0*i*dt*gfs(k)
            gfs(k)=a(j+1,k)
*            a(j+1,k)=a(j-1,k)-2.0*dt*fs(j,k)
  35  continue
  30  continue
*-- DFT の計算 ---
*-- グラフへの描画 ---
*-- グラフ ---
      write(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
      read(*,*) IWS

      CALL GROPN( IWS )

      do 200 l=1,lmax
         if(mod((l-1),10).eq.0)then
         do 210 j=-mmax,mmax
            gts(j)=a(l,j)
  210 continue
         CALL IDFT(nmax-1,2*mmax,gts(-mmax:mmax),y(0:nmax-1))
         CALL GRFRM

         CALL USSPNT( nmax, x(0:nmax-1), y(0:nmax-1) )
         CALL GRSWND( xmin, xmax, 0.0, 1.0 )
*         CALL USPFIT
         CALL GRSVPT( 0.2, 0.8, 0.2, 0.8 )
         CALL GRSTRF
         CALL USDAXS
         CALL UULIN( nmax, x(0:nmax-1), y(0:nmax-1) )
      end if
  200 continue
      CALL GRCLS
      stop
      end