SUBROUTINE DAF_s(sigma, delta, mval, dafx, kmax, hermit, ndaf, nd_daf)
USE Numeric_Kinds_Module
!     CALL daf(sigma, delta, mval, dafx, kmax, hermit, ndaf, nd_daf, nd_daf)
!==============================================================================
! 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 
! ndaf is a fixed parameter
!  kdaf is order of derivative from 0 to 2
!-----------------------------------------------------------------------
!this subroutine calculates the daf deltas
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER, PARAMETER:: Msg_Unit=13
INTEGER lval, mval, kmax, mval2, kval, kmaxx, kdaf, ndaf, khalf, idaf, nd_daf
REAL(dp) :: sigma, delta, twopi, harg, expfac, fac
REAL(dp) :: dafx(0:nd_daf,0:ndaf), hermit(0:mval+ndaf+1)
REAL(dp),PARAMETER :: tol=1.d-50
!----------------------------------------------------------------------- 
! calculate the daf deltas.
!-----------------------------------------------------------------------
twopi=8.d0*atan(1.d0)
mval2=mval/2  ! mval is a fixed parameter
kmaxx=kmax    ! kmax is calculated as = jhalf-1 = nmax/2
DO kval=0,kmax 
   harg=(kval*delta)/(SQRT(2.d0)*sigma)
   expfac=EXP(-harg**2)
   CALL hermit_ex_s(harg, mval+ndaf, hermit, expfac)
   DO kdaf=0,ndaf  
      khalf=kdaf/2
      dafx(kval,kdaf)=0.d0
      fac=1.d0
      DO idaf=1,khalf
         fac=fac*idaf  ! I think this is the factorial.
      ENDDO
      DO lval=0,mval2
         IF(lval.ne.0)fac=(lval+khalf)*fac/lval
         dafx(kval,kdaf)=dafx(kval,kdaf)+hermit(2*lval+kdaf)*fac 
      ENDDO
      dafx(kval,kdaf)=(-1.d0/(SQRT(2.d0)*sigma))**kdaf*dafx(kval,kdaf)*delta/(SQRT(twopi)*sigma)*(-4.d0)**khalf
   ENDDO
   IF(ABS(dafx(kval,ndaf)).gt.tol)THEN
      kmaxx=kval
   ELSE IF(ABS(dafx(kval,ndaf)).lt.tol) THEN
      GO TO 5
   ENDIF
ENDDO
!-----------------------------------------------------------------------
!write any important information to file
!-----------------------------------------------------------------------
WRITE(Msg_Unit,*)'Program must be modified to increase dimension'
WRITE(Msg_Unit,*)'on array dafx'
WRITE(Msg_Unit,*)'kmaxx=',kmaxx
WRITE(Msg_Unit,*)'mval=',mval
WRITE(Msg_Unit,*)'1  nd_daf=',nd_daf
STOP 'DAF'

5  CONTINUE
IF(kmaxx.gt.nd_daf)THEN
   WRITE(Msg_Unit,*)'Program must be modified to increase dimension'
   WRITE(Msg_Unit,*)'on array dafx'
   WRITE(Msg_Unit,*)'kmaxx=',kmaxx
   WRITE(Msg_Unit,*)'mval=',mval
   WRITE(Msg_Unit,*)'2  nd_daf=',nd_daf
   STOP 'DAF'
ENDIF

kmax=kmaxx

RETURN
END SUBROUTINE DAF_s