spint1 Subroutine

public subroutine spint1(x, avh, y, dy, avdy, n, a, c, ic, r, t, ier)

FIXME : Add documentation

Arguments

Type IntentOptional 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

Contents

Source Code


Source Code

  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