*表題  赤道波の分散曲線の理論解を表示するプログラム 2
*      (3次方程式の解の公式から導いた理論解を表示)
*
*履歴  2001/03/19 H.Taniguchi

      parameter( NMAX = 50 )
      integer n
      real i
      real omega01(NMAX)
      real omega02(NMAX)
      real omega03(NMAX)
      real omegan1(NMAX)
      real omegan2(NMAX)
      real omegan3(NMAX)
      real k
      real a
      real theta
      real kx(NMAX)
      character csgi*1
      character aaa
      character title00*3
      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.)

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

      n = 0

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

      in = 0

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

         a = sqrt( (1./3.)*(k**2+2.*n+1.) )
         theta
     &     = (1./3.)*(atan(sqrt((4./27.)*(k**2+2.*n+1.)**3-k**2)/k ))
         
         omega01(in) = 2.*a*cos(theta)
         omega02(in) = a*(-cos(theta)-sqrt(3.)*sin(theta))
         omega03(in) = a*(-cos(theta)+sqrt(3.)*sin(theta))

         write(*,*) 'kx(',in,')=',kx(in)
         write(*,*) 'omega01(',in,')=',omega01(in)
         write(*,*) 'omega02(',in,')=',omega02(in)
         write(*,*) 'omega03(',in,')=',omega03(in)

      end do

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

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

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


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

      do n = 1,3

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

         in = 0

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

         a = sqrt( (1./3.)*(k**2+2*n+1) )
         theta
     &        = (1./3.)*(atan(sqrt((4./27.)*(k**2+2*n+1)**3-k**2)/k ))
         
         omegan1(in) = 2.*a*cos(theta)
         omegan2(in) = a*(-cos(theta)-sqrt(3.)*sin(theta))
         omegan3(in) = a*(-cos(theta)+sqrt(3.)*sin(theta))

         write(*,*) 'kx(',in,')=',kx(in)
         write(*,*) 'omegan01(',in,')=',omega01(in)
         write(*,*) 'omegan02(',in,')=',omega02(in)
         write(*,*) 'omegan03(',in,')=',omega03(in)

         end do

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

*--------- ラベル文字列の文字高設定
      call sgspls(0.01)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( title1 )
*--------- 折れ線の線種の設定
      call uuslnt(1)
*--------- 分散曲線を描く
      call uulin( NMAX, kx, omegan1 )
      call uulin( NMAX, kx, omegan2 )
      call uulin( NMAX, kx, omegan3 )
*--------- ラベル文字列の消去
      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
