newton-cotes integration of entire interval at once

1 files changed, 45 insertions, 4 deletions

diff --git a/src/quadratures.f90 b/src/quadratures.f90
index 90bba6e..a3f869f 100644
--- a/src/quadratures.f90
+++ b/src/quadratures.f90
@@ -21,13 +21,22 @@ MODULE quadratures
   END INTERFACE
   
   INTERFACE
+     FUNCTION quadrature(qbeg, qend, fun) result(val)
+       IMPLICIT none
+       real(kind=8), intent(in) :: qbeg, qend
+       procedure(funint) :: fun
+       real(kind=8) :: val
+     END FUNCTION quadrature
+  END INTERFACE
+  
+  INTERFACE
      FUNCTION funint(x) result(y)
        IMPLICIT none
        real(kind=8), intent(in) :: x
        real(kind=8) :: y
      END FUNCTION funint
   END INTERFACE
-  
+
 CONTAINS
   
   FUNCTION newton_cotes(ibeg, iend, fun, p) result(val)
@@ -38,19 +47,51 @@ CONTAINS
     integer(kind=4), intent(in) :: p
     real(kind=8) :: val
 
+    procedure(quadrature), pointer :: quad
+    
     SELECT CASE (p)
     CASE (:-1)
        STOP "negative interpolationg polynomial order passed"
     CASE (0)
-       val = 0 ! rectangle(ibeg, iend, fun)
+       quad => rectangle
     CASE (1)
-       val = 0 ! trapeze(ibeg, iend, fun)
+       quad => trapeze
     CASE (2)
-       val = 0 ! simpson_1_third(ibeg, iend, fun)
+       quad => simpson_1_third
     CASE default
        STOP "Newton-Cotes quadratures only implemented for order < 3"
     END SELECT
+
+    ! compute integral using quadrature pointed by quad
+    val = quad(ibeg, iend, fun)
   END FUNCTION newton_cotes
 
+  FUNCTION rectangle(qbeg, qend, fun) result(val)
+    IMPLICIT none
+    real(kind=8), intent(in) :: qbeg, qend
+    procedure(funint) :: fun
+    real(kind=8) :: val
+
+    val = (qend - qbeg) * fun((qend + qbeg) * 0.5)
+  END FUNCTION rectangle
+  
+  FUNCTION trapeze(qbeg, qend, fun) result(val)
+    IMPLICIT none
+    real(kind=8), intent(in) :: qbeg, qend
+    procedure(funint) :: fun
+    real(kind=8) :: val
+
+    val = (qend - qbeg) * 0.5 * (fun(qbeg) + fun(qend))
+  END FUNCTION trapeze
   
+  FUNCTION simpson_1_third(qbeg, qend, fun) result(val)
+    IMPLICIT none
+    real(kind=8), intent(in) :: qbeg, qend
+    procedure(funint) :: fun
+    real(kind=8) :: val
+
+    val = (qend - qbeg) * (1/6.0) *                               &
+         (fun(qbeg) + 4 * fun ((qbeg + qend) * 0.5) + fun(qend))
+  END FUNCTION simpson_1_third
+
 END MODULE quadratures