fitp_curvi Function

private function fitp_curvi(xl, xu, n, x, y, yp, sigma)

FIXME : Add documentation

Arguments

Type IntentOptional Attributes Name
real, intent(in) :: xl
real, intent(in) :: xu
integer, intent(in) :: n
real, intent(in), dimension(n) :: x
real, intent(in), dimension(n) :: y
real, intent(in), dimension(n) :: yp
real, intent(in) :: sigma

Return Value real


Contents

Source Code


Source Code

  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, dummy, dels, sum, del1, del2

    ! denormalize tension factor
    sigmap = abs(sigma)*real(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
        call fitp_snhcsh (dummy,cu1,sigmap*delu1,2)
        call fitp_snhcsh (dummy,cu2,sigmap*delu2,2)
        call fitp_snhcsh (dummy,cl1,sigmap*dell1,2)
        call fitp_snhcsh (dummy,cl2,sigmap*dell2,2)
        call fitp_snhcsh (ss,dummy,sigmap*dels,-1)
        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
          call fitp_snhcsh (dummy,c1,sigmap*del1,2)
          call fitp_snhcsh (dummy,c2,sigmap*del2,2)
          call fitp_snhcsh (ss,cs,sigmap*dels,3)
          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
            call fitp_snhcsh (ss,cs,sigmap*dels,3)
            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
          call fitp_snhcsh (dummy,c1,sigmap*del1,2)
          call fitp_snhcsh (dummy,c2,sigmap*del2,2)
          call fitp_snhcsh (ss,cs,sigmap*dels,3)
          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