FIXME : Add documentation
Type | Intent | Optional | 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 |
pure 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
!
! 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) * 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