mod_spline.f90

Go to the documentation of this file.

module mod_spline

  use itm_types
  use mod_profiles

  implicit none

  type(spline_coefficients) :: q_spline, p_spline, rbphi_spline

contains

  subroutine evaluate_spline(p, x, x_new, f_new, error, fprime_new)
!-----------------------------------------------------------------------
!     input:
!
!     p           arrays of spline coefficients (from spline)
!     x           array of knots for splines
!     x_new       array of knots at which the spline should be evaluated
!
!     output:
!
!     f_new       array of interpolated values
!     fprime_new  array of interpolated first derivatives
!     error       error code (0 = no error)
!-----------------------------------------------------------------------

    use itm_types

    implicit none

    type(spline_coefficients), intent(in) :: p
    real(r8), intent(in) :: x(:), x_new(:)
    real(r8), intent(inout) :: f_new(:)
    real(r8), intent(inout), optional :: fprime_new(:)
    integer(itm_i4), intent(out) :: error

    real(r8) :: x_loc(size(x))
    real(r8) :: x_new_loc(size(x_new))
    real(r8) :: xx, xwert_loc
    integer(itm_i4) ::  i, j, k, m
    logical :: inverted

    error = 1

!-- dimension checks
    if (size(x) /= size(p%sp1) .or. size(x_new) /= size(f_new)) then
      error = 1
      return
    end if

    if (present(fprime_new)) then
      if (size(x_new) /= size(fprime_new)) then
        error = 1
        return
      end if
    end if

!-- check whether abscissa has to be inverted
    if (x(size(x)) < x(1)) then
      inverted = .true.
    else
      inverted = .false.
    end if

!-- invert abscissa if necessary
    if (inverted) then
      x_loc = -x
      x_new_loc = -x_new
    else
      x_loc = x
      x_new_loc = x_new
    end if

    xwert_loc = x_new_loc(1)

!-- find matching interval (binary search)
    i = 1
    k = size(x_loc)

    m = (i + k) / 2

    do
      if (m == i) exit
      if (xwert_loc < x_loc(m)) then
        k = m
      else 
        i = m
      end if
      m = (i + k) / 2
    end do

    xx = xwert_loc - x_loc(i)
    f_new(1) = ((p%sp4(i) * xx + p%sp3(i)) * xx + p%sp2(i)) * xx + p%sp1(i)

    if (present(fprime_new)) &
     fprime_new(1) = (3.0_r8 * p%sp4(i) * xx + 2.0_r8 * p%sp3(i)) * xx &
     + p%sp2(i)

    do j = 2, size(x_new)
      xwert_loc = x_new_loc(j)

      do while (x_loc(k) <= xwert_loc)
        if (k == size(x_loc)) exit
        k = k + 1
        i = i + 1
      end do

      xx = xwert_loc - x_loc(i)
      f_new(j) = ((p%sp4(i) * xx + p%sp3(i)) * xx + p%sp2(i)) * xx + p%sp1(i)

      if (present(fprime_new)) &
       fprime_new(j) = (3.0_r8 * p%sp4(i) * xx + 2.0_r8 * p%sp3(i)) * xx &
       + p%sp2(i)
    end do

!-- adapt derivative if necessary
    if (inverted .and. present(fprime_new)) then
      fprime_new = -fprime_new
    end if

    return

  end subroutine evaluate_spline

end module mod_spline