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

c The following produces piecewise monotonic interpolants rather than monotonicity preserving
c the shape of the original data.
c
c Steffen (1990) limits interpolated values so they never exceed the given values, 
c even in places where the given values are not monotonic.
c
c The downside is perhaps reduced accuracy on smooth data.
c
c Steffen, M., 1990: A simple method for monotonic interpolation in one dimension. 
c Astron. Astrophys., 239, 443-450.
c

      implicit none
      integer nobs, nint, i, L
      real x(nobs),y(nobs),xint(nint),yint(nint)
      real a(3)
      real xp(4),yp(4)
      do L=1,nint
          do i=2,nobs-2
             if(xint(L).gt.x(i).and.xint(L).le.x(i+1)) then
             xp(1)=x(i-1)
             xp(2)=x(i)
             xp(3)=x(i+1)
             xp(4)=x(i+2)
             yp(1)=y(i-1)
             yp(2)=y(i)
             yp(3)=y(i+1)
             yp(4)=y(i+2)
             call getcoeff(xp,yp,a)
             yint(L)=y(i)+a(1)*(xint(L)-x(i))+a(2)*(xint(L)-x(i))**2
     &                                     +a(3)*(xint(L)-x(i))**3
             goto 40
             endif
          enddo

 40       continue
      enddo
      return
      end

c23456789012346578901234657890123465789012346578901234657890123465789012


      subroutine getcoeff(x, y, a)

c returns coefficients a(3) for monotonic cubic interpolation from x(2) to x(3)
      real x(4)    ! junction points, strictly monotonic
      real y(4)    ! data values at x's
      real a(3)    ! coefficients
      real h1, h2, h3, s1, s2, s3, p2, p3, as2, ss2, yp2, yp3

c for x(2) <= xint <= x(3) and dx = xint-x(2), 
c y(x) = y(2) + dx*(a(1) + dx*(a(2) + dx*a(3)))
c where xint is the x value yint is being interpolated to

      h1 = x(2)-x(1)
      h2 = x(3)-x(2)
      h3 = x(4)-x(3)
      s1 = (y(2)-y(1))/h1
      s2 = (y(3)-y(2))/h2
      s3 = (y(4)-y(3))/h3
      p2 = (s1*h2+s2*h1)/(h1+h2)
      p3 = (s2*h3+s3*h2)/(h2+h3)
      as2 = abs(s2)
      ss2 = sign(1.0, s2)
      yp2 = (sign(1.0, s1)+ss2)*min(abs(s1), as2, 0.5*abs(p2))
      yp3 = (ss2+sign(1.0, s3))*min(as2, abs(s3), 0.5*abs(p3))
      a(1) = yp2
      a(2) = (3*s2-2*yp2-yp3)/h2
      a(3) = (yp2+yp3-2*s2)/h2**2
      return
      end