fitp_ceez Subroutine

  pure subroutine fitp_ceez (del1,del2,sigma,c1,c2,c3,n)
    implicit none
    integer, intent(in) :: n
    real, intent(in) :: del1, del2, sigma
    real, intent(out) :: c1, c2, c3
!
!                                 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 determines the coefficients c1, c2, and c3
! used to determine endpoint slopes. specifically, if
! function values y1, y2, and y3 are given at points x1, x2,
! and x3, respectively, the quantity c1*y1 + c2*y2 + c3*y3
! is the value of the derivative at x1 of a spline under
! tension (with tension factor sigma) passing through the
! three points and having third derivative equal to zero at
! x1. optionally, only two values, c1 and c2 are determined.
!
! on input--
!
!   del1 is x2-x1 (.gt. 0.).
!
!   del2 is x3-x1 (.gt. 0.). if n .eq. 2, this parameter is
!   ignored.
!
!   sigma is the tension factor.
!
! and
!
!   n is a switch indicating the number of coefficients to
!   be returned. if n .eq. 2 only two coefficients are
!   returned. otherwise all three are returned.
!
! on output--
!
!   c1, c2, and c3 contain the coefficients.
!
! none of the input parameters are altered.
!
! this subroutine references package module snhcsh.
!
!-----------------------------------------------------------
    real :: delm, delp, sinhmp, denom, sinhmm, del, coshm2, coshm1

    if (n == 2) then
      ! two coefficients
      c1 = -1./del1
      c2 = -c1
    else if (is_zero(sigma)) then
      ! tension .eq. 0.
      del = del2-del1
      c1 = -(del1+del2)/(del1*del2)
      c2 = del2/(del1*del)
      c3 = -del1/(del2*del)
   else
      ! tension .ne. 0.
      coshm1 = coshm_fun(sigma * del1)
      coshm2 = coshm_fun(sigma * del2)
      delp = sigma*(del2+del1)/2.
      delm = sigma*(del2-del1)/2.
      sinhmp = sinhm_fun(delp)
      sinhmm = sinhm_fun(delm)
      denom = coshm1*(del2-del1)-2.*del1*delp*delm*(1.+sinhmp)*(1.+sinhmm)
      c1 = 2.*delp*delm*(1.+sinhmp)*(1.+sinhmm)/denom
      c2 = -coshm2/denom
      c3 = coshm1/denom
   end if
  end subroutine fitp_ceez