1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
|
MODULE quadratures
IMPLICIT none
integer(kind=8) :: subintervals = 100
INTERFACE
FUNCTION integrate(ibeg, iend, myfun, p) result(val)
IMPLICIT none
! beginning of integration interval
real(kind=8), intent(in) :: ibeg
! ending of integration interval
real(kind=8), intent(in) :: iend
! function to be integrated
procedure(funint) :: myfun
! polynomial order
integer(kind=4), intent(in) :: p
! result of integration
real(kind=8) :: val
END FUNCTION integrate
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)
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
END MODULE quadratures
|
