subroutine daf_hermite(x0, nx0, x, nx, eqsp, f, dn, v0, v1, v2)
use numeric_kinds_module

!=========================================================================================
! This routine uses the DAF method to provide the value of a function at an arbitrary 
! point, or the functions first or second derivative. 
!
! INPUT:
! x0 - an array of length nx0 containing the values of x for which the function value
!     is desired. 
! nx0 - the size of array x0, if only one point is desired, then input as 1
! x - the array of given abscissa points for the known function values. These need not
!     be equally spaced, but it saves calculation time if they are. 
! eq - logical - true if the x points are equally spaced. False if not. 
! nx - the size of x
! f - the given function values, also an array of size nx
! dn =0 just calculate the function value
! dn = 1 calculate the first derivative and the function value
! dn = 2 calculate the first & second derivative and the function value
! 
! OUTPUT:

! v0(nx0) - the values of the function f at the points x0
! v1(nx0) - the first derivatives (0 if not requested)
! v2(nx0) - the second derivatives (0 if not requested)
!
! Daniel Brue 2006
!=========================================================================================

logical :: eqsp

integer :: i, j, k, l
integer :: nx0, nx, der
integer :: M      ! determines the number of polynomials in the sum 

real(kind=WP_Kind) :: x0(nx0), x(nx), f(nx)
real(kind=WP_Kind) :: v0(nx0) ! function values
real(kind=WP_Kind) :: v1(nx0) ! first derivative values
real(kind=WP_Kind) :: v2(nx0) ! second derivative values

real(kind=WP_Kind) :: d0(nx0,nx) ! delta function matrix
real(kind=WP_Kind) :: d1(nx0,nx) ! delta prime matrix
real(kind=WP_Kind) :: d2(nx0,nx) ! delta double prime matrix

real(kind=WP_Kind) :: center  ! used to center the distribution
real(kind=WP_Kind) :: sigma   ! determines diffuseness of Hermite Polynomials

real(kind=WP_Kind) :: w(nx)   ! integration weights

real(kind=WP_Kind) :: tmp

! CENTER THE X VALUES OF NOT ALREADY DONE
if ( abs(x(nx)) /= abs(x(1)) ) then
   ! the domain is not centered about zero. Move it for now and move it back later. 
   center=(x(nx)+x(1))/2.d0
   x=x-center
   x0=x0-center
else
   center = 0.d0
endif


! Get DAF matrices
call daf_deltas(x0, nx0, x, nx, M, sigma, eqsp, 0, d0)
v0=matmul(d0,f)
call daf_deltas(x0, nx0, x, nx, M, sigma, eqsp, 1, d0)
v1=matmul(d0,f)
call daf_deltas(x0, nx0, x, nx, M, sigma, eqsp, 2, d0)
v2=matmul(d0,f)


! MOVE THE X VALUES BACK
x=x+center
x0=x0+center


end subroutine daf_hermite