Line data    Source code

       1              : module test_integrate
       2              : 
       3              :    use math_lib
       4              :    use num_def
       5              :    use num_lib
       6              :    use const_def
       7              : 
       8              :    implicit none
       9              : 
      10              : contains
      11              : 
      12            1 :    subroutine run_test_integrate()
      13              : 
      14            1 :       call test_basic()
      15            1 :       call test_sine()
      16            1 :       call test_exp()
      17            1 :       call test_box()
      18              : 
      19            1 :    end subroutine run_test_integrate
      20              : 
      21            1 :    subroutine test_basic
      22              :       real(dp), parameter :: xlow = 0, xhigh = 1
      23              :       real(dp), parameter :: expected = 0.5d0
      24              :       real(dp) :: res
      25              :       integer :: ierr
      26              : 
      27            1 :       res = integrate(linear, xlow, xhigh, [1d0], 1d-3, 1d-3, 10, ierr)
      28              : 
      29            1 :       call check_result('linear', expected, res, ierr)
      30              : 
      31              :    contains
      32              : 
      33           12 :       real(dp) function linear(x, args, ierr)
      34              :          real(dp), intent(in) :: x
      35              :          real(dp), intent(in) :: args(:)
      36              :          integer, intent(inout) :: ierr
      37              : 
      38           12 :          ierr = 0
      39           12 :          linear = x
      40              : 
      41           12 :       end function linear
      42              : 
      43              :    end subroutine test_basic
      44              : 
      45            1 :    subroutine test_sine
      46              :       real(dp) :: res
      47              :       integer :: ierr
      48              : 
      49            1 :       res = integrate(sine, 0d0, pi, [1d0], 1d-8, 1d-8, 10, ierr)
      50              : 
      51            1 :       call check_result('sine', 2d0, res, ierr)
      52              : 
      53            1 :       res = integrate(sine, 0d0, 2*pi, [1d0], 1d-8, 1d-8, 10, ierr)
      54              : 
      55            1 :       call check_result('sine', 0d0, res, ierr)
      56              : 
      57              :    contains
      58              : 
      59          720 :       real(dp) function sine(x, args, ierr)
      60              :          real(dp), intent(in) :: x
      61              :          real(dp), intent(in) :: args(:)
      62              :          integer, intent(inout) :: ierr
      63              : 
      64          720 :          ierr = 0
      65          720 :          sine = sin(x)
      66              : 
      67          720 :       end function sine
      68              : 
      69              :    end subroutine test_sine
      70              : 
      71            1 :    subroutine test_exp
      72              :       real(dp) :: res
      73              :       integer :: ierr
      74              : 
      75            1 :       res = integrate(iexp, 0d0, 2d0, [1d0], 1d-8, 1d-8, 50, ierr)
      76              : 
      77            1 :       call check_result('exp', exp(2d0) - 1d0, res, ierr)
      78              : 
      79            1 :       res = integrate(iexp, 0d0, 10d0, [1d0], 1d-8, 1d-8, 50, ierr)
      80              : 
      81            1 :       call check_result('exp', exp(10d0) - 1d0, res, ierr)
      82              : 
      83              :    contains
      84              : 
      85         3648 :       real(dp) function iexp(x, args, ierr)
      86              :          real(dp), intent(in) :: x
      87              :          real(dp), intent(in) :: args(:)
      88              :          integer, intent(inout) :: ierr
      89              : 
      90         3648 :          ierr = 0
      91         3648 :          iexp = exp(x)
      92              : 
      93         3648 :       end function iexp
      94              : 
      95              :    end subroutine test_exp
      96              : 
      97            1 :    subroutine test_box
      98              :       real(dp) :: res
      99              :       integer :: ierr
     100              : 
     101            1 :       res = integrate(box, 0d0, 2d0, [1d0], 1d-8, 1d-8, 50, ierr)
     102              : 
     103            1 :       call check_result('box', 1d0, res, ierr)
     104              : 
     105            1 :       res = integrate(box, 0.99d0, 1.5d0, [1d0], 1d-8, 1d-8, 50, ierr)
     106              : 
     107            1 :       call check_result('box', 0.5d0, res, ierr)
     108              : 
     109              :    contains
     110              : 
     111         1152 :       real(dp) function box(x, args, ierr)
     112              :          real(dp), intent(in) :: x
     113              :          real(dp), intent(in) :: args(:)
     114              :          integer, intent(inout) :: ierr
     115              : 
     116         1152 :          ierr = 0
     117              : 
     118         1152 :          if (x < 1) then
     119              :             box = 0d0
     120          487 :          else if (x >= 1d0 .and. x <= 2d0) then
     121              :             box = 1d0
     122              :          else
     123          665 :             box = 0d0
     124              :          end if
     125              : 
     126         1152 :       end function box
     127              : 
     128              :    end subroutine test_box
     129              : 
     130            7 :    subroutine check_result(name, tgt, val, ierr)
     131              :       character(len=*), intent(in) :: name
     132              :       real(dp), intent(in) :: tgt, val
     133              :       integer, intent(in) :: ierr
     134              : 
     135              :       write (*, '(a40, 1pd26.16, a7, 1pd26.16, i4)') &
     136            7 :          'integrate '//trim(name)//' expected', &
     137           14 :          tgt, 'got', val, ierr
     138              : 
     139            7 :    end subroutine check_result
     140              : 
     141              : end module test_integrate