!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  T2PACK/ISPACK-F90 テストプログラム
!      渦度から流線を求める. 
!
!履歴  2002/01/25  竹広真一
!
program t2pack_test

  use t2pack

  implicit none
  integer, parameter  :: im=32, jm=32        ! 格子点の設定(X,Y)
  integer, parameter  :: km=10,  lm=21       ! 切断波数の設定(X,Y)
  real(8), parameter  :: xmin=0,  xmax=1    ! X 座標範囲
  real(8), parameter  :: ymin=0, ymax=1     ! Y 座標範囲

  real(8), parameter  :: pi=3.1415926535897932385D0
  real(8), parameter  :: eps=1D-10

  real(8), dimension(0:jm,0:im-1) :: gg_zeta     ! 格子点データ
  real(8), dimension(0:jm,0:im-1) :: gg_zeta1    ! 格子点データ
  real(8), dimension(0:jm,0:im-1) :: gg_psi      ! 格子点データ
  real(8), dimension(0:jm,0:im-1) :: gg_psi1     ! 格子点データ
  real(8), dimension(-km:km,0:lm) :: et_zeta     ! スペクトルデータ
  real(8), dimension(-km:km,0:lm) :: et_psi      ! スペクトルデータ

  integer :: k, l, i, j

  call t2initial(im,jm,km,lm,xmin,xmax,ymin,ymax)

  ! 応力なし条件用テスト
  !gg_zeta = cos(4*pi*gg_x)*sin(pi*gg_y)
  !gg_zeta = cos(2*pi*gg_x)*sin(pi*gg_y)
  !gg_psi1  = -1/(5*pi**2)*cos(2*pi*gg_x)*sin(pi*gg_y)    ! 正解
  !gg_zeta  = sin(pi*gg_y)

  !gg_zeta = cos(4*pi*gg_x)*sin(2*pi*gg_y)
  !gg_psi1  = -1/(20*pi**2)*cos(4*pi*gg_x)*sin(2*pi*gg_y)    ! 正解

  ! 粘着条件用テスト
  gg_psi1   = cos(2*pi*gg_x)*gg_y**2*(1-gg_y)**2     ! 正解
  gg_zeta   = cos(2*pi*gg_x)*( -4*pi**2*gg_y**2*(1-gg_y)**2 &
                               + 2*(1 - 6*gg_y + 6*gg_y**2 ))
  !gg_zeta   = gg_et(et_lapla_et(et_gg(gg_psi1)))

  ! δ応答テスト
  !gg_zeta = 0
  !gg_zeta(jm/2,:) = cos(2*pi*g_x) 

  et_zeta = et_gg(gg_zeta)

  !et_psi = et_laplainv_et(et_zeta)
  !et_psi = et_vor2strm_et(et_zeta,rigid=(/.false.,.false./))
  !et_psi = et_vor2strm_et(et_zeta,rigid=(/.true.,.true./))
  et_psi = et_vor2strm1_et(et_zeta,rigid=(/.false.,.false./))
  !et_psi = et_vor2strm1_et(et_zeta,rigid=(/.true.,.true./))
  gg_psi = gg_et(et_psi)
  gg_zeta1 = gg_et(et_lapla_et(et_psi))

  gg_zeta = gg_et(et_zeta)

  call gropn(1)
  call grfrm
  call grsvpt( 0.2, 0.8, 0.2, 0.8 )
  call grswnd( real(xmin), real(xmax), real(ymin), real(ymax) )
  call grstrf
  call uwsgxa(real(g_x),im)
  call uwsgya(real(g_y),jm+1)
  call usdaxs
  call udcntr(transpose(real(gg_psi)),im,im,jm+1)

  call grfrm
  call grsvpt( 0.2, 0.8, 0.2, 0.8 )
  call grswnd( real(ymin), real(ymax), real(minval(gg_psi(:,im/2))), real(maxval(gg_psi(:,im/2))) )
  call grstrf
  call uwsgya(real(g_y),jm+1)
  call usgrph(jm+1,real(g_y),real(gg_psi(:,im/2)))

  call grfrm
  call grsvpt( 0.2, 0.8, 0.2, 0.8 )
  call grswnd( real(xmin), real(xmax), real(ymin), real(ymax) )
  call grstrf
  call uwsgxa(real(g_x),im)
  call uwsgya(real(g_y),jm+1)
  call usdaxs
  call udcntr(transpose(real(gg_psi1)),im,im,jm+1)

  call grfrm
  call grsvpt( 0.2, 0.8, 0.2, 0.8 )
  call grswnd( real(xmin), real(xmax), real(ymin), real(ymax) )
  call grstrf
  call uwsgxa(real(g_x),im)
  call uwsgya(real(g_y),jm+1)
  call usdaxs
  call udcntr(transpose(real(gg_zeta)),im,im,jm+1)

  call grfrm
  call grsvpt( 0.2, 0.8, 0.2, 0.8 )
  call grswnd( real(xmin), real(xmax), real(ymin), real(ymax) )
  call grstrf
  call uwsgxa(real(g_x),im)
  call uwsgya(real(g_y),jm+1)
  call usdaxs
  call udcntr(transpose(real(gg_zeta1)),im,im,jm+1)

  call grcls

end program t2pack_test
