diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/quadratures.f90 | 49 |
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 |