********************************************************************************
*                                                                              *
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### prepdiff ###
c
c     Initializes arrays used for differentiation.
c
      subroutine prepdiff (borb,d)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension borb(nni,mxnmu),d(nni,mxnmu)
      dimension aa1(9),aa2(9),a1(9),a2(9)
      dimension dc1(9,7,9),dc2(9,7,9)

c     derivative coeff. from 8th-order Stirling interpolation formula

      data aa1/ 3.d0, -32.d0, 168.d0, -672.d0, 0.d0, 672.d0,
     &	     -168.d0,  32.d0,	-3.d0 /
      data aa2/  -9.d0,  128.d0, -1008.d0, 8064.d0, -14350.d0,
     &	       8064.d0,-1008.d0,   128.d0,   -9.d0 /

c     initialization of the array performing differentiation
c     based on the 9-point interpolating formula

      w1=1.d0/(840.d0)
      w2=1.d0/(5040.d0)

c     initialize dni1, dni2, dmu1 and dmu2
      
      do ig=1,ngrids
         do  i=1,4
            dmu2(i,ig)=w2*aa2(i)/(hmu(ig)*hmu(ig))
            dmu1(i,ig)=w1*aa1(i)/hmu(ig)
         enddo
      enddo

      do i=1,4
         dni2(i)=w2*aa2(i)/(hni*hni)
         dni1(i)=w1*aa1(i)/hni
      enddo

c  initialize arrays used in matrix differentiation

c  initialize dmu

      ib=1
      do ig=1,ngrids
         do k=1,9
            a2(k)=w2*aa2(k)/(hmu(ig)*hmu(ig))
            a1(k)=w1*aa1(k)/hmu(ig)
         enddo 

         ie=iemu(ig)
         do i=ib,ie
            do k=1,9
               dmu(k,i)=a2(k)+borb(1,i)*a1(k)
            enddo
         enddo
         ib=ie+1
      enddo

c     initialize dni 
      
      do k=1,9
         a2(k)=w2*aa2(k)/(hni*hni)
         a1(k)=w1*aa1(k)/hni
      enddo

      do i=1,nni
         do k=1,9
            dni(k,i)=a2(k)+d(i,1)*a1(k)
         enddo
      enddo

      if (ngrids.ne.1) then
         
c        initialize arrays dc1 and dc2

c        in order to check lagrpol use two ajasent grids with equal step size 
c        and compare coefficients with those of the Stirling formulae

         call lagrpolq (dc1,dc2)

c        modify dmu array at the grid boundary region
c        derivative coeff. from 8th-order Lagrange interpolation formula
c        are stored in dc2(ngbound,imu,k) and dc1(ngbound,imu,k),
c        where ngbound is the number of 
c        the grid boundaries (1 for 1-2, 2 for 2-3 etc), imu is one of the 7 
c        intergrid mu values and k points to a derivative coefficient.

         do ig=1,ngrids-1
            ib=iemu(ig)-3
            ie=iemu(ig)+3
            imu=0
            do i=ib,ie
               imu=imu+1
               do k=1,9
c                 a2(k)=dc2(ig,imu,k)/(hmu(ig)*hmu(ig))
c                 a1(k)=dc1(ig,imu,k)/hmu(ig)
                  a2(k)=dc2(ig,imu,k)
                  a1(k)=dc1(ig,imu,k)
               enddo 
               
               do k=1,9
                  dmu(k,i)=a2(k)+borb(1,i)*a1(k)
               enddo
            enddo
         enddo
      endif

      return
      end