diff options
| -rw-r--r-- | src/quadratures.f90 | 45 |
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 |