program spectrum
      parameter( nmax=50 )   ! 実空間の格子点数
      parameter( mmax=50 )    ! 波数空間の格子点数
      parameter( lmax=1000 )   ! 時間ステップの回数
      real x(0:2*nmax-1), y(0:2*nmax-1)
      real xmin, xmax, dx
      real dt
      real pi, gamma
      complex i
      complex a(lmax,-mmax:mmax), W, V
      complex b(0:2*nmax-1), gf(-mmax:mmax)
      integer iws
*-- 出力ファイルをダンプ (保存) するための命令 ---
      CALL SWpSET( 'LDUMP', .TRUE. )
*      CALL SWpSET( 'LWAIT', .FALSE. )
      CALL SWpSET( 'LWAIT1', .FALSE. )
      i=(0.0,1.0)
      pi=3.14159
      gamma=1.0                 ! 移流速度
*-- 時間ステップの決定 ---
      dt=1.0/(mmax*mmax*gamma)
*-- 描画領域の決定 ---
      xmin=0.0
      xmax=2.0*pi
      dx=(xmax-xmin)/(2*nmax-1)
*-- 実空間での初期値 ---
      do 10 j=0,2*nmax-1
         x(j)=xmin+dx*j
         y(j)=exp(-(x(j)-pi)**2)
  10  continue
*-- フーリエ変換 ---
      CALL DFT(2*nmax-1,2*mmax,y(0:2*nmax-1),gf(-mmax:mmax))
      do 25 k=-mmax,mmax
         a(1,k)=gf(k)
  25  continue
*-- 常微分方程式 ---
      do 30 j=1,lmax-1
         do 35 k=-mmax+1,mmax-1
            a(j+1,k)=(1.0-k*k*dt)*a(j,k)
*            a(j+1,k)=((1.0-i*k*dt)/(1.0+k*k*dt*dt))*a(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),100).eq.0)then
         do 210 j=-mmax,mmax
            gf(j)=a(l,j)
  210 continue
         CALL IDFT(2*nmax-1,2*mmax,gf(-mmax:mmax),y(0:2*nmax-1))
         CALL GRFRM

         CALL USSPNT( 2*nmax, x(0:2*nmax-1), y(0:2*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( 2*nmax, x(0:2*nmax-1), y(0:2*nmax-1) )
      end if
  200 continue
      CALL GRCLS
      stop
      end