*表題 u-v-time-graph.f
*
*     n=2, k=0.0 の慣性振動の速度ベクトルとその成分, 圧力傾度力の
*     一周期分時間変化の図示
*
*履歴 2001/11/20 H.Taniguchi 
      
      integer xgr, ygr
      real ymin, ymax
      parameter (xgr=100, ygr=1, numv=34) !" グリッド数の設定
      parameter (ymin=-3.0, ymax=3.0)

      real t1(0:xgr), u1(0:xgr), v1(0:xgr), p1(0:xgr), y1
      real t2(0:xgr), u2(0:xgr), v2(0:xgr), p2(0:xgr), y2
      real t3(0:xgr), u3(0:xgr), v3(0:xgr), p3(0:xgr), y3
      real t4(0:xgr), u4(0:xgr), v4(0:xgr), p4(0:xgr), y4
      real uve1(1:numv, 3)
      real vve1(1:numv, 3)
      real uve2(1:numv, 3)
      real vve2(1:numv, 3)
      real uve3(1:numv, 3)
      real vve3(1:numv, 3)
      real uve4(1:numv, 3)
      real vve4(1:numv, 3)
      real y(3)

      character csgi*1

      omega = sqrt(5.)

      do i1 = 0, xgr
         y1   = 1.4
         t1(i1)=i1/35.4
         u1(i1)=(2./omega)*y1*(2.*y1**2-1.)
     &        *exp(-(1./2.)*y1**2)*sin(omega*t1(i1))
         v1(i1)=(4.*y1**2-2.)*exp( -(1./2.)*y1**2 )*cos(omega*t1(i1))
         p1(i1)=-(4./omega)*(-y1**4+5.*y1**2-2.)
     &        *exp(-(1./2.)*y1**2)*sin(omega*t1(i1))
      end do

      do ix1 = 1, numv
         uve1(ix1, 1) = 0.
         vve1(ix1, 1) = 0.
         uve1(ix1, 2) = u1(1+int((xgr-1)/(numv-1))*(ix1-1))
         vve1(ix1, 2) = v1(1+int((xgr-1)/(numv-1))*(ix1-1))
         uve1(ix1, 3) = 0.
         vve1(ix1, 3) = 0.
      end do

      y(1) = -0.1
      y(2) = 0.
      y(3) = 0.1

      do i2 = 0, xgr
         y2    =2.236
         t2(i2)=i2/35.4
         u2(i2)=(2./omega)*y2*(2.*y2**2-1.)
     &        *exp(-(1./2.)*y2**2)*sin(omega*t2(i2))
         v2(i2)=(4.*y2**2-2.)*exp( -(1./2.)*y2**2 )*cos(omega*t2(i2))
         p2(i2)= 0.
c         p2(i2)=-(4./omega)*(-y2**4+5.*y2**2-2.)
c     &        *exp(-(1./2.)*y2**2)*sin(omega*t2(i2))
      end do

      do ix1 = 1, numv
         uve2(ix1, 1) = 0.
         vve2(ix1, 1) = 0.
         uve2(ix1, 2) = u2(1+int((xgr-1)/(numv-1))*(ix1-1))
         vve2(ix1, 2) = v2(1+int((xgr-1)/(numv-1))*(ix1-1))
         uve2(ix1, 3) = 0.
         vve2(ix1, 3) = 0.
      end do

      do i3 = 0, xgr
         y3    =3.0
         t3(i3)=i3/35.4
         u3(i3)=(2./omega)*y3*(2.*y3**2-1.)
     &        *exp(-(1./2.)*y3**2)*sin(omega*t3(i3))
         v3(i3)=(4.*y3**2-2.)*exp( -(1./2.)*y3**2 )*cos(omega*t3(i3))
         p3(i3)=-(4./omega)*(-y3**4+5.*y3**2-2.)
     &        *exp(-(1./2.)*y3**2)*sin(omega*t3(i3))
       end do

      do ix1 = 1, numv
         uve3(ix1, 1) = 0.
         vve3(ix1, 1) = 0.
         uve3(ix1, 2) = u3(1+int((xgr-1)/(numv-1))*(ix1-1))
         vve3(ix1, 2) = v3(1+int((xgr-1)/(numv-1))*(ix1-1))
         uve3(ix1, 3) = 0.
         vve3(ix1, 3) = 0.
      end do

      do i4 = 0, xgr
         y4    =0.6
         t4(i4)=i4/35.4
         u4(i4)=(2./omega)*y4*(2.*y4**2-1.)
     &        *exp(-(1./2.)*y4**2)*sin(omega*t4(i4))
         v4(i4)=(4.*y4**2-2.)*exp( -(1./2.)*y4**2 )*cos(omega*t4(i4))
         p4(i4)=-(4./omega)*(-y4**4+5.*y4**2-2.)
     &        *exp(-(1./2.)*y4**2)*sin(omega*t4(i4))
       end do

      do ix1 = 1, numv
         uve4(ix1, 1) = 0.
         vve4(ix1, 1) = 0.
         uve4(ix1, 2) = u4(1+int((xgr-1)/(numv-1))*(ix1-1))
         vve4(ix1, 2) = v4(1+int((xgr-1)/(numv-1))*(ix1-1))
         uve4(ix1, 3) = 0.
         vve4(ix1, 3) = 0.
      end do

      write(*,*) ' workstation id (i) ? ;'
      call sgpwsn
      read(*,*) iws

*--------- 出力装置オープン      
      call sgopn(-abs(iws))

*--------- コーナマーカを消す
      call sglset('lcorner', .FALSE.)

*--------- 第一レベルのマージン
      call slmgn( 0., 0., 0., 0. )

*--------- 画面分割 (縦'T' 横'Y', 縦の分割数, 横の分割数)
      call sldiv( 'S', 1, 1 )
      call slmgn( 0., 0., 0., 0. )
      call sgsvpt( 0., 0., 0., 0. )

*--------- 1フレーム目の設定
      call grfrm
      call glrget('RUNDEF',RUNDEF)
      call grswnd( RUNDEF, RUNDEF, ymin, ymax )
      call grsvpt( 0.15, 0.90, 0.73, 0.92 )
c      call grsvpt( 0.15, 0.90, 0.76, 0.92 )
*--------- 書きたいデータの中から最大値と最小値をみつける
      call usspnt(xgr,t3,u3)
      call usspnt(xgr,t3,v3)
      call usspnt(xgr,t3,p3)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 上下の座標軸を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('B')
      call usxaxs('T')
*--------- 左右の座標軸を描く
      call usyaxs('l')
      call usyaxs('r')

*--------- X軸のタイトル指定
      call uzfact( 0.7 )
c      call uxpttl('t',1,'LAT 2.4',0.0)
c      call uxsttl('b','Time',0.0)
*--------- Y軸のタイトル指定
      call uysttl('l','Amp.',0.0)
      call uzfact( 1/0.7 )

*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( 'u' )
*--------- 折れ線を描く
      call uulinz(xgr,t3,u3,1,3)
*--------- 折れ線のラベルの設定
      call sgsplc( 'v' )
      call uulinz(xgr,t3,v3,1,1)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )
*--------- 折れ線のラベルの設定
c      call sgsplc( 'p' )
      call uulinz(xgr,t3,p3,3,1)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )

*--------- 1.5フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grsvpt( 0.15, 0.90, 0.885, 0.995 )
      call grswnd( t1(0), t1(xgr), y(1), y(3))
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call uwsgxb( t1(0), t1(xgr), numv )
      call uwsgya( y, 3 )

*--------- ベクトル図を描くルーチン
      call uglset( 'LMSG', .FALSE.)
      call uglset( 'LUMSG', .FALSE.)
      call ugiset( 'ICENT', 0) !" ベクトルの始点の位置設定
      call ugvect( uve3, numv, vve3, numv, numv, 3 )

*--------- 2フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grswnd( RUNDEF, RUNDEF, ymin, ymax )
      call grsvpt( 0.15, 0.90, 0.50, 0.69 )
c      call grsvpt( 0.15, 0.90, 0.52, 0.68 )
*--------- 書きたいデータの中から最大値と最小値をみつける
      call usspnt(xgr,t2,u2)
      call usspnt(xgr,t2,v2)
      call usspnt(xgr,t2,p2)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 上下の座標軸を描く
      call usxaxs('B')
      call usxaxs('T')
*--------- 左右の座標軸を描く
      call usyaxs('l')
      call usyaxs('r')

*--------- X軸のタイトル指定
      call uzfact( 0.7 )
c      call uxpttl('t',1,'LAT '//csgi(210)//'3',0.0)
c      call uxsttl('b','Time',0.0)
*--------- Y軸のタイトル指定
      call uysttl('l','Amp.',0.0)
      call uzfact( 1/0.7 )

*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( 'u' )
*--------- 折れ線を描く
      call uulinz(xgr,t2,u2,1,3)
*--------- 折れ線のラベルの設定
      call sgsplc( 'v' )
      call uulinz(xgr,t2,v2,1,1)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )
*--------- 折れ線のラベルの設定
c      call sgsplc( 'p' )
      call uulinz(xgr,t2,p2,3,1)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )

*--------- 2.5フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grsvpt( 0.15, 0.90, 0.66, 0.75 )
      call grswnd( t1(0), t1(xgr), y(1), y(3))
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call uwsgxb( t1(0), t1(xgr), numv )
      call uwsgya( y, 3 )

*--------- ベクトル図を描くルーチン
      call uglset( 'LMSG', .FALSE.)
      call uglset( 'LUMSG', .FALSE.)
      call ugiset( 'ICENT', 0) !" ベクトルの始点の位置設定
      call ugvect( uve2, numv, vve2, numv, numv, 3 )

*--------- 3フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grswnd( RUNDEF, RUNDEF, ymin, ymax)
      call grsvpt( 0.15, 0.90, 0.27, 0.46 )
*--------- 書きたいデータの中から最大値と最小値をみつける
      call usspnt(xgr,t1,u1)
      call usspnt(xgr,t1,v1)
      call usspnt(xgr,t1,p1)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 上下の座標軸を描く
      call usxaxs('B')
      call usxaxs('T')
*--------- 左右の座標軸を描く
      call usyaxs('l')
      call usyaxs('r')

*--------- X軸のタイトル指定
*--------- トップのタイトル指定
      call uzfact( 0.7 )
c      call uxpttl('t',1,'LAT 1',0.0)
c      call uxsttl('b','Time',0.0)
*--------- Y軸のタイトル指定
      call uysttl('l','Amp.',0.0)
      call uzfact( 1/0.7 )

*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( 'u' )
*--------- 折れ線を描く
      call uulinz(xgr,t1,u1,1,3)
*--------- 折れ線のラベルの設定
      call sgsplc( 'v' )
      call uulinz(xgr,t1,v1,1,1)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )
*--------- 折れ線のラベルの設定
c      call sgsplc( 'p' )
      call uulinz(xgr,t1,p1,3,1)

*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )

*--------- 3.5フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grsvpt( 0.15, 0.90, 0.43, 0.53 )
      call grswnd( t1(0), t1(xgr), y(1), y(3))
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call uwsgxb( t1(0), t1(xgr), numv )
      call uwsgya( y, 3 )

*--------- ベクトル図を描くルーチン
      call uglset( 'LMSG', .FALSE.)
      call uglset( 'LUMSG', .FALSE.)
      call ugiset( 'ICENT', 0) !" ベクトルの始点の位置設定
      call ugvect( uve1, numv, vve1, numv, numv, 3 )

*--------- 4フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grswnd( RUNDEF, RUNDEF, ymin, ymax)
      call grsvpt( 0.15, 0.90, 0.04, 0.23 )
*--------- 書きたいデータの中から最大値と最小値をみつける
      call usspnt(xgr,t4,u4)
      call usspnt(xgr,t4,v4)
      call usspnt(xgr,t4,p4)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 上下の座標軸を描く
      call uzlset('LABELXB', .TRUE.)
      call usxaxs('B')
      call usxaxs('T')
*--------- 左右の座標軸を描く
      call usyaxs('l')
      call usyaxs('r')

*--------- X軸のタイトル指定
*--------- トップのタイトル指定
      call uzfact( 0.7 )
c      call uxpttl('t',1,'LAT 1',0.0)
      call uxsttl('b','Time',0.0)
*--------- Y軸のタイトル指定
      call uysttl('l','Amp.',0.0)
      call uzfact( 1/0.7 )

*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( 'u' )
*--------- 折れ線を描く
      call uulinz(xgr,t4,u4,1,3)
*--------- 折れ線のラベルの設定
      call sgsplc( 'v' )
      call uulinz(xgr,t4,v4,1,1)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )
*--------- 折れ線のラベルの設定
c      call sgsplc( 'p' )
      call uulinz(xgr,t4,p4,3,1)


*--------- 4.5フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grsvpt( 0.15, 0.90, 0.20, 0.30 )
      call grswnd( t1(0), t1(xgr), y(1), y(3))
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call uwsgxb( t1(0), t1(xgr), numv )
      call uwsgya( y, 3 )

*--------- ベクトル図を描くルーチン
      call uglset( 'LMSG', .FALSE.)
      call uglset( 'LUMSG', .FALSE.)
      call ugiset( 'ICENT', 0) !" ベクトルの始点の位置設定
      call ugvect( uve4, numv, vve4, numv, numv, 3 )

*--------- ベクトルが 45 度の方向を向く時間
      call sglnzv( 0.31, 0.04, 0.22, 0.92, 1) 
      call sglnzv( 0.68, 0.04, 0.59, 0.92, 1) 
      call sgtxzv( 0.74, 0.07, 'LAT 0.6',0.02,0,-1,3) 
      call sgtxzv( 0.74, 0.30, 'LAT 1.4',0.02,0,-1,3) 
      call sgtxzv( 0.74, 0.53, 'LAT '//csgi(210)//'5',0.02,0,-1,3) 
      call sgtxzv( 0.74, 0.76, 'LAT 3.0',0.02,0,-1,3) 

*---------出力装置のクローズ
      call grcls
      end