aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorWojtek Kosior <kwojtus@protonmail.com>2019-06-21 15:33:18 +0200
committerWojtek Kosior <kwojtus@protonmail.com>2019-06-21 15:33:18 +0200
commit9ca933fd7b44d8913d13a94b54dd7fd1635a717a (patch)
tree272621c32629e64e5105d0d3bcb6e917bb8012d0 /src
parente5a89cb63036f343badf1ebbb4ed07b06b4eed0a (diff)
downloadfortran-assignment3-9ca933fd7b44d8913d13a94b54dd7fd1635a717a.tar.gz
fortran-assignment3-9ca933fd7b44d8913d13a94b54dd7fd1635a717a.zip
add gauss quadratures, generic way (untestedgit add quadratures.f90)
Diffstat (limited to 'src')
-rw-r--r--src/quadratures.f9072
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