76 lines
1.7 KiB
Fortran
76 lines
1.7 KiB
Fortran
module Integration
|
|
implicit none
|
|
|
|
contains
|
|
|
|
! function, lower limit, upper limit, steps, method
|
|
function integrate(f, a, b, in, method)
|
|
real :: integrate
|
|
real, intent(in) :: a, b
|
|
integer, optional, intent(in) :: in
|
|
character(len=*), intent(in), optional :: method
|
|
interface
|
|
elemental function f(ra)
|
|
real :: f
|
|
real, intent(in) :: ra
|
|
end function f
|
|
end interface
|
|
|
|
integer :: n, i, m
|
|
real :: h
|
|
real, dimension(:), allocatable :: xpoints
|
|
real, dimension(:), target, allocatable :: fpoints
|
|
real, dimension(:), pointer :: fleft, fmid, fright
|
|
|
|
if ( present(in) ) then
|
|
n = in
|
|
else
|
|
n = 20
|
|
end if
|
|
|
|
if ( present(method) ) then
|
|
select case (method)
|
|
case ('leftrect')
|
|
m = 1
|
|
case ('midrect')
|
|
m = 2
|
|
case ('rightrect')
|
|
m = 3
|
|
case ( 'trapezoid' )
|
|
m = 4
|
|
case default
|
|
m = 0
|
|
end select
|
|
else
|
|
m = 0
|
|
end if
|
|
|
|
h = (b - a) / n
|
|
|
|
allocate(xpoints(0:2*n), fpoints(0:2*n))
|
|
|
|
xpoints = (/ (a + h*i/2, i = 0,2*n) /)
|
|
|
|
fpoints = f(xpoints)
|
|
fleft => fpoints(0 : 2*n-2 : 2)
|
|
fmid => fpoints(1 : 2*n-1 : 2)
|
|
fright => fpoints(2 : 2*n : 2)
|
|
|
|
select case (m)
|
|
case (0) ! simpson
|
|
integrate = h / 6.0 * sum(fleft + fright + 4.0*fmid)
|
|
case (1) ! leftrect
|
|
integrate = h * sum(fleft)
|
|
case (2) ! midrect
|
|
integrate = h * sum(fmid)
|
|
case (3) ! rightrect
|
|
integrate = h * sum(fright)
|
|
case (4) ! trapezoid
|
|
integrate = h * sum(fleft + fright) / 2
|
|
end select
|
|
|
|
deallocate(xpoints, fpoints)
|
|
end function integrate
|
|
|
|
end module Integration
|