Line data    Source code

       1             : MODULE m_trapz
       2             :    !General Purpose trapezian method integration
       3             :    !Used in green's function calculations because the
       4             :    !integrands are very spiky
       5             : 
       6             :    IMPLICIT NONE
       7             : 
       8             :    INTERFACE trapz
       9             :       PROCEDURE :: trapzr, trapzc
      10             :    END INTERFACE
      11             : 
      12             :    CONTAINS
      13             : 
      14         788 :    PURE REAL FUNCTION trapzr(y,h,n)
      15             : 
      16             :       REAL,          INTENT(IN)     :: y(:)
      17             : 
      18             :       INTEGER,       INTENT(IN)     :: n
      19             :       REAL,          INTENT(IN)     :: h
      20             : 
      21             :       INTEGER i
      22             : 
      23         788 :       trapzr = y(1)
      24     4784064 :       DO i = 2, n-1
      25     4784064 :          trapzr = trapzr + 2*y(i)
      26             :       ENDDO
      27         788 :       trapzr = trapzr + y(n)
      28             : 
      29         788 :       trapzr = trapzr*h/2.0
      30             : 
      31         788 :    END FUNCTION trapzr
      32             : 
      33           0 :    PURE COMPLEX FUNCTION trapzc(y,h,n)
      34             : 
      35             :       COMPLEX,       INTENT(IN)     :: y(:)
      36             : 
      37             :       INTEGER,       INTENT(IN)     :: n
      38             :       REAL,          INTENT(IN)     :: h
      39             : 
      40             :       INTEGER i
      41             : 
      42           0 :       trapzc = y(1)
      43           0 :       DO i = 2, n-1
      44           0 :          trapzc = trapzc + 2*y(i)
      45             :       ENDDO
      46           0 :       trapzc = trapzc + y(n)
      47             : 
      48           0 :       trapzc = trapzc*h/2.0
      49             : 
      50           0 :    END FUNCTION trapzc
      51             : 
      52             : 
      53             : END MODULE m_trapz