From 00e4d345ac3a40e7d1d18df7cf751833cd608612 Mon Sep 17 00:00:00 2001
From: Wojtek Kosior <kwojtus@protonmail.com>
Date: Fri, 21 Jun 2019 13:22:58 +0200
Subject: compozite newton-cotes integration

---
 src/quadratures.f90 | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/src/quadratures.f90 b/src/quadratures.f90
index a3f869f..9be7517 100644
--- a/src/quadratures.f90
+++ b/src/quadratures.f90
@@ -45,9 +45,13 @@ CONTAINS
     real(kind=8), intent(in) :: iend
     procedure(funint) :: fun
     integer(kind=4), intent(in) :: p
-    real(kind=8) :: val
+    
+    real(kind=8) :: val, subinterval_width, qbeg, qend
 
     procedure(quadrature), pointer :: quad
+    integer(kind=8) :: i
+
+    subinterval_width = (iend - ibeg) / subintervals
     
     SELECT CASE (p)
     CASE (:-1)
@@ -63,7 +67,13 @@ CONTAINS
     END SELECT
 
     ! compute integral using quadrature pointed by quad
-    val = quad(ibeg, iend, fun)
+    val = 0
+    
+    DO i = 1, subintervals
+       qend = ibeg + i * subinterval_width
+       qbeg = ibeg + (i - 1) * subinterval_width
+       val = val + quad(qbeg, qend, fun)
+    END DO
   END FUNCTION newton_cotes
 
   FUNCTION rectangle(qbeg, qend, fun) result(val)
-- 
cgit v1.2.3