!----------------------------------------------------------------------
!     Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wsc_module_svpack テストプログラム
!
!      トロイダルポロイダル関係微分関数のテスト
!  
!履歴  2017/05/20  竹広真一
!
program wsc_module_svpack_derivative2

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

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

  ! 判定誤差設定
  integer, parameter :: check_digits = 8
  integer, parameter :: ignore = -9

  real(8), dimension(0:im-1,jm,0:km)     :: xyz_DData
  real(8), dimension(0:im-1,jm,0:km)     :: xyz_Psi
  real(8), dimension((nm+1)**2,0:lm)     :: wc_Psi
  real(8), dimension((nm+1)**2,lm)       :: ws_Phi

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

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

  real(8), parameter  :: pi=3.1415926535897932385D0
  
  call MessageNotify('M','wsc_module_svpack_derivative_test2', &
       'wsc_module_svpack derivative function test #2')

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

  xyz_xi = (xyz_Rad - ri)/(ro-ri)
  
  !========== wc_KxRGrad_wc, ws_KxRGrad_ws ==========

  ! ----------------- 例 1 --------------------
  xyz_Psi = sin(Pi*xyz_xi) * cos(xyz_lat)*sin(xyz_lon)   ! Y_1^-1 sin(z)

  xyz_DData = sin(Pi*xyz_xi) *  cos(xyz_lat)*cos(xyz_lon)
  call AssertEqual(&
       message='ws_KxRGrad_ws with Y_1^-1 sin(z)',                   &
       answer = xyz_DData,                                           &
       check = xyz_ws(ws_KxRGrad_ws(ws_xyz(xyz_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 例 2 --------------------
  xyz_Psi = cos(xyz_lat)*sin(xyz_lat) * sin(xyz_lon) * cos(PI*xyz_xi)  ! Y_2^1 cos(z)

  xyz_DData = cos(xyz_lat)*sin(xyz_lat) * cos(xyz_lon) * cos(PI*xyz_xi) ! Y_2^1
  call AssertEqual(&
       message='wc_KxRGrad_wc with Y_2^1 cos(z)',                    &
       answer = xyz_DData,                                           &
       check = xyz_wc(wc_KxRGrad_wc(wc_xyz(xyz_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !========== wsc_Potential2Vector, wsc_Potential2Rotation, wc_RadRot, ws_RadRotRot ==========

  !----------------- 剛体回転場 --------------------
  ! 
  wc_Psi = wc_xyz(sin(xyz_Lat))
  ws_Phi = 0.0D0

  call wsc_Potential2Vector(&
       xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi )

  xyz_DData = cos(xyz_Lat)
  call AssertEqual(&
      message='wsc_Potential2Vector (Lon) with rigid rotation field(E-W)', &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2Vector (Lat) with rigid rotation field(E-W)', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2Vector (Rad) with rigid rotation field(E-W)', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_RadRot = xyz_wc(wc_L2_wc(wc_Psi))
  call AssertEqual(&
       message='wc_RadRot_xyz_xyz with with rigid rotation field(E-W)', &
       answer = xyz_RadRot,                                          &
       check = xyz_wc(wc_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_RadRotRot = 0.0D0
  call AssertEqual(&
       message='ws_RadRotRot_xyz_xyz with rigid rotation field(E-W)', &
       answer = xyz_RadRotRot,                                       &
       check = xyz_ws(ws_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wsc_Potential2Rotation(&
       xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi )

  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2Rotation (Lon) with rigid rotation field(E-W)', &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2Rotation (Lat) with rigid rotation field(E-W)', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = xyz_RadRot/Ro
  call AssertEqual(&
      message='wsc_Potential2Rotation (Rad) with rigid rotation field(E-W)', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !----------------- 剛体回転場(南北流) --------------------
  ! 
  wc_Psi = wc_xyz(cos(xyz_Lat) * sin(xyz_Lon))
  ws_Phi = 0.0D0

  call wsc_Potential2Vector(&
       xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi )

  xyz_DData = -sin(xyz_Lat)*sin(xyz_Lon)
  call AssertEqual(&
      message='wsc_Potential2Vector (Lon) with rigid rotation field(N-S)', &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = -cos(xyz_Lon)
  call AssertEqual(&
      message='wsc_Potential2Vector (Lat) with rigid rotation field(N-S)', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2Vector (Rad) with rigid rotation field(N-S)', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_RadRot = xyz_wc(wc_L2_wc(wc_Psi))
  call AssertEqual(&
       message='wc_RadRot_xyz_xyz with with rigid rotation field(N-S)', &
       answer = xyz_RadRot,                                          &
       check = xyz_wc(wc_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_RadRotRot = 0.0D0
  call AssertEqual(&
       message='ws_RadRotRot_xyz_xyz with rigid rotation field(N-S)', &
       answer = xyz_RadRotRot,                                       &
       check = xyz_ws(ws_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wsc_Potential2Rotation(&
       xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi )

  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2Rotation (Lon) with rigid rotation field(N-S)', &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2Rotation (Lat) with rigid rotation field(N-S)', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = xyz_RadRot/Ro
  call AssertEqual(&
      message='wsc_Potential2Rotation (Rad) with rigid rotation field(N-S)', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 渦無し場 --------------------
  wc_Psi = 0.0D0
  ws_Phi = ws_xyz(sin(PI*xyz_xi) * sin(xyz_Lat))

  call wsc_Potential2Vector(&
       xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi )

  xyz_DData = 0
  call AssertEqual(&
       message='wsc_Potential2Vector (Lon) with non-vortical field', &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = PI * cos(xyz_Lat) * cos(PI*xyz_xi) 
  call AssertEqual(&
       message='wsc_Potential2Vector (Lat) with non-vortical field', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 2 * sin(xyz_Lat) * sin(PI*xyz_xi)/Ro
  call AssertEqual(&
      message='wsc_Potential2Vector (Rad) with non-vortical field',  &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_RadRot = 0.0D0
  call AssertEqual(&
       message='wc_RadRot_xyz_xyz with non-vortical field',          &
       answer = xyz_RadRot,                                          &
       check = xyz_wc(wc_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_RadRotRot = - xyz_ws(ws_L2_ws(ws_Lapla_ws((ws_Phi))))
  call AssertEqual(&
       message='ws_RadRotRot_xyz_xyz with non-vortical field',       &
       answer = xyz_RadRotRot,                                       &
       check = xyz_ws(ws_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wsc_Potential2Rotation(&
       xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi )

  xyz_DData = (2/Ro**2 + PI**2) * cos(xyz_Lat) * sin(PI*xyz_xi)
  call AssertEqual(&
      message='wsc_Potential2Rotation (Lon) with non-vortical field', &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2Rotation (Lat) with non-vortical field', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = xyz_RadRot/Ro
  call AssertEqual(&
      message='wsc_Potential2Rotation (Rad) with non-vortical field', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 例 4 --------------------
  ! ポロイダル速度場

  wc_Psi = 0.0D0
  ws_Phi = ws_xyz(sin(2*PI*xyz_xi) * cos(xyz_Lat)*sin(xyz_Lon))

  call wsc_Potential2Vector(&
       xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi )

  xyz_DData = 2*PI * cos(2*PI*xyz_xi) * cos(xyz_Lon)
  call AssertEqual(&
      message='wsc_Potential2Vector (Lon) with poloidal field',   &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = - 2*PI* cos(2*PI*xyz_xi) * sin(xyz_Lat)*sin(xyz_Lon)
  call AssertEqual(&
      message='wsc_Potential2Vector (Lat) with poloidal field', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 2*sin(2*PI*xyz_xi) * cos(xyz_Lat)*sin(xyz_Lon)/Ro
  call AssertEqual(&
      message='wsc_Potential2Vector (Rad) with poloidal field', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_RadRot = 0.0D0
  call AssertEqual(&
       message='wc_RadRot_xyz_xyz with with poloidal field',         &
       answer = xyz_RadRot,                                          &
       check = xyz_wc(wc_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_RadRotRot = - xyz_ws(ws_L2_ws(ws_Lapla_ws((ws_Phi))))
  call AssertEqual(&
       message='ws_RadRotRot_xyz_xyz with with poloidal field',      &
       answer = xyz_RadRotRot,                                       &
       check = xyz_ws(ws_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wsc_Potential2Rotation(&
       xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi )

  xyz_DData = -(2/Ro**2 + 4*PI**2)*sin(2*PI*xyz_xi)*sin(xyz_Lat)*sin(xyz_Lon)
  call AssertEqual(&
      message='wsc_Potential2Rotation (Lon) with poloidal field',    &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = -(2/Ro**2+4*PI**2) * sin(2*PI*xyz_xi) * cos(xyz_Lon)
  call AssertEqual(&
      message='wsc_Potential2Rotation (Lat) with poloidal field', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = xyz_RadRot/Ro
  call AssertEqual(&
      message='wsc_Potential2Rotation (Rad) with poloidal field',    &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','wsc_module_svpack_derivative_test2', &
       'wsc_module_svpack derivative function test #2 succeeded!')

end program wsc_module_svpack_derivative2
