fitp_curvd Function

  real function fitp_curvd (t,n,x,y,yp,sigma)
    implicit none
    integer, intent(in) :: n
    real, dimension(n), intent(in) :: x, y, yp
    real, intent(in) :: t, sigma
!
!                                 coded by alan kaylor cline
!                           from fitpack -- january 26, 1987
!                        a curve and surface fitting package
!                      a product of pleasant valley software
!                  8603 altus cove, austin, texas 78759, usa
!
! this function differentiates a curve at a given point
! using a spline under tension. the subroutine curv1 should
! be called earlier to determine certain necessary
! parameters.
!
! on input--
!
!   t contains a real value at which the derivative is to be
!   determined.
!
!   n contains the number of points which were specified to
!   determine the curve.
!
!   x and y are arrays containing the abscissae and
!   ordinates, respectively, of the specified points.
!
!   yp is an array of second derivative values of the curve
!   at the nodes.
!
! and
!
!   sigma contains the tension factor (its sign is ignored).
!
! the parameters n, x, y, yp, and sigma should be input
! unaltered from the output of curv1.
!
! on output--
!
!   curvd contains the derivative value.
!
! none of the input parameters are altered.
!
! this function references package modules intrvl and
! snhcsh.
!
!-----------------------------------------------------------
    integer :: i, im1
    real :: ss, sigdel, dummy, c1, c2, sum, sigmap
    real :: del1, del2, dels
!
! determine interval
!
    im1 = fitp_intrvl(t,x,n)
    i = im1+1
!
! denormalize tension factor
!
    sigmap = abs(sigma)*real(n-1)/(x(n)-x(1))
!
! set up and perform differentiation
!
    del1 = t-x(im1)
    del2 = x(i)-t
    dels = x(i)-x(im1)
    sum = (y(i)-y(im1))/dels
    if (is_zero(sigmap)) then
      fitp_curvd = sum+(yp(i)*(2.*del1*del1-del2*(del1+dels))- &
           yp(im1)*(2.*del2*del2-del1*(del2+dels))) &
           /(6.*dels)
      return
    end if
    sigdel = sigmap*dels
    call fitp_snhcsh (ss,dummy,sigdel,-1)
    call fitp_snhcsh (dummy,c1,sigmap*del1,1)
    call fitp_snhcsh (dummy,c2,sigmap*del2,1)
    fitp_curvd = sum+(yp(i)*(c1-ss)-yp(im1)*(c2-ss))/(sigdel*sigmap*(1.+ss))
    return
  end function fitp_curvd