SUBROUTINE get_daf_s(delta, mval, dafx, hermit, ndaf, nd_daf, toe, usys2, kmax)
! CALL get_daf      (nmax+1, rvals, h, mval, dafx, hermit, ndaf, nmax, toeplitz, 2.d0*mu, kmax)
!-----------------------------------------------------------------------
! This routine was written by G. A. Parker
! If you find an error or have an improvement please send a messge to
! parker@phyast.nhn.uoknor.edu 
!-----------------------------------------------------------------------
!this subroutine determines the sigmas to be used throughout the daf 
!method and also determines the daf deltas for later use
!-----------------------------------------------------------------------
USE Numeric_Kinds_Module
USE CommonInfo_Module
IMPLICIT NONE
CHARACTER(len=17) :: filename
INTEGER  :: mval, jhalf, ndaf, kmax, idaf, karg, nd_daf
REAL(dp) :: dafx(0:nd_daf,0:ndaf), hermit(0:mval+ndaf+1), toe(0:nd_daf)
REAL(dp) :: delta, sigma, usys2, aterm, bterm, cterm, sigma_ratio
 
!----------------------------------------------------------------------- 
! This section gets a guess for the best sigma. One could refine
! this guess but is probably not very useful.
!-----------------------------------------------------------------------
 aterm=10.4893d0
 bterm=-3.5095d0
 cterm=-0.4744d0-mval   ! mval is a fixed parameter.
 sigma_ratio=(-bterm+sqrt(bterm**2-4.d0*aterm*cterm))/(2.d0*aterm) ! quadratic solution
 sigma=sigma_ratio*delta  ! delta is coordinate grid spacing
 jhalf=nd_daf/2+1  ! nd_daf = nmax (total dimension - 1)
 kmax=(jhalf-1)
!----------------------------------------------------------------------- 
! calculate the daf deltas.
!-----------------------------------------------------------------------
 kmax=MAX(kmax,nd_daf)
 dafx=0.d0  ! Initialize dafx
 CALL daf_s(sigma, delta, mval, dafx, kmax, hermit, ndaf, nd_daf)
 DO karg=0,MAX(kmax,nd_daf)
   IF(karg.le.kmax)THEN
      toe(karg)=-dafx(karg,2)/usys2  ! apply schrodinger kinetic coeff.
   ELSE
      toe(karg)=0.d0
   ENDIF
 ENDDO

RETURN
END SUBROUTINE get_daf_s