fitp_curv1 Subroutine

public pure subroutine fitp_curv1(n, x, y, slp1, slpn, islpsw, yp, temp, sigma, ierr)

FIXME : Add documentation (or tidyup above)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n
real, intent(in), dimension(n) :: x
real, intent(in), dimension(n) :: y
real, intent(in) :: slp1
real, intent(in) :: slpn
integer, intent(in) :: islpsw
real, intent(out), dimension(n) :: yp
real, intent(out), dimension(n) :: temp
real, intent(in) :: sigma
integer, intent(out) :: ierr

Contents

Source Code


Source Code

  pure subroutine fitp_curv1 (n,x,y,slp1,slpn,islpsw,yp,temp,sigma,ierr)
    implicit none
    integer, intent(in) :: n, islpsw
    integer, intent(out) :: ierr
    real, dimension(n), intent(in) :: x, y
    real, dimension(n), intent(out) :: yp, temp
    real, intent(in) :: slp1,slpn,sigma
!
!                                 coded by alan kaylor cline
!                           from fitpack -- january 26, 1987
!                        a curve and surface fitting package
!                      a product of pleasant valley software
!                  8603 altus cove, austin, texas 78759, usa
!
! this subroutine determines the parameters necessary to
! compute an interpolatory spline under tension through
! a sequence of functional values. the slopes at the two
! ends of the curve may be specified or omitted.  for actual
! computation of points on the curve it is necessary to call
! the function curv2.
!
! on input--
!
!   n is the number of values to be interpolated (n.ge.2).
!
!   x is an array of the n increasing abscissae of the
!   functional values.
!
!   y is an array of the n ordinates of the values, (i. e.
!   y(k) is the functional value corresponding to x(k) ).
!
!   slp1 and slpn contain the desired values for the first
!   derivative of the curve at x(1) and x(n), respectively.
!   the user may omit values for either or both of these
!   parameters and signal this with islpsw.
!
!   islpsw contains a switch indicating which slope data
!   should be used and which should be estimated by this
!   subroutine,
!          = 0 if slp1 and slpn are to be used,
!          = 1 if slp1 is to be used but not slpn,
!          = 2 if slpn is to be used but not slp1,
!          = 3 if both slp1 and slpn are to be estimated
!              internally.
!
!   yp is an array of length at least n.
!
!   temp is an array of length at least n which is used for
!   scratch storage.
!
! and
!
!   sigma contains the tension factor. this value indicates
!   the curviness desired. if abs(sigma) is nearly zero
!   (e.g. .001) the resulting curve is approximately a
!   cubic spline. if abs(sigma) is large (e.g. 50.) the
!   resulting curve is nearly a polygonal line. if sigma
!   equals zero a cubic spline results.  a standard value
!   for sigma is approximately 1. in absolute value.
!
! on output--
!
!   yp contains the values of the second derivative of the
!   curve at the given nodes.
!
!   ierr contains an error flag,
!        = 0 for normal return,
!        = 1 if n is less than 2,
!        = 2 if x-values are not strictly increasing.
!
! and
!
!   n, x, y, slp1, slpn, islpsw and sigma are unaltered.
!
! this subroutine references package modules ceez, terms,
! and snhcsh.
!
!-----------------------------------------------------------
    integer :: i, ibak, nm1, np1
    real :: sdiag1, diag1, delxnm, dx1, diag, sdiag2, dx2, diag2
    real :: delxn, slpp1, delx1, sigmap, c3, c2, c1, slppn, delx2
    
    nm1 = n-1
    np1 = n+1
    ierr = 0
    if (n .le. 1) go to 8
    if (x(n) .le. x(1)) go to 9
!
! denormalize tension factor
!
    sigmap = abs(sigma) * (n - 1) / (x(n) - x(1))
!
! approximate end slopes
!
    if (islpsw .ge. 2) go to 1
    slpp1 = slp1
    go to 2
1   delx1 = x(2)-x(1)
    delx2 = delx1+delx1
    if (n .gt. 2) delx2 = x(3)-x(1)
    if (delx1 .le. 0. .or. delx2 .le. delx1) go to 9
    call fitp_ceez (delx1,delx2,sigmap,c1,c2,c3,n)
    slpp1 = c1*y(1)+c2*y(2)
    if (n .gt. 2) slpp1 = slpp1+c3*y(3)
2   if (islpsw .eq. 1 .or. islpsw .eq. 3) go to 3
    slppn = slpn
    go to 4
3   delxn = x(n)-x(nm1)
    delxnm = delxn+delxn
    if (n .gt. 2) delxnm = x(n)-x(n-2)
    if (delxn .le. 0. .or. delxnm .le. delxn) go to 9
    call fitp_ceez (-delxn,-delxnm,sigmap,c1,c2,c3,n)
    slppn = c1*y(n)+c2*y(nm1)
    if (n .gt. 2) slppn = slppn+c3*y(n-2)
!
! set up right hand side and tridiagonal system for yp and
! perform forward elimination
!
4   delx1 = x(2)-x(1)
    if (delx1 .le. 0.) go to 9
    dx1 = (y(2)-y(1))/delx1
    call fitp_terms (diag1,sdiag1,sigmap,delx1)
    yp(1) = (dx1-slpp1)/diag1
    temp(1) = sdiag1/diag1
    if (n .eq. 2) go to 6
    do i = 2,nm1
       delx2 = x(i+1)-x(i)
       if (delx2 .le. 0.) go to 9
       dx2 = (y(i+1)-y(i))/delx2
       call fitp_terms (diag2,sdiag2,sigmap,delx2)
       diag = diag1+diag2-sdiag1*temp(i-1)
       yp(i) = (dx2-dx1-sdiag1*yp(i-1))/diag
       temp(i) = sdiag2/diag
       dx1 = dx2
       diag1 = diag2
       sdiag1 = sdiag2
    end do
6   diag = diag1-sdiag1*temp(nm1)
    yp(n) = (slppn-dx1-sdiag1*yp(nm1))/diag
!
! perform back substitution
!
    do i = 2,n
       ibak = np1-i
       yp(ibak) = yp(ibak)-temp(ibak)*yp(ibak+1)
    end do
    return
!
! too few points
!
8   ierr = 1
    return
!
! x-values not strictly increasing
!
9   ierr = 2
    return
  end subroutine fitp_curv1