fitp_curvp1 Subroutine

public subroutine fitp_curvp1(n, x, y, p, yp, temp, sigma, ierr)

FIXME : Add documentation

Arguments

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

Contents

Source Code


Source Code

  subroutine fitp_curvp1 (n,x,y,p,yp,temp,sigma,ierr)
    implicit none
    integer, intent(in) :: n
    integer, intent(out) :: ierr
    real, intent(in) :: sigma, p
    real, dimension(:), intent(in) :: x, y
    real, dimension(:), intent(out) :: yp, temp
!!    real x(n),y(n),p,yp(n),temp(2*n),sigma
!      real x(n),y(n),p,yp(n),temp(1),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 a periodic interpolatory spline under tension
! through a sequence of functional values. for actual ends
! of the curve may be specified or omitted.  for actual
! computation of points on the curve it is necessary to call
! the function curvp2.
!
! 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) ).
!
!   p is the period (p .gt. x(n)-x(1)).
!
!   yp is an array of length at least n.
!
!   temp is an array of length at least 2*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 p is less than or equal to x(n)-x(1),
!        = 3 if x-values are not strictly increasing.
!
! and
!
!  n, x, y, and sigma are unaltered.
!
! this subroutine references package modules terms and
! snhcsh.
!
!-----------------------------------------------------------
    integer :: i, npibak, npi, ibak, nm1, np1
    real :: diag, diag2, sdiag2, ypn, dx2, sigmap, delx1
    real :: sdiag1, delx2, dx1, diag1

    nm1 = n-1
    np1 = n+1
    ierr = 0
    if (n .le. 1) go to 6
    if (p .le. x(n)-x(1) .or. p .le. 0.) go to 7
!
! denormalize tension factor
!
    sigmap = abs(sigma)*real(n)/p
!
! set up right hand side and tridiagonal system for yp and
! perform forward elimination
!
    delx1 = p-(x(n)-x(1))
    dx1 = (y(1)-y(n))/delx1
    call fitp_terms (diag1,sdiag1,sigmap,delx1)
    delx2 = x(2)-x(1)
    if (delx2 .le. 0.) go to 8
    dx2 = (y(2)-y(1))/delx2
    call fitp_terms (diag2,sdiag2,sigmap,delx2)
    diag = diag1+diag2
    yp(1) = (dx2-dx1)/diag
    temp(np1) = -sdiag1/diag
    temp(1) = sdiag2/diag
    dx1 = dx2
    diag1 = diag2
    sdiag1 = sdiag2
    if (n .eq. 2) go to 2
    do i = 2,nm1
       npi = n+i
       delx2 = x(i+1)-x(i)
       if (delx2 .le. 0.) go to 8
       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(npi) = -temp(npi-1)*sdiag1/diag
       temp(i) = sdiag2/diag
       dx1 = dx2
       diag1 = diag2
       sdiag1 = sdiag2
    end do
2   delx2 = p-(x(n)-x(1))
    dx2 = (y(1)-y(n))/delx2
    call fitp_terms (diag2,sdiag2,sigmap,delx2)
    yp(n) = dx2-dx1
    temp(nm1) = temp(2*n-1)-temp(nm1)
    if (n .eq. 2) go to 4
!
! perform first step of back substitution
    !
    do i = 3,n
       ibak = np1-i
       npibak =n+ibak
       yp(ibak) = yp(ibak)-temp(ibak)*yp(ibak+1)
       temp(ibak) =temp(npibak)-temp(ibak)*temp(ibak+1)
    end do
4   yp(n) = (yp(n)-sdiag2*yp(1)-sdiag1*yp(nm1))/ &
         (diag1+diag2+sdiag2*temp(1)+sdiag1*temp(nm1))
!
! perform second step of back substitution
!
    ypn =   yp(n)
    do i = 1,nm1
       yp(i) = yp(i)+temp(i)*ypn
    end do
    return
!
! too few points
!
6   ierr = 1
    return
!
! period too small
!
7   ierr = 2
    return
!
! x-values not strictly increasing
!
8   ierr = 3
    return
  end subroutine fitp_curvp1