From 7d28fe599e1a5498a7e86dc04fa207cab11559ac Mon Sep 17 00:00:00 2001
From: Wojtek Kosior <kwojtus@protonmail.com>
Date: Fri, 21 Jun 2019 13:08:31 +0200
Subject: newton-cotes integration of entire interval at once

---
 src/quadratures.f90 | 49 +++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 45 insertions(+), 4 deletions(-)

(limited to 'src')

diff --git a/src/quadratures.f90 b/src/quadratures.f90
index 90bba6e..a3f869f 100644
--- a/src/quadratures.f90
+++ b/src/quadratures.f90
@@ -20,6 +20,15 @@ MODULE quadratures
      END FUNCTION integrate
   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
@@ -27,7 +36,7 @@ MODULE quadratures
        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
-- 
cgit v1.2.3