!----------------------------------------------------------------------
!     Copyright (c) 2024 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wa_zonal_module テストプログラム :: 微分関数のテスト
!
!履歴  2024/02/17  竹広真一
!
program wa_zonal_module_deriv_test

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

  integer, parameter :: im=1, jm=16, nm=10, km=2

  real(8), dimension(0:im-1,1:jm,km)     ::  xya_data          ! 元の関数
  real(8), dimension(0:im-1,1:jm,km)     ::  xya_ddata         ! 微分の正解
  real(8), dimension(0:im-1,1:jm)        ::  xy_mu             ! μ=sinφ

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

  integer :: j

  call MessageNotify('M','wa_zonal_module_deriv_test', &
                         'wa_zonal_module derivative function tests') 

  call wa_Initial( nm, im, jm, km )

  !---- P_1^0 P_2^2 のテスト ----
  xya_data(:,:,1) = sqrt(3.0D0)*sin(xy_Lat)                          ! P_1^0
  xya_data(:,:,2) = sqrt(5.0D0)*(3.0D0*sin(xy_Lat)**2-1.0D0)/2.0D0   ! P_2^0

  xya_ddata(:,:,1) = -2*xya_data(:,:,1)  ! cla_Lapla_cla
  xya_ddata(:,:,2) = -6*xya_data(:,:,2)  ! 

  call AssertEqual(&
    message='P_1^0, P_2^0 Test of wa_Lapla_wa',                 &
    answer = xya_ddata,                                           &
    check = xya_wa(wa_Lapla_wa(wa_xya(xya_data))),            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = -1.0D0/2.0*xya_data(:,:,1)      ! wa_LaplaInv_wa
  xya_ddata(:,:,2) = -1.0D0/6.0*xya_data(:,:,2)

  call AssertEqual(&
    message='P_1^0, P_2^2 Test of wa_LaplaInv_wa',              &
    answer = xya_ddata,                                           &
    check = xya_wa(wa_LaplaInv_wa(wa_xya(xya_data))),         &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata = 0.0D0    ! wa_DLon_wa
  call AssertEqual(&
    message='P_1^0, P_2^0 Test of wa_DLon_wa',                  &
    answer = xya_ddata,                                         &
    check = xya_wa(wa_DLon_wa(wa_xya(xya_data))),               &
    significant_digits = check_digits, ignore_digits = ignore   &
    )

  xya_ddata = 0.0D0    ! xya_GradLon_wa
  call AssertEqual(&
    message='P_1^0, P_2^0 Test of xya_GradLon_wa',                &
    answer = xya_ddata,                                           &
    check = xya_GradLon_wa(wa_xya(xya_data)),                     &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = sqrt(3.0D0)*cos(xy_Lat)         ! xya_GradLat_wa
  xya_ddata(:,:,2) = sqrt(5.0D0)*3.0D0*sin(xy_Lat)*cos(xy_Lat)
  call AssertEqual(&
    message='P_1^0, P_2^0 Test of xya_GradLat_wa',               &
    answer = xya_ddata,                                          &
    check = xya_GradLat_wa(wa_xya(xya_data)),                    &
    significant_digits = check_digits, ignore_digits = ignore    &
    )

  !---- P_1^0 cosφ, P_2^0 cosφ のテスト ----
  xya_data(:,:,1) = sqrt(3.0D0)*sin(xy_Lat)*cos(xy_Lat)        ! P_1 cosφ
  xya_data(:,:,2) =  sqrt(5.0D0)*(3.0D0*sin(xy_Lat)**2-1.0D0)/2.0D0*cos(xy_Lat) ! P_2^0 cosφ

  xya_ddata = 0.0D0
  call AssertEqual(&
    message='P_1^0 cos(phi), P_2^0 cos(phi) Test of wa_DivLon_xya', &
    answer = xya_ddata,                                             &
    check = xya_wa(wa_DivLon_xya(xya_data)),                        &
    significant_digits = check_digits, ignore_digits = ignore       &
    )

  xya_ddata(:,:,1) = sqrt(3.0D0)*(1.0D0-3.0D0*sin(xy_Lat)**2) ! wa_DivLat_xya
  xya_ddata(:,:,2) = sqrt(5.0D0)*sin(xy_Lat)*(3.0D0*cos(2*xy_Lat)+1.0D0) 

  call AssertEqual(&
    message='P_1^0 cos(phi), P_2^0 cos(phi) Test of wa_DivLat_xya', &
    answer = xya_ddata,                                             &
    check = xya_wa(wa_DivLat_xya(xya_data)),                        &
    significant_digits = check_digits, ignore_digits = ignore       &
    )

  !============== 微分計算 (λ,μ座標系用) のテスト ==============
  xy_mu = sin(xy_Lat)

  !----- P_2^0, P_3^0 のテスト -----
  xya_data(:,:,1) = (3.0D0*xy_mu**2-1.0D0)/2.0D0            ! P_2^0
  xya_data(:,:,2) = (5*xy_mu**3-3*xy_mu)/2.0D0              ! P_3^0

  xya_ddata(:,:,1) = (1-xy_mu**2)*3*xy_mu                  ! (1-μ^2)∂/∂μ 
  xya_ddata(:,:,2) = (1-xy_mu**2)*(15*xy_mu**2-3)/2.0D0

  call AssertEqual(&
    message='P_2^0,P_3^0 Test of xya_GradMu_wa',                 &
    answer = xya_ddata,                                          &
    check = xya_GradMu_wa(wa_xya(xya_data)),                     &
    significant_digits = check_digits, ignore_digits = ignore    &
    )

  !----- P_2 (1-μ^2), P_3 (1-μ^2) のテスト -----
  xya_data(:,:,1) = (3.0D0*xy_mu**2-1.0D0)/2.0D0*(1-xy_mu**2)  ! P_2^0 (1-μ^2)
  xya_data(:,:,2) = (5*xy_mu**3-3*xy_mu)/2.0D0*(1-xy_mu**2)    ! P_3 (1-μ^2)

  xya_ddata(:,:,1) = -6*xy_mu**3+4*xy_mu                       ! ∂/∂μ
  xya_ddata(:,:,2) = (-25*xy_mu**4+24*xy_mu**2-3)/2.0D0

  call AssertEqual(&
    message='P_2 (1-mu^2), P_3 (1-mu^2) Test of wa_DivMu_xya',   &
    answer = xya_ddata,                                          &
    check = xya_wa(wa_DivMu_xya(xya_data)),                      &
    significant_digits = check_digits, ignore_digits = ignore    &
    )

  call MessageNotify('M','wa_zonal_module_deriv_test', &
                         'wa_zonal_module derivative function tests succeeded!') 

end program wa_zonal_module_deriv_test