make quadratures a module, sketch the integrate fuction for newton-cotes quadratures

1 files changed, 30 insertions, 15 deletions

diff --git a/src/quadratures.f90 b/src/quadratures.f90
index 9e6bd38..90bba6e 100644
--- a/src/quadratures.f90
+++ b/src/quadratures.f90
@@ -1,7 +1,11 @@
-PROGRAM quadratures
+MODULE quadratures
 
+  IMPLICIT none
+  
+  integer(kind=8) :: subintervals = 100
+  
   INTERFACE
-     FUNCTION integrate(ibeg, iend, myfun, p) result(value)
+     FUNCTION integrate(ibeg, iend, myfun, p) result(val)
        IMPLICIT none
        ! beginning of integration interval
        real(kind=8), intent(in) :: ibeg
@@ -12,7 +16,7 @@ PROGRAM quadratures
        ! polynomial order
        integer(kind=4), intent(in) :: p
        ! result of integration
-       real(kind=8) :: value
+       real(kind=8) :: val
      END FUNCTION integrate
   END INTERFACE
   
@@ -23,19 +27,30 @@ PROGRAM quadratures
        real(kind=8) :: y
      END FUNCTION funint
   END INTERFACE
-
-  procedure(funint), pointer :: ptr
-
-  ptr => my_exp
-
-  write(*,*) my_exp(2.0_8)
   
 CONTAINS
+  
+  FUNCTION newton_cotes(ibeg, iend, fun, p) result(val)
+    IMPLICIT none
+    real(kind=8), intent(in) :: ibeg
+    real(kind=8), intent(in) :: iend
+    procedure(funint) :: fun
+    integer(kind=4), intent(in) :: p
+    real(kind=8) :: val
+
+    SELECT CASE (p)
+    CASE (:-1)
+       STOP "negative interpolationg polynomial order passed"
+    CASE (0)
+       val = 0 ! rectangle(ibeg, iend, fun)
+    CASE (1)
+       val = 0 ! trapeze(ibeg, iend, fun)
+    CASE (2)
+       val = 0 ! simpson_1_third(ibeg, iend, fun)
+    CASE default
+       STOP "Newton-Cotes quadratures only implemented for order < 3"
+    END SELECT
+  END FUNCTION newton_cotes
 
-  FUNCTION my_exp(x) result(y)
-    real(kind=8), intent(in) :: x
-    real(kind=8) :: y
-    y = exp(x)
-  END FUNCTION my_exp
   
-END PROGRAM quadratures
+END MODULE quadratures