!----------------------------------------------------------------------
!     Copyright (c) 2010--2012 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wtq_module テストプログラム
!
!      磁場トロイダルポテンシャルの境界値問題
!
!履歴  2010/04/18  竹広真一  wtq_test_polmagbcgrid.f90 より改変
!      2012/04/03  竹広真一  wtq_module_sjpack_polmagbcgrid_test.f90 より改変
!
program wtq_module_polmagbcgrid_test

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use wtq_module

  implicit none

  integer,parameter  :: im=32, jm=16         ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: kmo=16, kmi=16       ! 格子点の設定(球殻動径, 球動径)
  integer,parameter  :: nm=10                ! 切断波数の設定(水平)
  integer,parameter  :: lmo=10, lmi=21       ! 切断波数の設定(球殻動径, 球動径)
  real(8),parameter  :: ri=0.5, ro=1.5       ! 内外半径

  real(8), dimension(im,jm,0:kmo)           :: xyz_POLMAG
  real(8), dimension((nm+1)*(nm+1),0:kmo)   :: wz_POLMAG
  real(8), dimension(im,jm,kmi)             :: xyr_POLMAG
  real(8), dimension((nm+1)*(nm+1),kmi)     :: wr_POLMAG
  real(8), dimension(im,jm,0:kmo)           :: xyz_DPOLDR
  real(8), dimension(im,jm,kmi)             :: xyr_DPOLDR

  real(8), dimension((nm+1)*(nm+1),0:kmo)   :: wz_TopBoundary
  real(8), dimension((nm+1)*(nm+1),0:kmo)   :: wz_n   ! 全波数
  real(8), dimension((nm+1)*(nm+1))         :: w_Null=0.0D0

  real(8), parameter  :: pi=3.1415926535897932385D0
  ! 判定誤差設定
  integer, parameter :: check_digits = 11
  integer, parameter :: ignore = -12

  integer :: n,k, nn(2)

  call MessageNotify('M','wtq_test_tormagbcgrid', &
       'wtq_module  wtq_TormagBoundariesGrid subroutine test')

  call wtq_Initial(im,jm,kmi,kmo,nm,lmi,lmo,ri,ro)

  do k=0,kmo
     do n=1,(nm+1)**2
        nn=nm_l(n)
        wz_n(n,k) = nn(1)
     enddo
  enddo

  ! P_10

  xyz_POLMAG = sin(xyz_lat) * cos( pi*(xyz_rad-ri)/(ro-ri) )
  xyr_POLMAG = sin(xyr_lat) * cos( pi*(xyr_rad-ri)/ri ) * xyr_Rad

  wz_POLMAG = wz_xyz(xyz_POLMAG)
  wr_POLMAG = wr_xyr(xyr_POLMAG)
  call wtq_PolmagBoundariesGrid(wz_POLMAG,wr_POLMAG)
  xyz_POLMAG = xyz_wz(wz_POLMAG)
  xyr_POLMAG = xyr_wr(wr_POLMAG)

  xyz_DPOLDR = xyz_wt(wt_DRad_wt(wt_wz(wz_POLMAG)))
  xyr_DPOLDR = xyr_wq(wq_RadDRad_wq(wq_wr(wr_POLMAG)))/xyr_Rad
  wz_TopBoundary = wz_wt(wt_DRad_wt(wt_wz(wz_POLMAG))) &
                     + (wz_n +1)*wz_POLMAG/wz_RAD

  call checkresults

  ! P_1_1
  xyz_POLMAG = cos(xyz_lat)*cos(xyz_lon)* cos( pi*(xyz_rad-ri)/(ro-ri) )
  xyr_POLMAG = cos(xyr_lat)*cos(xyr_lon)* cos( pi*(xyr_rad-ri)/ri )*xyr_Rad

  wz_POLMAG = wz_xyz(xyz_POLMAG)
  wr_POLMAG = wr_xyr(xyr_POLMAG)
  call wtq_PolmagBoundariesGrid(wz_POLMAG,wr_POLMAG)
  xyz_POLMAG = xyz_wz(wz_POLMAG)
  xyr_POLMAG = xyr_wr(wr_POLMAG)

  xyz_DPOLDR = xyz_wt(wt_DRad_wt(wt_wz(wz_POLMAG)))
  xyr_DPOLDR = xyr_wq(wq_RadDRad_wq(wq_wr(wr_POLMAG)))/xyr_Rad

  call checkresults

  call MessageNotify('M','wtq_test_polmagbcgrid', &
       'wtq_module  wtq_TormagBoundariesGrid subroutine test succeded')


contains
  subroutine checkresults

    call AssertEqual(&
      message='Top B.C',                                            &
      answer = wz_TopBoundary(:,0),                                 &
      check = w_Null,                                               &
      significant_digits = check_digits, ignore_digits = ignore     &
      )

    call AssertEqual(&
      message='Inner B.C (value conituity) ',                       &
      answer = xyz_POLMAG(:,:,kmo),                                 &
      check = xyr_POLMAG(:,:,kmi),                                  &
      significant_digits = check_digits, ignore_digits = ignore     &
      )

    call AssertEqual(&
      message='Inner B.C (derivative conituity) ',                  &
      answer = xyz_DPOLDR(:,:,kmo),                                 &
      check = xyr_DPOLDR(:,:,kmi),                                  &
      significant_digits = check_digits, ignore_digits = ignore     &
      )
  end subroutine checkresults

end program wtq_module_polmagbcgrid_test