fitp_curvpi Function

  real function fitp_curvpi (xl,xu,n,x,y,p,yp,sigma)
    implicit none
    integer, intent(in) :: n
    real, intent(in) :: xl, xu, p, 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 periodic
! spline under tension between two given limits. the
! subroutine curvp1 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, curvpi (xl,xu,...) .eq. -curvpi (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.
!
!   p contains the period.
!
!   yp is an array from subroutine curvp1 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, p, yp, and sigma should be input
! unaltered from the output of curvp1.
!
! on output--
!
!
!   curvpi contains the integral value.
!
! none of the input parameters are altered.
!
! this function references package modules intrvp and
! snhcsh.
!
!--------------------------------------------------------------
    integer :: np1, im1, ii, iup1, ilp1, ideltp
    real :: s7, s6, s3, c2, c1, s5, s4, cl1, cu2, cu1, si, so
    real :: cl2, delu2, delu1, s8, deli, dell2, dell1
    integer :: ium1, il, isave, isign, lper, ilm1, iu, i
    real :: xxu, xsave, x1pp, sigmap, xxl, xil, del1, s2, cs
    real :: t2, t1, del2, s1, xiu, ss, dels

    integer :: uper
    logical :: bdy
!
! denormalize tension factor
!
    sigmap = abs(sigma) * n / p
!
! determine actual upper and lower bounds
!
    x1pp = x(1)+p
    isign = 1
    xxl = fitp_periodic_wrap_value(xl, x(1), p)
    ilm1 = fitp_intrvp (xxl, x, n)
    lper = int((xl-x(1))/p)
    if (xl .lt. x(1)) lper = lper-1
    xxu = fitp_periodic_wrap_value(xu, x(1), p)
    ium1 = fitp_intrvp (xxu, x, n)
    uper = int((xu-x(1))/p)
    if (xu .lt. x(1)) uper = uper-1
    ideltp = uper-lper
    bdy = ideltp * (xxu-xxl) .lt. 0.
    if ((ideltp .eq. 0 .and. xxu .lt. xxl) .or. ideltp .lt. 0) isign = -1
    if (bdy) ideltp = ideltp-isign
    if (xxu .ge. xxl) go to 1
    xsave = xxl
    xxl = xxu
    xxu = xsave
    isave = ilm1
    ilm1 = ium1
    ium1 = isave
1   il = ilm1+1
    if (ilm1 .eq. n) il = 1
    xil = x(il)
    if (ilm1 .eq. n) xil = x1pp
    iu = ium1+1
    if (ium1 .eq. n) iu = 1
    xiu = x(iu)
    if (ium1 .eq. n) xiu = x1pp
    s1 = 0.
    if (ilm1 .eq. 1 .or. (ideltp .eq. 0 .and..not. bdy)) go to 4
!
! integrate from x(1) to x(ilm1), store in s1
!
    do i = 2,ilm1
       dels = x(i)-x(i-1)
       s1 = s1+(y(i)+y(i-1))*dels/2.
       if (is_not_zero(sigma)) then
          ss = sinhm_fun(sigmap * dels)
          cs = coshmm_fun(sigmap * dels)
          s1 = s1+(yp(i)+yp(i-1))*dels*(cs-ss/2.)/(sigmap*sigmap*(1.+ss))
       else
          s1 = s1-(yp(i)+yp(i-1))*dels*dels*dels/24.
       end if
    end do
4   s2 = 0.
    if (x(ilm1) .ge. xxl .or. (ideltp .eq. 0 .and. .not. bdy)) go to 6
!
! integrate from x(ilm1) to xxl, store in s2
!
    del1 = xxl-x(ilm1)
    del2 = xil-xxl
    dels = xil-x(ilm1)
    t1 = del1*del1/(2.*dels)
    t2 = (del2+dels)*del1/(2.*dels)
    s2 = t1*y(il)+t2*y(ilm1)
    if (is_not_zero(sigma)) then
       c1 = coshmm_fun(sigmap * del1)
       c2 = coshmm_fun(sigmap * del2)
       ss = sinhm_fun(sigmap * dels)
       cs = coshmm_fun(sigmap * dels)
       s2 = s2+(yp(il)*del1*del1*(c1-ss/2.)+yp(ilm1)* &
            (dels*dels*(cs-ss/2.)-del2*del2*(c2-ss/2.))) &
            /(sigmap*sigmap*dels*(1.+ss))
    else
       s2 = s2-t1*(del2*(del1+dels) &
            +dels*dels)*yp(il)/12. &
            -t2*t2*dels*yp(ilm1)/6.
    end if

6   s3 = 0.
    if (xxl .ge. xil .or. (ideltp .eq. 0 .and. bdy).or. ilm1 .eq. ium1) go to 8


!
! integrate from xxl to xil, store in s3
!
    del1 = xxl-x(ilm1)
    del2 = xil-xxl
    dels = xil-x(ilm1)
    t1 = (del1+dels)*del2/(2.*dels)
    t2 = del2*del2/(2.*dels)
    s3 = t1*y(il)+t2*y(ilm1)
    if (is_not_zero(sigma)) then
       c1 = coshmm_fun(sigmap * del1)
       c2 = coshmm_fun(sigmap * del2)
       ss = sinhm_fun(sigmap * dels)
       cs = coshmm_fun(sigmap * dels)
       s3 = s3+((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
       s3 = s3-t1*t1*dels*yp(il)/6. &
            -t2*(del1*(del2+dels)+dels*dels)* &
            yp(ilm1)/12.
    end if
8   s4 = 0.
    if (ilm1 .ge. ium1-1 .or. (ideltp .eq. 0 .and. bdy)) go to 11
!
! integrate from xil to x(ium1), store in s4
!
    ilp1 = il+1
    do i = ilp1,ium1
       dels = x(i)-x(i-1)
       s4 = s4+(y(i)+y(i-1))*dels/2.
       if (is_not_zero(sigma)) then
          ss = sinhm_fun(sigmap * dels)
          cs = coshmm_fun(sigmap * dels)
          s4 = s4+(yp(i)+yp(i-1))*dels*(cs-ss/2.)/(sigmap*sigmap*(1.+ss))
       else
          s4 = s4-(yp(i)+yp(i-1))*dels*dels*dels/24.
       end if
    end do
11  s5 = 0.
    if (x(ium1) .ge. xxu .or. (ideltp .eq. 0 .and. bdy) .or. ilm1 .eq. ium1) go to 13
!
! integrate from x(ium1) to xxu, store in s5
!
    del1 = xxu-x(ium1)
    del2 = xiu-xxu
    dels = xiu-x(ium1)
    t1 = del1*del1/(2.*dels)
    t2 = (del2+dels)*del1/(2.*dels)
    s5 = t1*y(iu)+t2*y(ium1)
    if (is_not_zero(sigma)) then
       c1 = coshmm_fun(sigmap * del1)
       c2 = coshmm_fun(sigmap * del2)
       ss = sinhm_fun(sigmap * dels)
       cs = coshmm_fun(sigmap * dels)
       s5 = s5+(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
       s5 = s5-t1*(del2*(del1+dels)+dels*dels)*yp(iu)/12.-t2*t2*dels*yp(ium1)/6.
    end if
13  s6 = 0.
    if (xxu .ge. xiu .or. (ideltp .eq. 0 .and. .not. bdy)) go to 15
!
! integrate from xxu to xiu, store in s6
!
    del1 = xxu-x(ium1)
    del2 = xiu-xxu
    dels = xiu-x(ium1)
    t1 = (del1+dels)*del2/(2.*dels)
    t2 = del2*del2/(2.*dels)
    s6 = t1*y(iu)+t2*y(ium1)
    if (is_not_zero(sigma)) then
       c1 = coshmm_fun(sigmap * del1)
       c2 = coshmm_fun(sigmap * del2)
       ss = sinhm_fun(sigmap * dels)
       cs = coshmm_fun(sigmap * dels)
       s6 = s6+((dels*dels*(cs-ss/2.)-del1*del1*(c1-ss/2.)) &
            *yp(iu)+del2*del2*(c2-ss/2.)*yp(ium1))/ &
            (sigmap*sigmap*dels*(1.+ss))
    else
       s6 = s6-t1*t1*dels*yp(iu)/6.-t2*(del1*(del2+dels)+dels*dels)*yp(ium1)/12.
    end if
15  s7 = 0.
    if (iu .eq. 1 .or. (ideltp .eq. 0 .and. .not. bdy)) go to 18
!
! integrate from xiu to x1pp, store in s7
!
    np1 = n+1
    iup1 = iu+1
    do ii = iup1,np1
       im1 = ii-1
       i = ii
       if (i .eq. np1) i=1
       dels = x(i)-x(im1)
       if (dels .le. 0.) dels=dels+p
       s7 = s7+(y(i)+y(im1))*dels/2.
       if (is_not_zero(sigma)) then
          ss = sinhm_fun(sigmap * dels)
          cs = coshmm_fun(sigmap * dels)
          s7 = s7+(yp(i)+yp(im1))*dels*(cs-ss/2.)/(sigmap*sigmap*(1.+ss))
       else
          s7 = s7-(yp(i)+yp(im1))*dels*dels*dels/24.
       end if
    end do
18  s8 = 0.
    if (ilm1 .lt. ium1 .or. (ideltp .eq. 0 .and. bdy))go to 20

!
! integrate from xxl to xxu, store in s8
!
    delu1 = xxu-x(ium1)
    delu2 = xiu-xxu
    dell1 = xxl-x(ium1)
    dell2 = xiu-xxl
    dels = xiu-x(ium1)
    deli = xxu-xxl
    t1 = (delu1+dell1)*deli/(2.*dels)
    t2 = (delu2+dell2)*deli/(2.*dels)
    s8 = 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)
       s8 = s8+(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
       s8 = s8-t1*(delu2*(dels+delu1) &
            +dell2*(dels+dell1))*yp(iu)/12. &
            -t2*(dell1*(dels+dell2) &
            +delu1*(dels+delu2))*yp(ium1)/12.
    end if
20  so = s1+s2+s6+s7
    si = s3+s4+s5+s8
    if (bdy) then
      fitp_curvpi = ideltp * (so+si) + isign * so
    else
      fitp_curvpi = ideltp * (so+si) + isign * si
    end if
  end function fitp_curvpi