Performs 1D integration using Simpsons rule
###References * Simpson's rule
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=Float64), | intent(in), | dimension(:) | :: | f | Array of equally spaced values |
|
real(kind=Float64), | intent(in) | :: | dx | Spacing between x values |
function simpsons_rule(f, dx) result(I)
!+ Performs 1D integration using Simpsons rule
!+
!+ ###References
!+* [Simpson's rule](http://mathworld.wolfram.com/SimpsonsRule.html)
real(Float64), dimension(:), intent(in) :: f
!+ Array of equally spaced \(f(x)\) values
real(Float64), intent(in) :: dx
!+ Spacing between x values
real(Float64) :: I
integer :: s, ii
s = size(f)
I = 0.d0
if(mod(s,2).eq.1) then
write(*,'(a)') "Length of array must be even"
return
endif
I = f(1)
do ii=2,s-1
if(mod(ii,2).eq.1) then
I = I + 4.0*f(ii)
else
I = I + 2.0*f(ii)
endif
enddo
I = I + f(s)
I = (dx/3.0)*I
end function simpsons_rule