module quad

implicit none

contains

  real function quadintflg(x, y, xval, flag, xflg, yflg)

    ! Performs quadratic interpolation.
    ! flag returned as 'OK' for success, '1' if outside x-range,
    ! and hence extrapolated

    real, dimension(:), intent(in) :: x, y
    character(len=*), dimension(:), intent(in), optional :: xflg, yflg
    real, intent(in) :: xval
    character(len=*), intent(out) :: flag

    integer :: i, j
    real :: a, b, c

    ! Start as we hope to go on.
    flag='OK'

    if (size(x) == 2) then
      ! Only enough data for linear interpolation.
      quadintflg=y(1) + (y(2)-y(1))*(x(1)-xval)/(x(1)-x(2))
      if (present(xflg)) then
        do j=1, 2
          if (xflg(j) /= 'OK') flag=xflg(j)
        end do
      end if
      if (present(yflg)) then
        do j=1, 2
          if (yflg(j) /= 'OK') flag=yflg(j)
        end do
      end if
    else
      i=minloc(abs(x-xval),1) 
      ! Always interpolate between the two points either side and the
      ! point whose index is one more (this is for compatibility with 
      ! the numerical recipes routine we used to use, and to ensure 
      ! we do not change the interpolating polynomial between points).
      if ((xval-x(i))/(x(size(x))-x(1)) > 0.0) i=i+1
      ! Don't run off the top.
      if (i >= size(x)) i=size(x)-1
      ! Or the bottom.
      if (i <= 1) i=2
      a = (y(i-1)-y(i))/(x(i-1)-x(i)) - (y(i)-y(i+1))/(x(i)-x(i+1))
      a = a/((x(i-1)*x(i-1)-x(i)*x(i))/(x(i-1)-x(i)) - &
      (x(i)*x(i)-x(i+1)*x(i+1))/(x(i)-x(i+1)))
      b = (y(i-1)-y(i)-a*(x(i-1)*x(i-1)-x(i)*x(i)))/(x(i-1)-x(i))
      c = y(i-1) - a*x(i-1)*x(i-1) - b*x(i-1)
      quadintflg = a*xval*xval + b*xval + c
      if (present(xflg)) then
        do j=i-1, i+1
          if (xflg(j) /= 'OK') flag=xflg(j)
        end do
      end if
      if (present(yflg)) then
        do j=i-1, i+1
          if (yflg(j) /= 'OK') flag=yflg(j)
        end do
      end if
    end if

    ! Did we extrapolate?
    if (xval < minval(x)) then
      i=minloc(x,1)
      quadintflg=y(i)
      flag='1'
    else if (xval > maxval(x)) then
      i=maxloc(x,1)
      quadintflg=y(i)
      flag='1'
    end if

  end function quadintflg

  real function quadint(x, y, xval, flag)

    real, dimension(:), intent(in) :: x, y
    real, intent(in) :: xval
    character(len=*), intent(out) :: flag

    character(len=10), allocatable, dimension(:) :: xflg, yflg

    allocate(xflg(size(x)), yflg(size(y)))

    xflg='OK'
    yflg='OK'

    quadint=quadintflg(x, y, xval, flag, xflg, yflg)

    deallocate(xflg, yflg)

  end function quadint

  real function linint(x, y, xval, flag, xflg, yflg)

    real, dimension(:), intent(in) :: x, y
    ! Note these are interchangeable.  So if you only have one set of
    ! flags you can call it in just one set.
    character(len=*), dimension(:), intent(in), optional :: xflg, yflg
    real, intent(in) :: xval
    character(len=*), intent(out) :: flag
    
    integer :: i

    flag='OK'

    i=minloc(abs(x-xval),1) 
    ! Make i the data point above xval.
    if ((xval-x(i))/(x(size(x))-x(1)) > 0.0) i=i+1
    ! Don't run off the top.
    if (i >= size(x)) i=size(x)
    ! Or the bottom.
    if (i <= 1) i=2
    linint=y(i-1) + (y(i)-y(i-1))*(x(i-1)-xval)/(x(i-1)-x(i))

    ! Was either datapoint flagged?
    if (present(xflg)) then
      if (xflg(i-1) /= 'OK') flag=xflg(i-1)
      if (xflg(i) /= 'OK') flag=xflg(i)
    end if
    if (present(yflg)) then
      if (yflg(i-1) /= 'OK') flag=yflg(i-1)
      if (yflg(i) /= 'OK') flag=yflg(i)
    end if

    ! Did we extrapolate?
    if (xval < minval(x)) then
      flag='1'
    else if (xval > maxval(x)) then
      flag='1'
    end if

  end function linint

end module quad

