subroutine daf_deltas(x0, nx0, x1, nx1, sigma, M,  eqsp, kder, deltak)
use numeric_kinds_module

!=========================================================================================
! This subroutine returns the DAF delta matrices for the zeroth, first, and
! second derivatives. These matricies are not necessarily square, nor do the set
! of evaluation point need be the same. 
!
! It is understood that these matricies will be multiplied by vectors of length
! nx1 to perform an integration by use of the matrix multiply d(nx0,nx1)*f(nx1)
! as such, each point in the d matricies contains the appropriate weight factor
! for midpoint integration. If the points given are equally spaced, then this is
! a constant factor for the entire matrix. If not, then each column has it's own
! weight factor. How this is calculated is explained later on in the program
! comments. 
!
! INPUT:
! x0 - array of points at which the function value is desired. This is the array
!        of points defining the column indecies. 
! nx0 - number of x0 points
! x1 - array of points at whcih the function is known, or better said, this is
!        the set of points over which the identity integral is performed. 
! nx1 - number of x1 points
! sigma - standard deviation of gaussian. This is optional, but if it is 
!         given as 0.d0, then this subroutine calculates a guess for it. 
! M  - number (order) of hermite polynomials used.  This is optional, but if it is 
!         given as 0, then the subroutine choses a value. 
! eqsp - logical - True if the x1 points are equally spaced. False if not. 
! kder - integer - order of the derivative desired
!
! OUTPUT:
! deltak - contains the kder-th derivative of the delta function DAF matrix
!
! Daniel Brue 2006
!=========================================================================================


implicit none

logical :: eqsp  ! are all x spacings equal?
logical :: same  ! does x0=x1?

integer :: i, j, k, l, n
integer :: kder  ! derivative order
integer :: np    ! negative or positive, that is, +1 or -1
integer :: j0, j1
integer :: nx0, nx1
integer :: M ! Cutoff for Hermite expansion
integer :: dn

real(kind=WP_Kind) :: a, b, c, d  ! terms for calculating sigma
real(kind=WP_Kind) :: wc ! equally spaced weights
real(kind=WP_Kind) :: sigma  ! determines width of hermite polynomials
real(kind=WP_Kind) :: ksig   ! sigma**kder  -- needs only to be calulated once. 
real(kind=WP_Kind) :: y   ! y(i,j)=(x1(i)-x0(j))/sigma
real(kind=WP_Kind) :: dtmp, dtmp1
real(kind=WP_Kind) :: fact, exp2s, pi
real(kind=WP_Kind) :: x0(nx0), x1(nx1)
real(kind=WP_Kind) :: deltak(nx0,nx1)
!real(kind=WP_Kind) :: de(1604401)

real(kind=WP_Kind), allocatable :: dd(:,:)

!real(kind=WP_Kind) :: c1, c2, c3   ! terms for quadratic calculation of sigma

real(kind=WP_Kind) :: w(nx1) ! weights 
real(kind=WP_Kind),allocatable :: Q(:)  ! Coefficients of Hermite expansion
real(kind=WP_Kind),allocatable :: H(:)  ! Hermite values

allocate(dd(nx0,nx1))

!CHECK FOR CORRECT DATA
if (kder < 0) then
     stop "daf_delta: negative derivative"
endif

! INITIALIZE
pi=4.d0*atan(1.d0)
deltak=0.d0

!=========================================================================================
! DEFINE M AND ALLOCATE COEFFICIENTS
! M determins how many Hermite polynomials are in the sum
print*,M, nx0, nx1
if (M == 0) then 
   M=50
end if

if(.not.allocated(Q)) allocate(Q(0:M))
if(.not.allocated(H)) allocate(H(0:M+kder))
A=0.d0
H=0.d0

!=========================================================================================
! DETERMINE SIGMA
! sigma affects how broad the hermite polynomials are and the domain
! of compact support. This value should be affected by the domain of x1
! setting this to twice the range of x for now. 
! get sigma to the kder power
!
! Make a best guess for sigma, 
! a,b,c are empirical factors provided by study of convergence over sum of M

! If given a sigma from the calling routine, use that one, else, calculate
if (sigma == 0.d0) then 
   a=1.04893d01
   b=-3.5095d0
   c=-4.744d-1-(1.d0*M)
   d=(x1(nx1)-x1(1))/(1.d0*(nx1-1)) ! find average step size
   sigma=(-b+sqrt(b*b-4*a*c))*d/(2.d0*a)
endif
ksig=sigma**kder
print*,"SIGMA", sigma, "ksig", ksig

!=========================================================================================
! DETERMINE EXPANSION COEFFICIENTS
! only the even numbered coefficients contribute, the odd are all zero
! first define Q's to be the hermite polynomials evaluated at x=zero
Q=0.d0
Q(0)=1.d0
Q(2)=-2.d0
do i=4,M,2
   Q(i)=-2.d0*(i-1)*Q(i-2)
enddo

! now scale the coefficients according to the normalization. This is done
! in a separate step so not to effect the recursion relation in the above step.
j=1
exp2s=1.d0
fact=1.d0
do i=1,M  ! need not include zero; for i=0, all altering factors are 1
   exp2s=exp2s*2.d0  ! the 2**i factor in the norm. Must be REAL to accomodate range.
                   ! 2**53 is the highest and keep 15 digits of accuracy
   j=j*(-1)
   fact=fact*i
   if (j>0) Q(i)=Q(i)/(exp2s*fact)
enddo
! now include root(pi) and sigma factors
!Q=Q*(-1)**(kder)/(sqrt((2.d0**(1+kder)*pi))*sigma)
Q=Q/(sqrt((2.d0**(1+kder)*pi))*sigma)

if(m==160) then
do i=1,M
   write(70,*)i,Q(i)
enddo
write(70,*)'&'
endif
!=========================================================================================
! DETERMINE INTEGRATION WEIGHTS
if (eqsp) then  
        ! equally spaced weights, all spacings are the same
        wc=x1(2)-x1(1)
        w=wc
else
! Unequally spaced weights:
! for each point, the weight will be the distance between the midpoints between the 
! point before and the point after. The first weight will be twice the distance 
! between the first point and the midpoint of the first and second point, which works
! out to being the distance between the first and second point. Likewise the weight
! for the last point will be double half the distance between the last point and the
! second to last point. 
! for every other point, the weight being the sum of half the distance between
! the point and the point before and half the distance between the point and the
! point after works out to be half the distance between the point before and the
! point after. 
   do i=1,nx1
      if(i == 1) then
         w(1)=x1(2)-x1(1)
      else if (i == nx1) then
         w(nx1)=x1(nx1)-x1(nx1-1)
      else
         w(i)=(x1(i+1)-x1(i-1))/2.d0
     endif
  enddo
endif

!=========================================================================================
! CALCULATE DELTA FUNCTION MATRIX ELEMENTS

same=.true.
if (nx0==nx1) then 
   ! same number of points, matricies are square
   do i=1,nx0
      if( x0(i) /= x1(i) ) then
         same = .false.
         exit
      endif
   enddo
else
   same = .false.
endif


ksig=1.d0/(sigma**(kder+1)*Sqrt(2**(kder+1)*pi)) ! this is constant for all points and all n
do j1=1,nx1
   do j0=1,nx0
      y=(x1(j1)-x0(j0))/(sigma*sqrt(2.d0))
      H(0)=1.d0!*exp(-y*y)
      H(1)=2.d0*y*H(0)
      do n=2,M+kder
         H(n)=2.d0*y*H(n-1)-2*(n-1)*H(n-2)
      enddo
      dtmp=0.d0
      fact=1.d0  ! scale 1/factorial, include a term of 10^-100 later. 
      exp2s=1.d0
      !np=(-1)**(kder)
      np=-1
      do n=0,M/2
         if (n>0) fact=fact/n  ! should later implement a better factorial
         np=np*-1
         if (n>0) exp2s=exp2s/4.d0
         dtmp=dtmp+exp2s*np*fact*H(2*n+kder)!*ksig
      enddo
      deltak(j0,j1)=dtmp*exp(-y*y)
   enddo
enddo
deltak=deltak*ksig*wc

!!   if (same.and.eqsp) then
!   ! IF SAME AND EQSP ARE TRUE, THEN THE RESULTING DELTA MATRIX IS A SYMMETRIC 
!   ! TOEPLITZ MATRIX, AND ONLY THE FIRST ROW NEEDS TO BE CALCULATED. USING THE
!   ! SYMMETRY OF THE MATRIX CAN SAVE A LOT OF TIME IF THE MATRIX IS LARGE. 
!       deltak(j0,j1)=deltak(1,j1-j0+1) ! all subsequent rows are the first row shifted over
!
!   ! now fill out the upper triangle if it is multidimensional
!   if (nx0>1) then 
!      do j0=2,nx0
!         do j1=j0,nx1
!            ! now fill out the upper triangle
!            ! d(3,3)=d(2,2)=d(1,1); d(1,2)=d(2,3)=d(3,4)=d(4,5) and so on
!
!            ! now fill out lower triangle
!            deltak(j1,j0)=deltak(j0,j1)     ! symmetric
!         enddo
!      enddo
!   endif
!
do i=1,nx1
if (kder==1)   write(80,'(4(F15.12, 3x))') x1(i), deltak(511,i)
enddo

! calculation complete, now return

end subroutine daf_deltas