FIXME : Add documentation (or tidyup above)
Type | Intent | Optional | 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 |
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