*************************************
* x,z 平面で描画し, パラメータを
* a=10.0, b=28.0, c=8/3 で設定し,
* 初期条件を x=0.0, y=4.0, z=28.0
* で入力すると綺麗な絵が描ける
*************************************
      program Lorentz
      parameter( nmax=5000 )
      real t(nmax+1), x(nmax+1), y(nmax+1), z(nmax+1)
      real dt, dx
      real x_0, y_0, z_0
      real k1, k2, k3, k4
      real j1, j2, j3, j4
      real l1, l2, l3, l4
      real tmin, tmax
      real a, b, c
      integer IW
*-- 出力ファイルをダンプ (保存) するための命令 ---
*      CALL SWpSET( 'LDUMP', .TRUE. )
*      CALL SWpSET( 'LWAIT', .FALSE. )
*      CALL SWpSET( 'LWAIT1', .FALSE. )
*-- 係数の入力方式についての指定 ---
      write(*,*) "Which do you select?"
      write(*,*) "0: type, 1: pretype"
      read(*,*) IW
*-- 係数の決定 ---
      if(IW.eq.0)then
         write(*,*) "input the parameters"
         write(*,*) "a=, b=, c="
         read(*,*) a, b, c
      else
         a=10.0
         b=28.0
         c=8.0/3.0
      end if

*-- 時間ステップの設定 ---
      dt=0.01
*-- 初期条件の設定 ---
      t(1)=0.0
      if(IW.eq.0)then
         write(*,*) "input the initialized numbers"
         write(*,*) "x_0=, y_0=, z_0="
         read(*,*) x_0, y_0, z_0 
         x(1)=x_0
         y(1)=y_0
         z(1)=z_0
      else
         x(1)=0.0
         y(1)=4.0
         z(1)=28.0
      end if
*-- 実際の計算ここから ---
      do 10 i=1,nmax
*-- k の計算 ---
         t(i+1)=t(i)+dt
         k1=a*(-x(i)+y(i))*dt
         j1=(-x(i)*z(i)+b*x(i)-y(i))*dt
         l1=(x(i)*y(i)-c*z(i))*dt
         k2=a*(-(x(i)+0.5*k1)+(y(i)+0.5*j1))*dt
         j2=(-(x(i)+0.5*k1)*(z(i)+0.5*l1)
     &      +b*(x(i)+0.5*k1)-(y(i)+0.5*j1))*dt
         l2=((x(i)+0.5*k1)*(y(i)+0.5*j1)-c*(z(i)+0.5*l1))*dt
         k3=a*(-(x(i)+0.5*k2)+(y(i)+0.5*j2))*dt
         j3=(-(x(i)+0.5*k2)*(z(i)+0.5*l2)
     &      +b*(x(i)+0.5*k2)-(y(i)+0.5*j2))*dt
         l3=((x(i)+0.5*k2)*(y(i)+0.5*j2)-c*(z(i)+0.5*l2))*dt
         k4=a*(-(x(i)+k3)+(y(i)+j3))*dt
         j4=(-(x(i)+k3)*(z(i)+l3)
     &      +b*(x(i)+k3)-(y(i)+j3))*dt
         l4=((x(i)+k3)*(y(i)+j3)-c*(z(i)+l3))*dt
*-- 時間ステップの計算 ---
         x(i+1)=x(i)+(1.0/6.0)*(k1+2.0*k2+2.0*k3+k4)
         y(i+1)=y(i)+(1.0/6.0)*(j1+2.0*j2+2.0*j3+j4)
         z(i+1)=z(i)+(1.0/6.0)*(l1+2.0*l2+2.0*l3+l4)
  10  continue

*-- グラフの出力装置の指定 (IWS=4 は GTK 指定) ----
      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
*      IWS=4
      READ (*,*) IWS

*-- 出力装置の呼び出し
      CALL GROPN( IWS )
*-- 分割指定 ---
      CALL SLDIV( 'Y', 2, 2 )

*-- 1 枚目 (x,y) 座標 ---
      CALL GRFRM
      CALL USSPNT( nmax+1, x, y )

      CALL GRSWND(  -25.0,  25.0,  -25.0,  25.0 )
      CALL GRSVPT(  0.1,  0.9,  0.1,  0.9 )
      CALL USPFIT
      CALL GRSTRF

      CALL USSTTL( 'X-AXIS', '', 'Y-AXIS', '' )
      CALL USDAXS

      CALL UULIN(nmax+1, x, y )

*-- 2 枚目 ---
      CALL GRFRM
      CALL USSPNT( nmax+1, x, z )

      CALL GRSWND(  -25.0,  25.0,  0.0,  50.0 )
      CALL GRSVPT(  0.1,  0.9,  0.1,  0.9 )
      CALL USPFIT
      CALL GRSTRF

      CALL USSTTL( 'X-AXIS', '', 'Z-AXIS', '' )
      CALL USDAXS

      CALL UULIN(nmax+1, x, z )

*-- 3 枚目 ---
      CALL GRFRM
      CALL USSPNT( nmax+1, y, z )

      CALL GRSWND(  -25.0,  25.0,  0.0,  50.0 )
      CALL GRSVPT(  0.1,  0.9,  0.1,  0.9 )
      CALL USPFIT
      CALL GRSTRF

      CALL USSTTL( 'Y-AXIS', '', 'Z-AXIS', '' )
      CALL USDAXS

      CALL UULIN(nmax+1, y, z)
      CALL GRCLS

      stop
      end