spline Subroutine

public subroutine spline(N_, x_, y_, b_, c_, d_)

determines b, c and d coefficients of the cubic spline function y(x) = y_i + b_i * dx + c_i * dx^2 + d_i * dx^3, where dx = x - x_i, and x_i <= x < x_i+1. The algorithm is based on Gerald, C., and Wheatley, P., "Applied Numerical Analysis", Addison-Wesley, 1994.

Arguments

Type IntentOptional Attributes Name
integer(kind=int32), intent(in) :: N_

number of grid points
real(kind=dp), intent(in) :: x_(N_)

grid points (ascending order)
real(kind=dp), intent(in) :: y_(N_)

tabulated values
real(kind=dp), intent(out) :: b_(N_)

arrays with coefficients of the spline function
real(kind=dp), intent(out) :: c_(N_)

arrays with coefficients of the spline function
real(kind=dp), intent(out) :: d_(N_)

arrays with coefficients of the spline function

Calls

proc~~spline~~CallsGraph proc~spline spline proc~integer_to_character integer_to_character proc~spline->proc~integer_to_character proc~write_error write_error proc~spline->proc~write_error proc~write_message write_message proc~write_error->proc~write_message

Called by

proc~~spline~~CalledByGraph proc~spline spline proc~interpolate_radial_coupling_terms interpolate_radial_coupling_terms proc~interpolate_radial_coupling_terms->proc~spline program~scattering SCATTERING program~scattering->proc~interpolate_radial_coupling_terms

Contents

Variables

Source Code

spline

Variables

Type Visibility Attributes Name Initial
real(kind=dp), public :: diff_x_(N_-1)
integer(kind=int32), public :: i_
integer(kind=int32), public :: j_
real(kind=dp), public :: w_

Source Code

      subroutine spline (N_, x_, y_, b_, c_, d_)
         !! determines b, c and d coefficients of the cubic spline function
         !! y(x) = y_i + b_i * dx + c_i * dx^2 + d_i * dx^3,
         !! where dx = x - x_i, and x_i <= x < x_i+1.
         !! The algorithm is based on
         !! Gerald, C., and Wheatley, P.,
         !! "Applied Numerical Analysis", Addison-Wesley, 1994.
         !---------------------------------------------------------------------!
         integer(int32), intent(in) :: N_
            !! number of grid points
         real(dp), intent(in) :: x_(N_)
            !! grid points (ascending order)
         real(dp), intent(in) :: y_(N_)
            !! tabulated values
         real(dp), intent(out) :: b_(N_), c_(N_), d_(N_)
            !! arrays with coefficients of the spline function
         !---------------------------------------------------------------------!
         integer(int32) :: i_, j_
         real(dp) :: w_
         real(dp) :: diff_x_(N_-1)
         !---------------------------------------------------------------------!
         ! check if the number of points is larger than 4
         !---------------------------------------------------------------------!
         if (N_ < 4) then
            call write_error("spline function called with " //                 &
               trim(adjustl(integer_to_character(N_))) // " points")
         endif
         !---------------------------------------------------------------------!
         ! check if x is sorted in ascending order
         !---------------------------------------------------------------------!
         do i_ = 2, N_
            if (x_(i_) <= x_(i_-1)) then
               call write_error("spline: x values are not in ascending order " &
                  // " at index " // trim(adjustl(integer_to_character(i_))))
            endif
         end do
         !---------------------------------------------------------------------!
         diff_x_ = x_(2:N_) - x_(1:N_-1)
         !---------------------------------------------------------------------!
         b_(2:N_-1) = 2.0_dp * (diff_x_(1:N_-2) + diff_x_(2:N_-1))
         b_(1)  = - diff_x_(1)
         b_(N_) = - diff_x_(N_-1)

         c_(2:N_-1) = ( y_(3:N_) - y_(2:N_-1) ) / diff_x_(2:N_-1)              &
            - ( y_(2:N_-1) - y_(1:N_-2) ) / diff_x_(1:N_-2)

         c_(1)  = c_(3)/(x_(4)-x_(2)) - c_(2)/(x_(3)-x_(1))
         c_(N_) = c_(N_-1)/(x_(N_)-x_(N_-2)) - c_(N_-2)/(x_(N_-1)-x_(N_-3))

         c_(1)  = c_(1)/(x_(4)-x_(1))*diff_x_(1)**2
         c_(N_) = -c_(N_)/(x_(N_)-x_(N_-3))*diff_x_(N_-1)**2

         do i_ = 2, N_
            w_ = diff_x_(i_-1)/b_(i_-1)
            b_(i_) = b_(i_) - w_*diff_x_(i_-1)
            c_(i_) = c_(i_) - w_*c_(i_-1)
         end do

         c_(N_) = c_(N_) / b_(N_)
         do j_ = 1, N_-1
            i_ = N_-j_
            c_(i_) = (c_(i_) - diff_x_(i_)*c_(i_+1)) / b_(i_)
         end do

         b_(1:N_-1) = ( y_(2:N_) - y_(1:N_-1) ) / diff_x_(1:N_-1)              &
            - ( 2.0_dp * c_(1:N_-1) + c_(2:N_) ) * diff_x_(1:N_-1)

         d_(1:N_-1) = ( c_(2:N_) - c_(1:N_-1) ) / diff_x_(1:N_-1)

         c_ = c_ * 3.0_dp
           
      end subroutine spline