!----------------------------------------------------------------------
!     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=16        ! 格子点の設定(X,Y)
  integer, parameter  :: km=10,  lm=10        ! 切断波数の設定(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_data     ! 格子点データ
  real(8), dimension(0:jm,0:im-1) :: gg_data1    ! 格子点データ
  real(8), dimension(0:jm,0:im-1) :: gg_data2    ! 格子点データ
  real(8), dimension(-km:km,0:lm) :: et_data     ! スペクトルデータ

  integer :: k, l, i, j

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

  gg_data1 = cos(4*pi*gg_x)*sin(pi*gg_y)
  gg_data  = - 1/(17*pi**2) * cos(4*pi*gg_x)*sin(pi*gg_y)
  !gg_data1 = sin(2*pi*gg_x)*sin(pi*gg_y)
  !gg_data  = - 1/(5*pi**2) * sin(2*pi*gg_x)*sin(pi*gg_y)
  !gg_data1 = sin(pi*gg_y)
  !gg_data  = - 1/(pi**2) * sin(pi*gg_y)

  et_data = et_gg(gg_data1)
  gg_data2 = gg_et(et_laplainv_et(et_data)) 

  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_data)),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_data2)),im,im,jm+1)

  call grcls

  gg_data2 = gg_data2 - gg_data
  write(6,*)'Laplacian-inverse'
  do j=0,jm
     do i=0,im-1
        if ( abs(gg_data2(j,i)) > eps ) write(6,*) i,j, gg_data2(j,i)
     enddo
  enddo

end program t2pack_test
