!----------------------------------------------------------------------
!     Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wa_wave_module_svpack テストプログラム :: 微分関数のテスト
!
!履歴  2017/03/09  竹広真一
!
program wa_wave_module_svpack_deriv_test

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

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

  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_wave_module_svpack_deriv_test', &
                         'wa_wave_module_svpack derivative function tests') 

  call wa_Initial( nm, m, im, jm, km )

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

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

  call AssertEqual(&
    message='P_1^1, P_2^1 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^1, 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     &
    )

  do j=1,jm
     xya_ddata(0,j,1) = -m*sqrt(3.0D0/2.0D0)*cos(y_Lat(j))    ! wa_DLon_wa
     xya_ddata(1,j,1) =  m*sqrt(3.0D0/2.0D0)*cos(y_Lat(j))
     xya_ddata(0,j,2) = -m*sqrt(5.0D0/6.0D0)*(3*sin(y_Lat(j))*cos(y_Lat(j)))
     xya_ddata(1,j,2) =  m*sqrt(5.0D0/6.0D0)*(3*sin(y_Lat(j))*cos(y_Lat(j)))
  enddo
  call AssertEqual(&
    message='P_1^1, P_2^1 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   &
    )

  do j=1,jm
     xya_ddata(0,j,1) = -m*sqrt(3.0D0/2.0D0)
     xya_ddata(1,j,1) =  m*sqrt(3.0D0/2.0D0)
     xya_ddata(0,j,2) = -m*sqrt(5.0D0/6.0D0)*(3*sin(y_Lat(j)))
     xya_ddata(1,j,2) =  m*sqrt(5.0D0/6.0D0)*(3*sin(y_Lat(j)))
  enddo
  call AssertEqual(&
    message='P_1^1, P_2^1 Test of xya_GradLon_cla',               &
    answer = xya_ddata,                                           &
    check = xya_GradLon_wa(wa_xya(xya_data)),                     &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = -sqrt(3.0D0/2.D0)*sin(xy_Lat)         ! xya_GradLat_wa
  xya_ddata(:,:,2) = sqrt(5.0D0/6.0D0)*3.0*cos(2*xy_Lat)
  call AssertEqual(&
    message='P_1^1, P_2^2 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 cosφ, P_2 cosφ のテスト ----
  xya_data(:,:,1) = sqrt(3.0D0/2.0D0)*cos(xy_Lat)*cos(xy_Lat)        ! P_1 cosφ
  xya_data(:,:,2) = sqrt(5.0D0/6.0D0)*(3*sin(xy_Lat)*cos(xy_Lat))*cos(xy_Lat) ! P_2^1 cosφ

  do j=1,jm
     xya_ddata(0,j,1) = -m*sqrt(3.0D0/2.0D0)*cos(y_Lat(j))
     xya_ddata(1,j,1) =  m*sqrt(3.0D0/2.0D0)*cos(y_Lat(j))
     xya_ddata(0,j,2) = -m*sqrt(5.0D0/6.0D0)*(3*sin(y_Lat(j))*cos(y_Lat(j)))
     xya_ddata(1,j,2) =  m*sqrt(5.0D0/6.0D0)*(3*sin(y_Lat(j))*cos(y_Lat(j)))
  enddo
  call AssertEqual(&
    message='P_1^1 cos(phi), P_2^2 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) = -3*sqrt(3.0D0/2.0D0)*cos(xy_Lat)*sin(xy_Lat) ! wa_DivLat_xya
  xya_ddata(:,:,2) = sqrt(5.0D0/6.0D0) &
       *(3.0*(cos(xy_Lat)**3-3*sin(xy_Lat)**2*cos(xy_Lat)))    ! wa_DivLat_xya

  call AssertEqual(&
    message='P_1^1 cos(phi), P_2^2 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^1, P_3^1 のテスト -----
  xya_data(:,:,1) = xy_mu*sqrt(1-xy_mu**2)                  ! P_2^1
  xya_data(:,:,2) = (5*xy_mu**2-1)*sqrt(1-xy_mu**2)         ! P_3^1

  xya_ddata(:,:,1) = sqrt(1-xy_mu**2)*(1-2*xy_mu**2)        ! (1-μ^2)∂/∂μ 
  xya_ddata(:,:,2) = xy_mu*(11-15*xy_mu**2)*sqrt(1-xy_mu**2)

  call AssertEqual(&
    message='P_2^1,P_3^1 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_1 (1-μ^2) のテスト -----
  xya_data(:,:,1) = xy_mu*sqrt(1-xy_mu**2)*(1-xy_mu**2)           ! P_2^1 (1-μ^2)
  xya_data(:,:,2) = (5*xy_mu**2-1)*sqrt(1-xy_mu**2)*(1-xy_mu**2)  ! P_3 (1-μ^2)

  xya_ddata(:,:,1) = (1-4*xy_mu**2)*sqrt(1-xy_mu**2)              ! ∂/∂μ
  xya_ddata(:,:,2) = xy_mu*(13-25*xy_mu**2)*sqrt(1-xy_mu**2)


  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_wave_module_svpack_deriv_test', &
                         'wa_wave_module_svpack derivative function tests succeeded!') 

end program wa_wave_module_svpack_deriv_test