!----------------------------------------------------------------------
! COPYRIGHT (c) 2007-2010 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!
!表題  au_module テストプログラム
!
!履歴  2007/12/28  竹広真一
!      2008/01/05  竹広真一  ag_Dr2_au のテスト追加
!      2010/12/21  佐々木洋平 dc_test を使用するように修正
!
program au_test_base2d

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use au_module
  implicit none
  integer, parameter :: im=32, km=32, nm=4
  real(8), parameter :: ra=0.5D0
  
  real(8), dimension(0:nm,0:im) :: ag_y
  real(8), dimension(0:nm,0:im) :: ag_y_deriv
  real(8), dimension(0:nm,0:im) :: ag_y_deriv2

  real(8), dimension(0:nm,0:km) :: au_y_sol
  integer, dimension(0:nm)      :: nd=(/0,1,2,3,4/)

  real(8), dimension(0:im) :: g_x

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


  call MessageNotify('M','au_test_base2d','au_module 2-dim function tests')

  call au_initial(im,km,ra,nd)

  g_x = 2.0D0*g_R**2/ra**2-1.0D0
  
  ag_y(0,:) = 1.0d0                                        ! T_0
  ag_y(1,:) = g_R*g_x                                      ! T_1
  ag_y(2,:) = g_R**2*(2.0D0*g_x**2 - 1.0D0)                ! T_2
  ag_y(3,:) = g_R**3*(4.0D0*g_x**3 - 3.0D0*g_x)            ! T_3
  ag_y(4,:) = g_R**4*(8.0D0*g_x**4 - 8.0D0*g_x**2 + 1.0D0) ! T_4

  ag_y_deriv(0,:) = 0.0D0                
  ag_y_deriv(1,:) = g_x + 4*g_R**2/ra**2 
  ag_y_deriv(2,:) = 2.0D0*g_R*(2.0D0*g_x**2 - 1)+ g_R**2*4*g_x*4*g_R/ra**2
  ag_y_deriv(3,:) = 3.0D0*g_R**2.0D0*(4.0D0*g_x**3 - 3.0D0*g_x)            &
                  +(12.0D0*g_x**2 - 3.0D0)*4.0D0*g_R**4/ra**2
  ag_y_deriv(4,:) = 4.0D0*g_R**3.0D0*(8.0D0*g_x**4 - 8.0D0*g_x**2 + 1.0D0) &
                  +g_R**4 * (32.0D0*g_x**3 - 16.0D0*g_x)*4.0D0*g_R/ra**2

  ag_y_deriv2(0,:) = 0.0D0                
  ag_y_deriv2(1,:) = 12.0D0*g_R/ra**2
  ag_y_deriv2(2,:) = 2.0D0*(2.0D0*g_x**2 - 1.0D0)                        &
                   + 80.0D0*g_R**2*g_x/ra**2 +64.0D0*g_R**4/ra**4
  ag_y_deriv2(3,:) = 6.0D0*g_R*(4.0D0*g_x**3 - 3.0D0*g_x)                &
                   + 28.0D0*g_R**3/ra**2*(12.0D0*g_x**2-3.0D0)           &
                   + (24.0D0*g_x)*16.0D0*g_R**5/ra**4
  ag_y_deriv2(4,:) = 12.0D0*g_R**2*(8.0D0*g_x**4 - 8.0D0*g_x**2 + 1.0D0) &
                   + 36.0D0*g_R**4/ra**2*(32.0D0*g_x**3 - 16.0D0*g_x)    &
                   + 16.0D0*g_R**6/ra**4*(96.0D0*g_x**2 - 16.0D0)

  au_y_sol(0,:) = 0.0D0; au_y_sol(0,0) =1.0D0          ! T_0
  au_y_sol(1,:) = 0.0D0; au_y_sol(1,1) =1.0D0          ! T_1
  au_y_sol(2,:) = 0.0D0; au_y_sol(2,2) =1.0D0          ! T_2
  au_y_sol(3,:) = 0.0D0; au_y_sol(3,3) =1.0D0          ! T_3
  au_y_sol(4,:) = 0.0D0; au_y_sol(4,4) =1.0D0          ! T_4


  call MessageNotify('M','au_test_base2d', 'Test of ag_y(n) = T_n(x), n=0...4')
  
  call AssertEqual(&
    message='au_ag(ag_y)',                                        &
    answer = au_y_sol,                                            &
    check = au_ag(ag_y),                                          &
    significant_digits = check_digits, ignore_digits = ignore     &
    )


  call AssertEqual(&
    message='ag_au(au_ag(ag_y))',                                 &
    answer = ag_y,                                                &
    check = ag_au(au_ag(ag_y)),                                   &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call AssertEqual(&
    message='ag_Dr_au(au_ag(ag_y))',                              &
    answer = ag_y_deriv,                                          &
    check = ag_Dr_au(au_ag(ag_y)),                                &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call AssertEqual(&
    message='ag_Dr2_au(au_ag(ag_y))',                             &
    answer = ag_y_deriv2,                                         &
    check = ag_Dr2_au(au_ag(ag_y)),                               &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call MessageNotify('M','au_test_base2d','au_module 2-dim test suceeded!')

end program au_test_base2d