FIXME : Add documentation
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real, | dimension (:) | :: | x | |||
real | :: | avh | ||||
real, | dimension (:) | :: | y | |||
real, | dimension (:) | :: | dy | |||
real | :: | avdy | ||||
integer | :: | n | ||||
real, | dimension (:) | :: | a | |||
real, | dimension (:,:) | :: | c | |||
integer | :: | ic | ||||
real, | dimension (0:n+1,3) | :: | r | |||
real, | dimension (0:n+1,2) | :: | t | |||
integer | :: | ier |
subroutine spint1(x,avh,y,dy,avdy,n,a,c,ic,r,t,ier)
!
! initializes the arrays c, r and t for one dimensional cubic
! smoothing spline fitting by subroutine spfit1. the values
! df(i) are scaled so that the sum of their squares is n
! and the average of the differences x(i+1) - x(i) is calculated
! in avh in order to avoid underflow and overflow problems in
! spfit1.
!
! subroutine sets ier if elements of x are non-increasing,
! if n is less than 3, if ic is less than n-1 or if dy(i) is
! not positive for some i.
!
!---specifications for arguments---
integer n,ic,ier
! All variables should be double precision
real, dimension (:) :: x, y, dy, a
real, dimension (:,:) :: c
real, dimension (0:n+1,3) :: r
real, dimension (0:n+1,2) :: t
real :: avh, avdy
!
!---specifications for local variables---
integer i
real :: e,f,g,h
real :: zero = 0.
!
!---initialization and input checking---
ier = 0
if (n.lt.3) go to 60
if (ic.lt.n-1) go to 70
!
!---get average x spacing in avh---
g = zero
do i = 1,n - 1
h = x(i+1) - x(i)
if (h.le.zero) go to 80
g = g + h
end do
avh = g/ (n-1)
!
!---scale relative weights---
g = zero
do i = 1,n
if (dy(i).le.zero) go to 90
g = g + dy(i)*dy(i)
end do
avdy = sqrt(g/n)
dy = dy / avdy
!
!---initialize h,f---
h = (x(2)-x(1))/avh
f = (y(2)-y(1))/h
!
!---calculate a,t,r---
do i = 2,n - 1
g = h
h = (x(i+1)-x(i))/avh
e = f
f = (y(i+1)-y(i))/h
a(i) = f - e
t(i,1) = 2.0d0* (g+h)/3.0d0
t(i,2) = h/3.0d0
r(i,3) = dy(i-1)/g
r(i,1) = dy(i+1)/h
r(i,2) = -dy(i)/g - dy(i)/h
end do
!
!---calculate c = r'*r---
r(n,2) = zero
r(n,3) = zero
r(n+1,3) = zero
do i = 2,n - 1
c(i,1) = r(i,1)*r(i,1) + r(i,2)*r(i,2) + r(i,3)*r(i,3)
c(i,2) = r(i,1)*r(i+1,2) + r(i,2)*r(i+1,3)
c(i,3) = r(i,1)*r(i+2,3)
end do
return
!
!---error conditions---
60 ier = 130
return
70 ier = 129
return
80 ier = 131
return
90 ier = 132
return
end subroutine spint1