diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/quadratures.f90 | 72 |
1 files changed, 72 insertions, 0 deletions
diff --git a/src/quadratures.f90 b/src/quadratures.f90 index 43e6ac6..cc279ed 100644 --- a/src/quadratures.f90 +++ b/src/quadratures.f90 @@ -65,6 +65,33 @@ CONTAINS val = integrate_generic(ibeg, iend, fun, quad) END FUNCTION newton_cotes + FUNCTION gauss(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 + + procedure(quadrature), pointer :: quad + + SELECT CASE (p) + CASE (:0) + STOP "non-positive Legendre polynomial order passed" + CASE (1) + quad => gauss_n1 + CASE (2) + quad => gauss_n2 + CASE (3) + quad => gauss_n3 + CASE default + STOP "Gauss quadratures only implemented with roots of" & + // " Legendgre polynomial of order <= 3" + END SELECT + + val = integrate_generic(ibeg, iend, fun, quad) + END FUNCTION gauss + FUNCTION integrate_generic(ibeg, iend, fun, quad) result(val) IMPLICIT none real(kind=8), intent(in) :: ibeg @@ -140,4 +167,49 @@ CONTAINS (fun(qbeg) + 4 * fun ((qbeg + qend) * 0.5) + fun(qend)) END FUNCTION simpson_1_third + FUNCTION gauss_generic(mid, halfwidth, fun, points, weights) & + result(val) + IMPLICIT none + real(kind=8), intent(in) :: mid, halfwidth, points(:), weights(:) + procedure(funint) :: fun + real(kind=8) :: val + + integer(kind=4) :: i + + val = halfwidth * sum(weights * & + [(fun(points(i) * halfwidth + mid), i = 1, size(points))]) + END FUNCTION gauss_generic + + FUNCTION gauss_n1(qbeg, qend, fun) result(val) + IMPLICIT none + real(kind=8), intent(in) :: qbeg, qend + procedure(funint) :: fun + real(kind=8) :: val, weights(1) = [2], points(1) = [0] + + val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, & + fun, points, weights) + END FUNCTION gauss_n1 + + FUNCTION gauss_n2(qbeg, qend, fun) result(val) + IMPLICIT none + real(kind=8), intent(in) :: qbeg, qend + procedure(funint) :: fun + real(kind=8) :: val, weights(2) = [1, 1], & + points(2) = [1 / sqrt(3.0), -1 / sqrt(3.0)] + + val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, & + fun, points, weights) + END FUNCTION gauss_n2 + + FUNCTION gauss_n3(qbeg, qend, fun) result(val) + IMPLICIT none + real(kind=8), intent(in) :: qbeg, qend + procedure(funint) :: fun + real(kind=8) :: val, weights(3) = [8 / 9.0, 5 / 9.0, 5 / 9.0],& + points(3) = [0.0, sqrt(3 / 5.0), -sqrt(3 / 5.0)] + + val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, & + fun, points, weights) + END FUNCTION gauss_n3 + END MODULE quadratures |