#:include "common.fypp"
submodule (stdlib_quadrature) stdlib_quadrature_simps
use stdlib_error, only: check
implicit none
! internal use only
interface simps38
#:for k1, t1 in REAL_KINDS_TYPES
module procedure simps38_dx_${k1}$
module procedure simps38_x_${k1}$
#:endfor
end interface simps38
! internal use only
interface simps38_weights
#:for k1, t1 in REAL_KINDS_TYPES
module procedure simps38_weights_${k1}$
#:endfor
end interface simps38_weights
contains
#:for k1, t1 in REAL_KINDS_TYPES
pure recursive module function simps_dx_${k1}$(y, dx, even) result(integral)
${t1}$, dimension(:), intent(in) :: y
${t1}$, intent(in) :: dx
integer, intent(in), optional :: even
${t1}$ :: integral
integer :: n
n = size(y)
select case (n)
case (0:1)
integral = 0.0_${k1}$
case (2)
integral = 0.5_${k1}$*dx*(y(1) + y(2))
case (3)
integral = dx/3.0_${k1}$*(y(1) + 4*y(2) + y(3))
case (4)
integral = simps38(y, dx)
! case (5) not needed; handled by default
case (6) ! needs special handling because of averaged 3/8's rule case
if (present(even)) then
if (even < 0) then
! 3/8 rule on left
integral = simps38(y(1:4), dx) + simps(y(4:6), dx)
return
else if (even > 0) then
! 3/8 rule on right
integral = simps(y(1:3), dx) + simps38(y(3:6), dx)
return
else
! fall through
end if
end if
! either `even` not present or is zero
! equivalent to averaging left and right
integral = dx/48.0_${k1}$ * (17*(y(1) + y(6)) + 59*(y(2) + y(5)) + 44*(y(3) + y(4)))
case default
if (mod(n, 2) == 1) then
integral = dx/3.0_${k1}$*(y(1) + 4*sum(y(2:n-1:2)) + 2*sum(y(3:n-2:2)) + y(n))
else
if (present(even)) then
if (even < 0) then
! 3/8th rule on left
integral = simps38(y(1:4), dx) + simps(y(4:n), dx)
return
else if (even > 0) then
! 3/8 rule on right
integral = simps(y(1:n-3), dx) + simps38(y(n-3:n), dx)
return
else
! fall through
end if
end if
! either `even` not present or is zero
! equivalent to averaging left and right
integral = dx/48.0_${k1}$ * (17*(y(1) + y(n)) + 59*(y(2) + y(n-1)) &
+ 43*(y(3) + y(n-2)) + 49*(y(4) + y(n-3)) + 48*sum(y(5:n-4)))
end if
end select
end function simps_dx_${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
recursive module function simps_x_${k1}$(y, x, even) result(integral)
${t1}$, dimension(:), intent(in) :: y
${t1}$, dimension(:), intent(in) :: x
integer, intent(in), optional :: even
${t1}$ :: integral
integer :: i
integer :: n
${t1}$ :: h1, h2
${t1}$ :: a, b, c
n = size(y)
call check(size(x) == n, "simps: Arguments `x` and `y` must be the same size.")
select case (n)
case (0:1)
integral = 0.0_${k1}$
case (2)
integral = 0.5_${k1}$*(x(2) - x(1))*(y(1) + y(2))
case (3)
h1 = x(2) - x(1)
h2 = x(3) - x(2)
a = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
b = (h1+h2)**3/(6*h1*h2)
c = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
integral = a*y(1) + b*y(2) + c*y(3)
case (4)
integral = simps38(y, x)
! case (6) unneeded; handled by default
case default
if (mod(n, 2) == 1) then
integral = 0.0_${k1}$
do i = 1, n-2, 2
h1 = x(i+1) - x(i)
h2 = x(i+2) - x(i+1)
a = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
b = (h1+h2)**3/(6*h1*h2)
c = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
integral = integral + a*y(i) + b*y(i+1) + c*y(i+2)
end do
else
if (present(even)) then
if (even < 0) then
! 3/8 rule on left
integral = simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n))
return
else if (even > 0) then
! 3/8 rule on right
integral = simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n))
return
else
! fall through
end if
end if
! either `even` not present or is zero
integral = 0.5_${k1}$ * ( simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n)) &
+ simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n)) )
end if
end select
end function simps_x_${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
pure recursive module function simps_weights_${k1}$(x, even) result(w)
${t1}$, dimension(:), intent(in) :: x
integer, intent(in), optional :: even
${t1}$, dimension(size(x)) :: w
integer :: i, n
${t1}$ :: h1, h2
n = size(x)
select case (n)
case (0)
! no action needed
case (1)
w(1) = 0.0_${k1}$
case (2)
w = 0.5_${k1}$*(x(2) - x(1))
case (3)
h1 = x(2) - x(1)
h2 = x(3) - x(2)
w(1) = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
w(2) = (h1+h2)**3/(6*h1*h2)
w(3) = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
case (4)
w = simps38_weights(x)
case default
if (mod(n, 2) == 1) then
w = 0.0_${k1}$
do i = 1, n-2, 2
h1 = x(i+1) - x(i)
h2 = x(i+2) - x(i+1)
w(i) = w(i) + (2*h1**2 + h1*h2 - h2**2)/(6*h1)
w(i+1) = w(i+1) + (h1+h2)**3/(6*h1*h2)
w(i+2) = w(i+2) + (2*h2**2 + h1*h2 - h1**2)/(6*h2)
end do
else
if (present(even)) then
if (even < 0) then
! 3/8 rule on left
w = 0.0_${k1}$
w(1:4) = simps38_weights(x(1:4))
w(4:n) = w(4:n) + simps_weights(x(4:n)) ! position 4 needs both rules
return
else if (even > 0) then
! 3/8 rule on right
w = 0.0_${k1}$
w(1:n-3) = simps_weights(x(1:n-3))
w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n)) ! position n-3 needs both rules
return
else
! fall through
end if
end if
! either `even` not present or is zero
w = 0.0_${k1}$
! 3/8 rule on left
w(1:4) = simps38_weights(x(1:4))
w(4:n) = w(4:n) + simps_weights(x(4:n))
! 3/8 rule on right
w(1:n-3) = w(1:n-3) + simps_weights(x(1:n-3))
w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n))
! average
w = 0.5_${k1}$ * w
end if
end select
end function simps_weights_${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
pure function simps38_dx_${k1}$(y, dx) result (integral)
${t1}$, dimension(4), intent(in) :: y
${t1}$, intent(in) :: dx
${t1}$ :: integral
integral = 3.0_${k1}$*dx/8.0_${k1}$ * (y(1) + y(4) + 3*(y(2) + y(3)))
end function simps38_dx_${k1}$
#:endfor
#: for k1, t1 in REAL_KINDS_TYPES
pure function simps38_x_${k1}$(y, x) result(integral)
${t1}$, dimension(4), intent(in) :: y
${t1}$, dimension(4), intent(in) :: x
${t1}$ :: integral
${t1}$ :: h1, h2, h3
${t1}$ :: a, b, c, d
h1 = x(2) - x(1)
h2 = x(3) - x(2)
h3 = x(4) - x(3)
a = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2))
b = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3))
c = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2))
d = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3))
integral = a*y(1) + b*y(2) + c*y(3) + d*y(4)
end function simps38_x_${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
pure function simps38_weights_${k1}$(x) result(w)
${t1}$, intent(in) :: x(4)
${t1}$ :: w(size(x))
${t1}$ :: h1, h2, h3
h1 = x(2) - x(1)
h2 = x(3) - x(2)
h3 = x(4) - x(3)
w(1) = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2))
w(2) = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3))
w(3) = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2))
w(4) = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3))
end function simps38_weights_${k1}$
#:endfor
end submodule stdlib_quadrature_simps
