!----------------------------------------------------------------------
!     Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wsc_module_svpack テストプログラム
!
!      微分のテスト
!  
!履歴  2017/05/21  竹広真一
!
program wsc_module_svpack_derivative1

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use wsc_module_svpack
  use wa_module_svpack, only : wa_Lapla_wa
  implicit none

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

  real(8), dimension(0:im-1,1:jm,0:km)     :: xyz_Data
  real(8), dimension(0:im-1,1:jm,0:km)     :: xyz_DData

  real(8), dimension(0:im-1,1:jm,0:km)     :: xyz_VLon
  real(8), dimension(0:im-1,1:jm,0:km)     :: xyz_VLat
  real(8), dimension(0:im-1,1:jm,0:km)     :: xyz_VRad
  real(8), dimension(0:im-1,1:jm,0:km)     :: xyz_Div
  real(8), dimension(0:im-1,1:jm,0:km)     :: xyz_Psi

  real(8), dimension(0:im-1,1:jm,0:km)     ::  xyz_xi

  ! 判定誤差設定
  integer, parameter :: check_digits = 10
  integer, parameter :: ignore = -11

  integer, parameter :: nmin=1, nmax=9

  integer :: n

  real(8) :: pi

  pi = atan(1.0D0)*4.0D0

  call MessageNotify('M','wsc_module_svpack_derivative_test1', &
       'wsc_module_svpack derivative function test #1')

  call wsc_Initial(im,jm,km,nm,lm,ri,ro)

  xyz_xi = (xyz_Rad - ri)/(ro-ri)

  !==================== wc_DRad_ws, ws_DRad_wc ====================  
  xyz_Data  = sin(PI*xyz_xi)
  xyz_DData = PI*cos(PI*xyz_xi)
  !
  call AssertEqual(&
       message='wc_DRad_ws with sin(z)',                             &
       answer = xyz_DData,                                           &
       check = xyz_wc(wc_DRad_ws(ws_xyz(xyz_Data))),                 &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_Data  = cos(2*PI*xyz_xi)
  xyz_DData = -2*PI*sin(2*PI*xyz_xi)
  !
  call AssertEqual(&
       message='ws_DRad_wc with cos(2z)',                            &
       answer = xyz_DData,                                           &
       check = xyz_ws(ws_DRad_wc(wc_xyz(xyz_Data))),                 &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !==================== wc_DivRad_ws, ws_DivRad_wc ====================  

  xyz_Data = cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * cos(PI*xyz_xi) ! Y_2^-1 cos(z)
  xyz_DData = -PI* cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * sin(PI*xyz_xi)
  !
  call AssertEqual(&
       message='ws_DivRad_wc with Y_2^-1 cos(z)',                    &
       answer = xyz_DData,                                           &
       check = xyz_ws(ws_DivRad_wc(wc_xyz(xyz_Data))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_Data = cos(xyz_Lat)*sin(xyz_Lat) * cos(xyz_Lon) * sin(2*PI*xyz_xi) ! Y_2^1 sin(2z)
  xyz_DData =2*PI*cos(xyz_Lat)*sin(xyz_Lat) * cos(xyz_Lon) * cos(2*PI*xyz_xi)
  !
  call AssertEqual(&
       message='wc_DivRad_ws with Y_2^1 sin(2z)',                    &
       answer = xyz_DData,                                           &
       check = xyz_wc(wc_DivRad_ws(ws_xyz(xyz_Data))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !==================== wc_RotRad_ws, ws_RotRad_wc ====================  

  xyz_Data = sin(xyz_Lat) * cos(3*PI*xyz_xi) ! Y_2^0 cos(3z)
  xyz_DData = -3*PI* sin(xyz_Lat) * sin(3*PI*xyz_xi)
  !
  call AssertEqual(&
       message='ws_RotRad_wc with Y_2^0 cos(3z)',                    &
       answer = xyz_DData,                                           &
       check = xyz_ws(ws_RotRad_wc(wc_xyz(xyz_Data))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_Data = sin(xyz_Lat) * sin(4*PI*xyz_xi) ! Y_2^0 sin(4z)
  xyz_DData = 4*PI* sin(xyz_Lat) * cos(4*PI*xyz_xi)
  !
  call AssertEqual(&
       message='wc_RotRad_ws with Y_2^0 sin(4z)',                    &
       answer = xyz_DData,                                           &
       check = xyz_wc(wc_RotRad_ws(ws_xyz(xyz_Data))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !==================== ws_Lapla_ws, wc_Lapla_wc ====================  
  
  xyz_Data  = cos(xyz_Lat)*sin(xyz_Lon)*sin(3*PI*xyz_xi)  ! Y_1^-1 sin(3z)
  xyz_DData = -( 9*PI**2 + 1*2/Ro**2) * xyz_Data
  !
  call AssertEqual(&
       message='ws_Lapla_ws with Y_1^-1 sin(3z)',                    &
       answer = xyz_DData,                                           &
       check = xyz_ws(ws_Lapla_ws(ws_xyz(xyz_Data))),                &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_Data  = cos(xyz_Lat)*cos(xyz_Lon)*cos(4*PI*xyz_xi)  ! Y_1^1 cos(4z)
  xyz_DData = -( 16*PI**2 + 1*2/Ro**2) * xyz_Data
  !
  call AssertEqual(&
       message='wc_Lapla_wc with Y_1^1 cos(4z)',                     &
       answer = xyz_DData,                                           &
       check = xyz_wc(wc_Lapla_wc(wc_xyz(xyz_Data))),                &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !==================== ws_LaplaInv_ws, wc_LaplaInv_wc ====================  
  
  xyz_Data  = cos(xyz_Lat)*sin(xyz_Lon)*sin(3*PI*xyz_xi)  ! Y_1^-1 sin(3z)
  xyz_DData = -1.0D0/( 9*PI**2 + 1*2/Ro**2) * xyz_Data
  !
  call AssertEqual(&
       message='ws_LaplaInv_ws with Y_1^-1 sin(3z)',                 &
       answer = xyz_DData,                                           &
       check = xyz_ws(ws_LaplaInv_ws(ws_xyz(xyz_Data))),             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_Data  = cos(xyz_Lat)*sin(xyz_Lon)*sin(2*PI*xyz_xi)  ! Y_1^-1 sin(2z)
  xyz_DData = -1.0D0/( 4*PI**2 + 1*2/Ro**2) * xyz_Data
  !
  call AssertEqual(&
       message='ws_LaplaInv_ws with Y_1^-1 sin(2z)',                 &
       answer = xyz_DData,                                           &
       check = xyz_ws(ws_LaplaInv_ws(ws_xyz(xyz_Data))),             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_Data  = cos(xyz_Lat)*cos(xyz_Lon)*cos(4*PI*xyz_xi)  ! Y_1^1 cos(4z)
  xyz_DData = -1.0D0/( 16*PI**2 + 1*2/Ro**2) * xyz_Data
  !
  call AssertEqual(&
       message='wc_LaplaInv_wc with Y_1^1 cos(4z)',                  &
       answer = xyz_DData,                                           &
       check = xyz_wc(wc_LaplaInv_wc(wc_xyz(xyz_Data))),             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_Data  = cos(xyz_Lat)*cos(xyz_Lon)*cos(5*PI*xyz_xi)  ! Y_1^1 cos(4z)
  xyz_DData = -1.0D0/( 25*PI**2 + 1*2/Ro**2) * xyz_Data
  !
  call AssertEqual(&
       message='wc_LaplaInv_wc with Y_1^1 cos(5z)',                  &
       answer = xyz_DData,                                           &
       check = xyz_wc(wc_LaplaInv_wc(wc_xyz(xyz_Data))),             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !========== xyz_GradLon_ws, xyz_GradLat_ws ==========

  ! ------ Y_1^-1 sin(z) Toridal field -----

  xyz_Psi = cos(xyz_Lat)*sin(xyz_Lon) * sin(PI*xyz_xi)    ! Y_1^-1 sin(z)
  xyz_VRad = 0

  xyz_VLon =  cos(xyz_Lon) * sin(PI*xyz_xi) /Ro
  call AssertEqual(&
       message='xyz_GradLon_ws with Y_1^-1 sin(z)',                  &
       answer = xyz_VLon,                                            &
       check = xyz_GradLon_ws(ws_xyz(xyz_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_VLat = -sin(xyz_Lat)*sin(xyz_Lon) * sin(PI*xyz_xi) /Ro
  call AssertEqual(&
       message='xyz_GradLat_ws with Y_1^-1 sin(z)',                  &
       answer = xyz_VLat,                                            &
       check = xyz_GradLat_ws(ws_xyz(xyz_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !========== xyz_GradLon_wc, xyz_GradLat_wc, wc_Div_xyz_xyz_xyz, wc_RotRad_xyz_xyz ==========
  ! ------ Y_1^1 cos(2z) Toridal field -----

  xyz_Psi = cos(xyz_Lat)*cos(xyz_Lon) * cos(2*PI*xyz_xi)    ! Y_1^1 cos(2z)
  xyz_VRad = 0

  xyz_VLon =  - sin(xyz_Lon) * cos(2*PI*xyz_xi) /Ro
  call AssertEqual(&
       message='xyz_GradLon_wc with Y_1^1 cos(2z)',                  &
       answer = xyz_VLon,                                            &
       check = xyz_GradLon_wc(wc_xyz(xyz_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_VLat = -sin(xyz_Lat)*cos(xyz_Lon) * cos(2*PI*xyz_xi) /Ro
  call AssertEqual(&
       message='xyz_GradLat_wc with Y_1^1 cos(2z)',                  &
       answer = xyz_VLat,                                            &
       check = xyz_GradLat_wc(wc_xyz(xyz_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_Div = - 2 * xyz_Psi/Ro**2
  call AssertEqual(&
       message='wc_Div_xyz_xyz_xyz with Toroiral Y_1^-1 cos(2z)',     &
       answer = xyz_Div,                                              &
       check = xyz_wc(wc_Div_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)),&
       significant_digits = check_digits, ignore_digits = ignore      &
       )

  xyz_VLat = - xyz_Vlat
  call AssertEqual(&
       message='wc_RotRad_xyz_xyz_xyz with Toroiral Y_1^-1 cos(2z)',  &
       answer = xyz_Div,                                              &
       check = xyz_wc(wc_RotRad_xyz_xyz(xyz_VLon,xyz_VLat)),          &
       significant_digits = check_digits, ignore_digits = ignore      &
       )


  ! ------ Y_2^-1 sin(z) Poloidal field -----

  xyz_Psi = cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * sin(PI*xyz_xi) ! Y_2^-1 sin(z)

  xyz_VLon = xyz_GradLon_wc(wc_DRad_ws(ws_xyz(xyz_Psi)))*Ro ! PI*sin(lat)cos(lon))cos(z)
  xyz_VLat = xyz_GradLat_wc(wc_DRad_ws(ws_xyz(xyz_Psi)))*Ro ! PI*cos(2lat)sin(lon)cos(z)
  xyz_Vrad = -xyz_wz(wa_Lapla_wa(wz_xyz(xyz_Psi)))/Ro       ! (2*3/Ro) Psi

  xyz_DData = 0.0D0
  call AssertEqual(&
       message='wc_Div_xyz_xyz_xyz with  Y_2^-1 sin(z)',                &
       answer = xyz_DData,                                              &
       check =  xyz_wc(wc_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad)), &
       significant_digits = check_digits, ignore_digits = ignore        &
       )

  ! ------ Y_2^-1 cos(z) Poloidal field -----

  xyz_Psi = cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * cos(PI*xyz_xi) ! Y_2^-1 sin(z)

  xyz_VLon = xyz_GradLon_ws(ws_DRad_wc(wc_xyz(xyz_Psi)))*Ro ! PI*sin(lat)cos(lon))cos(z)
  xyz_VLat = xyz_GradLat_ws(ws_DRad_wc(wc_xyz(xyz_Psi)))*Ro ! PI*cos(2lat)sin(lon)cos(z)
  xyz_Vrad = -xyz_wz(wa_Lapla_wa(wz_xyz(xyz_Psi)))/Ro       ! (2*3/Ro) Psi

  xyz_DData = 0.0D0
  call AssertEqual(&
       message='ws_Div_xyz_xyz_xyz with  Y_2^-1 sin(z)',                &
       answer = xyz_DData,                                              &
       check =  xyz_ws(ws_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad)), &
       significant_digits = check_digits, ignore_digits = ignore        &
       )

!!$  !========== xyz_RotLon_ws_wc, xyz_RotLat_wc_ws ============
!!$
!!$  ! ------ Y_2^-1 sin(z) Poloidal field -----
!!$
!!$  xyz_Psi = cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * sin(PI*xyz_xi) ! Y_2^-1 sin(z)
!!$
!!$  xyz_VLon = xyz_GradLon_wc(wc_DRad_ws(ws_xyz(xyz_Psi)))*Ro ! PI*sin(lat)cos(lon))cos(z)
!!$  xyz_VLat = xyz_GradLat_wc(wc_DRad_ws(ws_xyz(xyz_Psi)))*Ro ! PI*cos(2lat)sin(lon)cos(z)
!!$  xyz_Vrad = -xyz_wz(wa_Lapla_wa(wz_xyz(xyz_Psi)))/Ro ! (2*3/Ro) Psi
!!$
!!$  xyz_DData = 0.0D0
!!$  call AssertEqual(&
!!$       message='wc_Div_xyz_xyz_xyz with  Y_2^-1 sin(z)',                &
!!$       answer = xyz_DData,                                              &
!!$       check =  xyz_wc(wc_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad)), &
!!$       significant_digits = check_digits, ignore_digits = ignore        &
!!$       )
!!$  
!!$  xyz_DData = (2*3/Ro**2 + PI**2) * cos(2*xyz_Lat) * sin(xyz_Lon) * sin(PI*xyz_xi)
  
!!$  call AssertEqual(&
!!$       message='xyz_RadLon_wt with  Y_2^-1 sin(z)',                   &
!!$       answer = xyz_DData,                                            &
!!$       check =  xyz_RotLon_ws_wc(ws_xyz(xyz_Vrad),wc_xyz(xyz_Vlat)),  &
!!$       significant_digits = check_digits, ignore_digits = ignore      &
!!$       )
!!$
!!$  xyz_DData = -(2*3/Ro**2 + PI**2) * sin(xyz_Lat)*cos(xyz_Lon)*sin(PI*xyz_xi)
!!$  call AssertEqual(&
!!$       message='xyz_RadLat_wt with  Y_2^-1 sin(z)',                   &
!!$       answer = xyz_DData,                                            &
!!$       check =  xyz_RotLat_wc_ws(wc_xyz(xyz_Vlon),ws_xyz(xyz_Vrad)),  &
!!$       significant_digits = check_digits, ignore_digits = ignore      &
!!$       )
!!$

!!$  xyz_VLat = xyz_Rad**(n-1)*cos(2*xyz_Lat)*sin(xyz_Lon)
!!$  call AssertEqual(&
!!$       message='xyz_GradLat_wt with r^2 Y_2^-1',                     &
!!$       answer = xyz_VLat,                                            &
!!$       check = xyz_GradLat_wt(wt_xyz(xyz_Psi)),                      &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
!!$  xyz_Div = - 6* xyz_Psi/xyz_Rad**2
!!$  call AssertEqual(&
!!$       message='wt_Div_xyz_xyz_xyz_wt with Toroiral r2 Y_1^-1',      &
!!$       answer = xyz_Div,                                             &
!!$       check = xyz_wt(wt_Div_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)),&
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
!!$  !==================== wt_DivLon_xyz, wt_DivLat_xyz ====================
!!$  n=2
!!$  write(6,*)'n =',n
!!$
!!$  xyz_VLon  = xyz_rad**n * cos(xyz_Lat)**2*sin(xyz_Lon)
!!$  xyz_DData = xyz_rad**(n-1)*cos(xyz_Lat)*cos(xyz_Lon)
!!$  call AssertEqual(&
!!$       message='wt_DivLon_xyz with r^n cos^2 φ sinλ',              &
!!$       answer = xyz_DData,                                           &
!!$       check = xyz_wt(wt_DivLon_xyz(xyz_VLon)),                      &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
!!$  xyz_VLat   = xyz_rad**n * cos(xyz_Lat)**2*sin(xyz_Lon)
!!$  xyz_DData = -3*xyz_rad**(n-1)*cos(xyz_Lat)*sin(xyz_Lat)*sin(xyz_Lon)
!!$  call AssertEqual(&
!!$       message='wt_DivLat_xyz with r^n cos^2 φ sinλ',              &
!!$       answer = xyz_DData,                                           &
!!$       check = xyz_wt(wt_DivLat_xyz(xyz_VLat)),                      &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
  call MessageNotify('M','wsc_module_svpack_derivative_test1', &
       'wsc_module_svpack derivative function test #1 succeeded!')

end program wsc_module_svpack_derivative1

