fitp_curvi Function

  real function fitp_curvi (xl,xu,n,x,y,yp,sigma)
    implicit none
    integer, intent(in) :: n
    real, intent(in) :: xl, xu, sigma
    real, dimension(n), intent(in) :: x, y, yp
!
!                                 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 integrates a curve specified by a spline
! under tension between two given limits. the subroutine
! curv1 should be called earlier to determine necessary
! parameters.
!
! on input--
!
!   xl and xu contain the upper and lower limits of inte-
!   gration, respectively. (sl need not be less than or
!   equal to xu, curvi (xl,xu,...) .eq. -curvi (xu,xl,...) ).
!
!   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 from subroutine curv1 containing
!   the values of the second derivatives 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--
!
!   curvi contains the integral value.
!
! none of the input parameters are altered.
!
! this function references package modules intrvl and
! snhcsh.
!
!-----------------------------------------------------------
    integer :: i, ilp1, ilm1, il, ium1, iu
    real :: delu1, delu2, c2, ss, cs, cu2, cl1, cl2, cu1
    real :: dell1, dell2, deli, c1, ssign, sigmap
    real :: xxl, xxu, t1, t2, dels, sum, del1, del2

    ! denormalize tension factor
    sigmap = abs(sigma) * (n - 1) / (x(n) - x(1))

    ! determine actual upper and lower bounds
    if (xl < xu) then
       xxl = xl
       xxu = xu
       ssign = 1.
    else if (xl > xu) then
       xxl = xu
       xxu = xl
       ssign = -1.
    else
       ! return zero if xl .eq. xu
       fitp_curvi = 0.
       return
    end if

    ! search for proper intervals
    ilm1 = fitp_intrvl (xxl,x,n)
    il = ilm1+1
    ium1 = fitp_intrvl (xxu,x,n)
    iu = ium1+1
    if (il == iu) then
       ! integrate from xxl to xxu
       delu1 = xxu-x(ium1)
       delu2 = x(iu)-xxu
       dell1 = xxl-x(ium1)
       dell2 = x(iu)-xxl
       dels = x(iu)-x(ium1)
       deli = xxu-xxl
       t1 = (delu1+dell1)*deli/(2.*dels)
       t2 = (delu2+dell2)*deli/(2.*dels)
       sum = t1*y(iu)+t2*y(ium1)
       if (is_not_zero(sigma)) then
          cu1 = coshmm_fun(sigmap*delu1)
          cu2 = coshmm_fun(sigmap*delu2)
          cl1 = coshmm_fun(sigmap*dell1)
          cl2 = coshmm_fun(sigmap*dell2)
          ss = sinhm_fun(sigmap*dels)

          sum = sum+(yp(iu)*(delu1*delu1*(cu1-ss/2.) &
               -dell1*dell1*(cl1-ss/2.)) &
               +yp(ium1)*(dell2*dell2*(cl2-ss/2.) &
               -delu2*delu2*(cu2-ss/2.)))/ &
               (sigmap*sigmap*dels*(1.+ss))
       else
          sum = sum-t1*(delu2*(dels+delu1)+dell2*(dels+dell1))* &
               yp(iu)/12. &
               -t2*(dell1*(dels+dell2)+delu1*(dels+delu2))* &
               yp(ium1)/12.
       end if
    else
       ! integrate from xxl to x(il)
       sum = 0.
       if (not_exactly_equal(xxl, x(il))) then
          del1 = xxl-x(ilm1)
          del2 = x(il)-xxl
          dels = x(il)-x(ilm1)
          t1 = (del1+dels)*del2/(2.*dels)
          t2 = del2*del2/(2.*dels)
          sum = t1*y(il)+t2*y(ilm1)
          if (is_not_zero(sigma)) then
             c1 = coshmm_fun(sigmap*del1)
             c2 = coshmm_fun(sigmap*del2)
             cs = coshmm_fun(sigmap*dels)
             ss = sinhm_fun(sigmap*dels)
             sum = sum+((dels*dels*(cs-ss/2.)-del1*del1*(c1-ss/2.)) &
                  *yp(il)+del2*del2*(c2-ss/2.)*yp(ilm1))/ &
                  (sigmap*sigmap*dels*(1.+ss))
          else
             sum = sum-t1*t1*dels*yp(il)/6. &
                  -t2*(del1*(del2+dels)+dels*dels)*yp(ilm1)/12.
          end if
       end if

       ! integrate over interior intervals
       if (iu-il /= 1) then
          ilp1 = il+1
          do i = ilp1,ium1
             dels = x(i)-x(i-1)
             sum = sum+(y(i)+y(i-1))*dels/2.
             if (is_not_zero(sigma)) then
                cs = coshmm_fun(sigmap*dels)
                ss = sinhm_fun(sigmap*dels)
                sum = sum+(yp(i)+yp(i-1))*dels*(cs-ss/2.)/(sigmap*sigmap*(1.+ss))
             else
                sum = sum-(yp(i)+yp(i-1))*dels*dels*dels/24.
             end if
          end do
       end if

       ! integrate from x(iu-1) to xxu
       if (not_exactly_equal(xxu, x(ium1))) then
          del1 = xxu-x(ium1)
          del2 = x(iu)-xxu
          dels = x(iu)-x(ium1)
          t1 = del1*del1/(2.*dels)
          t2 = (del2+dels)*del1/(2.*dels)
          sum = sum+t1*y(iu)+t2*y(ium1)
          if (is_not_zero(sigma)) then
             c1 = coshmm_fun(sigmap*del1)
             c2 = coshmm_fun(sigmap*del2)
             cs = coshmm_fun(sigmap*dels)
             ss = sinhm_fun(sigmap*dels)
             sum = sum+(yp(iu)*del1*del1*(c1-ss/2.)+yp(ium1)* &
                  (dels*dels*(cs-ss/2.)-del2*del2*(c2-ss/2.))) &
                  /(sigmap*sigmap*dels*(1.+ss))
          else
             sum = sum-t1*(del2*(del1+dels)+dels*dels)*yp(iu)/12.-t2*t2*dels*yp(ium1)/6.
          end if
       end if
    end if

    ! correct sign and return
    fitp_curvi = ssign*sum
  end function fitp_curvi