fitp_terms Subroutine

  pure subroutine fitp_terms (diag,sdiag,sigma,del)
    implicit none
    real, intent(in) :: sigma, del
    real, intent(out) :: diag, sdiag
!
!                                 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 subroutine computes the diagonal and superdiagonal
! terms of the tridiagonal linear system associated with
! spline under tension interpolation.
!
! on input--
!
!   sigma contains the tension factor.
!
! and
!
!   del contains the step size.
!
! on output--
!
!                sigma*del*cosh(sigma*del) - sinh(sigma*del)
!   diag = del*--------------------------------------------.
!                     (sigma*del)**2 * sinh(sigma*del)
!
!                   sinh(sigma*del) - sigma*del
!   sdiag = del*----------------------------------.
!                (sigma*del)**2 * sinh(sigma*del)
!
! and
!
!   sigma and del are unaltered.
!
! this subroutine references package module snhcsh.
!
!-----------------------------------------------------------
    real :: coshm, denom, sigdel, sinhm

    if (is_zero(sigma)) then
       diag = del/3.
       sdiag = del/6.
    else
       sigdel = sigma*del
       sinhm = sinhm_fun(sigdel)
       coshm = coshm_fun(sigdel)
       denom = sigma*sigdel*(1.+sinhm)
       diag = (coshm-sinhm)/denom
       sdiag = sinhm/denom
       ! Note we could use the cosh and sinh intrinsics to evalute the expression
       ! directly as below. This is found to agree very well with the existing form.
       !diag = del * (sigdel*cosh(sigdel) - sinh(sigdel))/(sigdel*sigdel*sinh(sigdel))
       !sdiag = del * (sinh(sigdel) - sigdel)/(sigdel*sigdel*sinh(sigdel))
    end if
  end subroutine fitp_terms