*表題  赤道波の分散曲線の理論解を表示するプログラム
*
*履歴  2001/03/19 H.Taniguchi

      parameter( NMAX = 51 )
      integer n
      real i
      real omega0a(NMAX)
      real omega0b(NMAX)
      real omega0c(NMAX)
      real omega1(NMAX)
      real omega2(NMAX)
      real omega3(NMAX)
      real k
      real kx(NMAX)
      character csgi*1
      character aaa
      character title0*3
      character title1*3

      title0 = 'n= '
      title1 = 'n= '

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

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

      call grfrm
      call grsvpt(  0.2,  0.8,  0.2,  0.8 )

*--------- ウインドウの設定
      call glrget('RUNDEF',RUNDEF)

      write(*,*) '横軸(東西波数)の最小値いくつにしますか?'
      read(*,*) xkmin
      write(*,*) '横軸(東西波数)の最大値いくつにしますか?'
      read(*,*) xkmax
      write(*,*) '縦軸(振動数)の最小値いくつにしますか?'
      read(*,*) ycmin
      write(*,*) '縦軸(振動数)の最大値いくつにしますか?'
      read(*,*) ycmax

      call grswnd( xkmin, xkmax, ycmin, ycmax )

*--------- 枠内にグラフを納めるための設定(クリッピング)
      call sglset('LCLIP',.TRUE.)

*---- n＝0 の分散関係式

      n = 0

      write(title0(3:3),'(i1)') n

      in = 0
         
      do i  = 0, 5, 0.1
         in = in + 1
         k  = real(i)
         kx(in) = k

         omega0a(in) =  0.5*( k + sqrt( k**2 + 4 ) )
         omega0b(in) =  0.5*( k - sqrt( k**2 + 4 ) )
         omega0c(in) = k

         write(*,*) 'kx(',in,')=',kx(in)
         write(*,*) 'omega0a(',in,')=',omega0a(in)
         write(*,*) 'omega0b(',in,')=',omega0b(in)         
      end do

*--------- 書きたいデータの中から最大値と最小値をみつける
      call usspnt(NMAX,kx,omega0a)
      call usspnt(NMAX,kx,omega0b)
      call usspnt(NMAX,kx,omega0c)

*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 折れ線の線種の設定
c      call uuslnt(5)
*--------- 分散曲線を描く
      call uulin( NMAX, kx, omega0c )
*--------- ラベル文字列の文字高設定
      call sgspls(0.01)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( title0 )
*--------- 折れ線の線種の設定
      call uuslnt(2)
*--------- 分散曲線を描く
      call uulin( NMAX, kx, omega0a )
*--------- 折れ線の線種の設定
      call uuslnt(4)
*--------- 分散曲線を描く
      call uulin( NMAX, kx, omega0b )

*---- n≠0 の分散関係式

      do n = 1,3

         write(title1(3:3),'(i1)') n

         in = 0

         do i  = 0, 5, 0.1
            in = in + 1
            k  = real(i)
            kx(in) = k

            omega1(in) =  sqrt(k**2 + 2*n + 1)
            omega2(in) = -sqrt(k**2 + 2*n + 1)
            omega3(in) = -k/(k**2 + 2*n + 1)

         write(*,*) 'omega1(',in,')=',omega1(in)
         write(*,*) 'omega2(',in,')=',omega2(in)         
         write(*,*) 'omega3(',in,')=',omega3(in)         

         end do

*--------- 書きたいデータの中から最大値と最小値をみつける
      call usspnt(NMAX,kx,omega1)
      call usspnt(NMAX,kx,omega2)
      call usspnt(NMAX,kx,omega3)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- ラベル文字列の文字高設定
      call sgspls(0.01)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( title1 )
*--------- 折れ線の線種の設定
      call uuslnt(1)
*--------- 分散曲線を描く
      call uulin( NMAX, kx, omega1 )
      call uulin( NMAX, kx, omega2 )
*--------- 折れ線の線種の設定
      call uuslnt(3)
      call uulin( NMAX, kx, omega3 )
*--------- ラベル文字列の消去
      call sglset( 'LCHAR' , .FALSE. )
*--------- 枠内にグラフを納めるための設定(クリッピング)
      call sglset('LCLIP',.TRUE.)

      end do

*--------- x 座標の間隔指定
      call usxaxs('B')
      call usxaxs('T')
*--------- y 座標の間隔指定
      call usyaxs('L')
      call usyaxs('R')
*--------- トップのタイトル指定
c      call uxsttl('t','dispersion relation',0.0)
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
      call uxsttl('b','zonal wavenumber k',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
      call uysttl('l','frequency '//csgi(175), 0.0)
*--------- ラベル文字列の消去
      call sglset( 'LCHAR' , .FALSE. )

      call grcls
      stop
      end
