!23456789012345678901234567890123456789012345678901234567890123456789012

      subroutine steffen(nobs, x, y, nint, xint, yint)
      implicit none

! STEFFEN 1-D Steffen interpolation
!
!    STEFFEN(NOBS,X,Y,NINT,XINT,YINT) interpolates to find YINT,
!    Input is NOBS, X, Y, NINT, and XINT. 
!
!    X and Y are observed data of length NOBS from which
!          a polynomial spline is formulated. X must be increasing
!          in the array. The code will stop with a warning if X(i) >= X(i+1).
!
!    XINT is an array of values with length NINT to compute the
!         inteprolated value YINT. XINT does not need to be increasing
!         in the array, but XINT cannot < X(1) or > X(NINT). In other words,
!         no extrapolation is allowed, and a missing value (MISS) is generated
!         with a warning message. Here, MISS=-999, but the user can change.
!
!    Steffen's interpolation method is based on a third-order polynomial.
!    The slope at each grid point is calculated in a way to guarantee
!    a monotonic behavior of the interpolating function. The curve 
!    is smooth up to the first derivative.

!    Steffen’s method guarantees the monotonicity of the interpolating
!    function between the given data points. Steffen's formulation is the
!    same as Akima (1970) for Eqs. 1-7. See:
!
!    https://quantmacro.wordpress.com/2015/09/01/akima-spline-interpolation-in-excel/
!
!    and Eq. 8, which computes the slope based on three points of known
!    y value and known slope value --- a technique known as the osculatory
!    interpolation method (see discussion on pg 592 in Akima (1970), his Eq. 10).
!
!    But, Steffen provides different constraints on the polynomial than Akima. 
!    Akima's constraint provided mostly smooth curves, but occassional wiggles
!    emerged, especially for adjacent data points of the same value. Steffen's
!    constrain imposes that minima and maxima can only occur exactly at
!    the data points, and there can never be spurious oscillations between
!    data points. The interpolated function is piecewise cubic in each
!    interval. The resulting curve and its first derivative are guaranteed
!    to be continuous, but the second derivative may be discontinuous.

!    This code is based on Joe Henning's original MATLAB program, written
!    in Summer 2014.

!    Yee Lau (lau@gri.msstate.edu) converted his code to FORTRAN
!    in January 2018. 
!
!    Pat Fitzpatrick (fitz@gri.msstate.edu) added a driver routine and
!    provided extensive comments in June 2018.
!
! References:
!
! Steffen, M., 1990: A simple method for monotonic interpolation
! in one dimension. Astron. Astrophys. 239, 443-450.
!
! Akima, H., 1970: A new method of interpolation and smooth curve fitting
! based on local procedure. J. Ass. Computing Machinery, 17, 589-602.
! 

      real x(nobs), y(nobs), yp(nobs)
      real xint(nint), yint(nint) 
      real hi, him1, si, sim1, pi
      real h, h1, h2, s, s1, s2, p1, hnm1, hnm2, snm1, snm2
      real t, t2, t3, a, b, c, d, pn, miss
      integer nobs, nint, iargc, i, k, isiny, klo, khi

      parameter(miss=-999)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Check input x array for monotonic increasing
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      do i = 1, nobs-1
         if (x(i) >= x(i+1)) then    
             print *, 'Error! Input x is not monotonic increasing. ' //
     &                'Program stop.'
             return
         endif
      enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The next section calculate slopes for array yp for internal points
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      do i = 2, nobs-1

! Eq. 6 in Steffen (1990)

         hi = x(i+1) - x(i)
         him1 = x(i) - x(i-1)

! Eq. 7 in Steffen (1990)

         si = (y(i+1) - y(i))/hi
         sim1 = (y(i) - y(i-1))/him1

! Eq. 8 in Steffen (1990).
!
! This computes the slope at an internal point
! based on the slope of the secant through points
! x(i-1), y(i-1); x(i), y(i); and x(i+1), y(i+1)
!
! This is known as osculatory interpolation, an old
! but still effective interpolation technique
! as discussed in Akima (1970).
!
! Steffen states it can be derived from Bessel
! intepolation, but the text lacks details.
!

         pi = (sim1*hi + si*him1)/(him1 + hi)

! Eq. 11 in Steffen (1990) applies the constraint
! to the slope, stating "the slope at i must not be
! greater than twice the minimum of the absolute values
! of the slopes given by the left- and right-handed
! finite differences. If these have different signs,
! (the slope) must be zero."
!
! Effectively, minima and maxima can only occur exactly
! at the data points, and there can never be spurious
! oscillations between data points. 

         if (sim1*si <= 0) then
            yp(i) = 0

         elseif ((abs(pi) > 2.*abs(sim1)) .or. 
     &          (abs(pi) > 2.*abs(si))) then
            if (sign(sim1,si) /= sim1) then
                print *, 'Error: Sign si not equal to sim1'
                stop
            endif
            if (sim1 > 0) then
                a = 1
            elseif (sim1 == 0) then
                a = 0
            else
                a = -1
            endif
            yp(i) = 2.*a*min1(abs(sim1),abs(si))
         else
            yp(i) = pi
         endif
      enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The next section calculate slopes for array yp at boundaries
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Apply boundary condition at the first point

      h1 = x(2) - x(1)
      h2 = x(3) - x(2)
      s1 = (y(2) - y(1))/h1
      s2 = (y(3) - y(2))/h2

! Eq. 24 in Steffen (1990)

      p1 = s1*(1. + h1/(h1 + h2)) - s2*h1/(h1 + h2)

! Eq. 26 in Steffen (1990)

      if (p1*s1 <= 0) then
         yp(1) = 0
      elseif (abs(p1) > 2*abs(s1)) then
         yp(1) = 2.*s1
      else
         yp(1) = p1
      endif

! Apply boundary condition at the last point

      hnm1 = x(nobs) - x(nobs-1)
      hnm2 = x(nobs-1) - x(nobs-2)
      snm1 = (y(nobs) - y(nobs-1))/hnm1
      snm2 = (y(nobs-1) - y(nobs-2))/hnm2

! Eq. 25 in Steffen (1990)

      pn = snm1*(1. + hnm1/(hnm1 + hnm2)) - snm2*hnm1/(hnm1 + hnm2)

! Eq. 27 in Steffen (1990)

      if (pn*snm1 <= 0) then
         yp(nobs) = 0
      elseif (abs(pn) > 2.*abs(snm1)) then
         yp(nobs) = 2*snm1
      else
         yp(nobs) = pn
      endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Perform interpolation at array xint
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      do i = 1, nint

! Check the current xint value is within range

         if (xint(i) < x(1) .or. xint(i) > x(nobs)) then
             print *, 'Warning! ' //
     &          'Input xint =',xint(i), ' is out of range. ' //
     &          'Set interpolated value to missing.'
             yint(i) = miss        ! go to next i in do loop
             cycle
         endif

! Find the interpolation interval using bisection

         klo = 1
         khi = nobs
         do while (khi-klo > 1) 
            k = floor((khi+klo)/2.0)
            if (x(k) > xint(i)) then
               khi = k
            else
               klo = k
            endif
         enddo

! Eq. 6 in Steffen (1990)

         h = x(khi) - x(klo)

         if (h == 0.0) then
             print *, 'Bad x input  ==> x values must be distinct. ' //
     &                'Set interpolated value to missing.'
             yint(i) = miss
             cycle                 ! go to next i in do loop
         endif

! If xint=x, then set yint=y. Also set isiny=1.
! Then go to next i in do loop

         isiny = 0
         do k = 1, nobs
           if (xint(i) == x(k)) then
              yint(i) = y(k)
              isiny = 1
              exit
           endif
         enddo

         if (isiny == 1) then
            cycle
         endif

! Perform interpolation

! Eq. 7 in Steffen (1990)

         s = (y(khi) - y(klo))/h

! Eqs. 2-5 in Steffen (1990)

         a = (yp(klo) + yp(khi) - 2.*s)/h/h
         b = (3.*s - 2.*yp(klo) - yp(khi))/h
         c = yp(klo)
         d = y(klo)

! Eq. 1 in Steffen (1990). This is the interpolation computation.
   
         t = xint(i) - x(klo)
         t2 = t*t
         t3 = t2*t
         yint(i) = a*t3 + b*t2 + c*t + d

      enddo

      return
      end