sll_m_trapz_integration.F90

Go to the documentation of this file.

 
module sll_m_trapz_integration
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_working_precision.h"
 
   implicit none
 
   public :: &
      sll_o_trapz_integrate_1d
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
#ifndef DOXYGEN_SHOULD_SKIP_THIS
   abstract interface
 
      function function_1d(x)
         use sll_m_working_precision ! can't pass a header file because the
         ! preprocessor prevents double inclusion.
         ! This is very rare.
         sll_real64             :: function_1d
         sll_real64, intent(in) :: x
      end function function_1d
   end interface
#endif
 
   interface sll_o_trapz_integrate_1d
      module procedure trapz_integral_1d
   end interface
 
contains
 
   function trapz_integral_1d(f, x, n)
      sll_real64                :: trapz_integral_1d
      procedure(function_1d)    :: f
      sll_int32, intent(in)    :: n
      sll_real64, dimension(n)  :: x
      sll_int32                 :: k
      sll_real64                :: ans
 
      ans = 0.0_f64
      do k = 1, n - 1
         ans = ans + 0.5_f64*(f(x(k)) + f(x(k + 1)))*(x(k + 1) - x(k))
      end do
      trapz_integral_1d = ans
 
   end function trapz_integral_1d
 
   function trapz_weights(n, x) result(w)
      sll_int32, intent(in)           :: n
      sll_real64, dimension(n)         :: x
      sll_real64, dimension(n)         :: w
      sll_int32                        :: k
 
      w(1) = 0.5_f64*(x(2) - x(1))
      do k = 2, n - 1
         w(k) = 0.5_f64*(x(k + 1) - x(k - 1))
      end do
      w(n) = 0.5_f64*(x(n) - x(n - 1))
 
   end function trapz_weights
 
end module sll_m_trapz_integration