C **********************************************************************
C *                                                                    *
C *                            NUMERIC (numeric.f)                     *
C *                                                                    *
C * Library module containing various numerical routines.              *
C *                                                                    *
C * Contains:                                                          *
C *   dbubble: sorts real numbers employing the bubble-sort algorithm. *
C *            lowest first                                            *
C *   zbubble: sorts complex numbers employing the bubble-sort         *
C *            algorithm. Lowest first.                                *
C *   dbubble1: same as "dbubble", except that a leading dimension can *
C *             be used. Lowest first.                                 *
C *   zbubble1: same as "zbubble", except that a leading dimension can *
C *             be used. Lowest first.                                 *
C *   abubble: sorts real numbers employing the bubble-sort algorithm. *
C *             Associated matrix to be sorted is type character       *
C *             Largest first.                                         *
C *   ibubble_up: sorts real numbers employing the bubble-sort         *
C *            algorithm lowest first. Also returns a vector of        *
C *            indices                                                 *
C *   ibubble_down: sorts real numbers employing the bubble-sort       *
C *            algorithm highest first. Also returns a vector of       *
C *            indices                                                 *
C *   dtrapez: integrates a function (given as vector) using the       *
C *            trapezoidal rule.                                       *
C *   dmatch:  finds optimal matches between two sorted lists of real  *
C *            numbers.                                                *
C *   zmatch:  finds optimal matches between two sorted lists of       *
C *            complex numbers.                                        *
C *   qqfit:     Performs a double square fit.                         *
C *   ortupr:    Sort a real array in increasing order.                *
C *   sortdowni: Sort an integer array in decreasing order.            *
C *   permutation: Generate permutations of order n.                   *
C *   ilookup:   scans a list of integers for an element.              *
C *   numder: Numerical derivatives using Ridders extrapolation        *
C *   dsinc: computes dsin(x)/x.  (MB, 8/96)                           *
C *   finddtuple: find a real*8 tuple in a list of real*8 tuples       *
C **********************************************************************


C ----------------------------------------------------------------------
C                           SUBROUTINE DBUBBLE
C
C PURPOSE:
C     Bubble-sort algorithm for double precision numbers.
C INPUT PARAMETERS:
C     vector: real vector containing the data to be sorted.
C     matrix: matrix the columns of which are to be arranged in the
C             same order as "vector" (can be a dummy array).
C     vecdim: length of "vector" and number of rows (i. e. first
C             dimension) of "matrix".
C     matdim: number of columns (i. e. second dimension) of matrix. If
C             "matrix" is a dummy array, set "matdim" to zero.
C OUTPUT PARAMETERS:
C     vector: contains the sorted data (from smallest to largest).
C     matrix: contains the "sorted" matrix.
C (See also "bsort" in propwf/getdavmat.F)
C ----------------------------------------------------------------------

      subroutine dbubble (vector, matrix, vecdim, matdim)

      implicit none

      integer vecdim,matdim,i,j,k
      real*8  vector(vecdim),matrix(vecdim,matdim),swap

C --- SORT THE DATA ---

      do i = 1,vecdim-1
         do j = i+1,vecdim
            if (vector(i) .gt. vector(j)) then
               swap = vector(i)
               vector(i) = vector(j)
               vector(j) = swap
               do k = 1,matdim
                  swap = matrix(i,k)
                  matrix(i,k) = matrix(j,k)
                  matrix(j,k) = swap
               enddo
            endif
         enddo
      enddo

      return
      end

C ----------------------------------------------------------------------
C                           SUBROUTINE ZBUBBLE
C
C PURPOSE:
C     Bubble-sort algorithm for complex numbers. The order is with
C     respect to the number's real parts.
C INPUT PARAMETERS:
C     vector: complex vector containing the data to be sorted.
C     matrix: matrix the columns of which are to be arranged in the
C             same order as "vector" (can be a dummy array).
C     vecdim: length of "vector" and number of rows (i. e. first
C             dimension) of "matrix".
C     matdim: number of columns (i. e. second dimension) of matrix. If
C             "matrix" is a dummy array, set "matdim" to zero.
C OUTPUT PARAMETERS:
C     vector: contains the sorted data (from smallest to largest).
C     matrix: contains the "sorted" matrix.
C ----------------------------------------------------------------------

      subroutine zbubble (vector, matrix, vecdim, matdim)

      implicit none

      integer    vecdim,matdim,i,j,k
      complex*16 vector(vecdim),matrix(vecdim,matdim),swap

C --- SORT THE DATA ---

      do i = 1,vecdim-1
         do j = i+1,vecdim
            if (dble(vector(i)) .gt. dble(vector(j))) then
               swap = vector(i)
               vector(i) = vector(j)
               vector(j) = swap
               do k = 1,matdim
                  swap = matrix(i,k)
                  matrix(i,k) = matrix(j,k)
                  matrix(j,k) = swap
               enddo
            endif
         enddo
      enddo

      return
      end

C ----------------------------------------------------------------------
C                           SUBROUTINE DBUBBLE1
C
C PURPOSE:
C     Bubble-sort algorithm for double precision numbers. A leading
C     dimension can be used.
C INPUT PARAMETERS:
C     vector: real vector containing the data to be sorted.
C     matrix: matrix the columns of which are to be arranged in the
C             same order as "vector" (can be a dummy array).
C     vecdim: length of "vector" and number of rows (i. e. first
C             dimension) of "matrix".
C     matdim1: leading dimension of "matrix". If "matrix" is a dummy
C              array, set "matdim" to zero.
C     matdim2: number of columns (i. e. second dimension) of matrix. If
C             "matrix" is a dummy array, set "matdim" to zero.
C OUTPUT PARAMETERS:
C     vector: contains the sorted data (from smallest to largest).
C     matrix: contains the "sorted" matrix.
C ----------------------------------------------------------------------

      subroutine dbubble1 (vector, matrix, vecdim, matdim1, matdim2)

      implicit none

      integer vecdim,matdim1,matdim2,i,j,k
      real*8  vector(vecdim),matrix(matdim1,matdim2),swap

C --- SORT THE DATA ---

      do i = 1,vecdim-1
         do j = i+1,vecdim
            if (vector(i) .gt. vector(j)) then
               swap = vector(i)
               vector(i) = vector(j)
               vector(j) = swap
               do k = 1,matdim2
                  swap = matrix(i,k)
                  matrix(i,k) = matrix(j,k)
                  matrix(j,k) = swap
               enddo
            endif
         enddo
      enddo

      return
      end

C ----------------------------------------------------------------------
C                           SUBROUTINE ZBUBBLE1
C
C PURPOSE:
C     Bubble-sort algorithm for complex numbers. A leading dimension can
C     be used. The order is with respect to the number's real parts.
C INPUT PARAMETERS:
C     vector: complex vector containing the data to be sorted.
C     matrix: matrix the columns of which are to be arranged in the
C             same order as "vector" (can be a dummy array).
C     vecdim: length of "vector" and number of rows (i. e. first
C             dimension) of "matrix".
C     matdim1: leading dimension of "matrix". If "matrix" is a dummy
C              array, set "matdim" to zero.
C     matdim2: number of columns (i. e. second dimension) of matrix. If
C             "matrix" is a dummy array, set "matdim" to zero.
C OUTPUT PARAMETERS:
C     vector: contains the sorted data (from smallest to largest).
C     matrix: contains the "sorted" matrix.
C ----------------------------------------------------------------------

      subroutine zbubble1 (vector, matrix, vecdim, matdim1, matdim2)

      implicit none

      integer    vecdim,matdim1,matdim2,i,j,k
      complex*16 vector(vecdim),matrix(matdim1,matdim2),swap

C --- SORT THE DATA ---

      do i = 1,vecdim-1
         do j = i+1,vecdim
            if (dble(vector(i)) .gt. dble(vector(j))) then
               swap = vector(i)
               vector(i) = vector(j)
               vector(j) = swap
               do k = 1,matdim2
                  swap = matrix(i,k)
                  matrix(i,k) = matrix(j,k)
                  matrix(j,k) = swap
               enddo
            endif
         enddo
      enddo

      return
      end

C ----------------------------------------------------------------------
C                           SUBROUTINE ABUBBLE
C
C PURPOSE:
C     Bubble-sort algorithm for double precision numbers.
C INPUT PARAMETERS:
C     vector: real vector containing the data to be sorted.
C     matrix: characetr matrix the columns of which are to be arranged 
C             in the same order as "vector" (can be a dummy array).
C     vecdim: length of "vector" and number of rows (i. e. first
C             dimension) of "matrix".
C     matdim: number of columns (i. e. second dimension) of matrix. If
C             "matrix" is a dummy array, set "matdim" to zero.
C OUTPUT PARAMETERS:
C     vector: contains the sorted data (from smallest to largest).
C     matrix: contains the "sorted" matrix.
C ----------------------------------------------------------------------

      subroutine abubble (vector, matrix, vecdim, matdim)

      implicit none

      integer vecdim,matdim,i,j,k
      real*8  vector(vecdim),swap
      character*(*)  matrix(vecdim,matdim),aswap*60

C --- SORT THE DATA ---

      do i = 1,vecdim-1
         do j = i+1,vecdim
            if (vector(i) .lt. vector(j)) then
               swap = vector(i)
               vector(i) = vector(j)
               vector(j) = swap
               do k = 1,matdim
                  aswap = matrix(i,k)
                  matrix(i,k) = matrix(j,k)
                  matrix(j,k) = aswap
               enddo
            endif
         enddo
      enddo

      return
      end

C ----------------------------------------------------------------------
C                           SUBROUTINE IBUBBLE_UP
C
C PURPOSE:
C     Bubble-sort algorithm for double precision numbers.
C INPUT PARAMETERS:
C     vector: real vector containing the data to be sorted.
C     ivector: integer vector 
C     vecdim: length of vectors
C
C OUTPUT PARAMETERS:
C     vector: contains the sorted data (from smallest to largest).
C     ivector: contains the indixes of the original vector
C (See also "bsort" in propwf/getdavmat.F)
C ----------------------------------------------------------------------

      subroutine ibubble_up (vector, ivector, vecdim)

      implicit none

      integer vecdim,i,j
      integer ivector(vecdim),iswap
      real*8  vector(vecdim),swap

      do i = 1,vecdim
         ivector(i)=i
      enddo

C --- SORT THE DATA ---

      do i = 1,vecdim-1
         do j = i+1,vecdim
            if (vector(i) .gt. vector(j)) then
               swap = vector(i)
               vector(i) = vector(j)
               vector(j) = swap

               iswap = ivector(i)
               ivector(i) = ivector(j)
               ivector(j) = iswap
            endif
         enddo
      enddo

      return
      end

C ----------------------------------------------------------------------
C                           SUBROUTINE IBUBBLE_DOWN
C
C PURPOSE:
C     Bubble-sort algorithm for double precision numbers.
C INPUT PARAMETERS:
C     vector: real vector containing the data to be sorted.
C     ivector: integer vector 
C     vecdim: length of vectors
C
C OUTPUT PARAMETERS:
C     vector: contains the sorted data (from largest to smallest).
C     ivector: contains the indixes of the original vector
C (See also "bsort" in propwf/getdavmat.F)
C ----------------------------------------------------------------------

      subroutine ibubble_down (vector, ivector, vecdim)

      implicit none

      integer vecdim,i,j
      integer ivector(vecdim),iswap
      real*8  vector(vecdim),swap

      do i = 1,vecdim
         ivector(i)=i
      enddo

C --- SORT THE DATA ---

      do i = 1,vecdim-1
         do j = i+1,vecdim
            if (vector(i) .lt. vector(j)) then
               swap = vector(i)
               vector(i) = vector(j)
               vector(j) = swap

               iswap = ivector(i)
               ivector(i) = ivector(j)
               ivector(j) = iswap
            endif
         enddo
      enddo

      return
      end

C ----------------------------------------------------------------------
C                          SUBROUTINE DTRAPEZ
C
C PURPOSE:
C     Evaluates an integral employing the trapezoidal rule.
C INPUT PARAMETERS:
C     func: vector containing the function values at equidistant points
C     dim: length of "func"
C     delta: distance between two consecutive grid points
C FUNCTION VALUE:
C     the computed integral
C ----------------------------------------------------------------------

      real*8 function dtrapez (func, dim, delta)

      implicit none

      integer dim,i
      real*8  func(dim),delta

C --- INTEGRATE FUNCTION ---

      dtrapez = 0.5d0*(func(1)+func(dim))
      do i = 2,dim-1
         dtrapez = dtrapez+func(i)
      enddo
      dtrapez = delta*dtrapez

      return
      end

C ----------------------------------------------------------------------
C                             SUBROUTINE DMATCH
C
C PURPOSE:
C     Finds optimal matches between two sorted lists of real numbers.
C INPUT PARAMETERS:
C     x,y:       lists (in ascending or descending order) to be matched
C     dimx,dimy: length of x and y
C OUTPUT PARAMETERS:
C     x2y:       integer array of length dimx; x2y(i)=j means that the
C                ith element of x matches best with element j of y;
C                0 means that no match has been found
C     y2x:       integer array of length dimy; y2x(j)=i means that the
C                jth element of y matches best with element i of x;
C                0 means that no match has been found
C ----------------------------------------------------------------------

      subroutine dmatch (x, y, dimx, dimy, x2y, y2x)

      implicit none

      integer dimx, dimy, x2y(1:dimx), y2x(1:dimy)
      real*8  x(1:dimx), y(1:dimy)

      logical violation
      integer posx, posy, i, j, k
      real*8  mindiff

C --- INITIALISE VARIABLES ---

      do i = 1,dimx
         x2y(i) = 0
      enddo
      do j = 1,dimy
         y2x(j) = 0
      enddo

C --- ASSIGN ALL VALUES IN X TO VALUES IN Y ---
      
      do k = 1,dimx

C    --- FIND SMALLEST DIFFERENCE BETWEEN VALUES NOT MATCHED SO FAR ---

         mindiff = abs(x(1)-y(dimy))+abs(x(dimx)-y(1))
         posx = 0
         posy = 0
         do i = 1,dimx
            if (x2y(i) .eq. 0) then
               do j = 1,dimy
                  if (y2x(j) .eq. 0) then
                     if (abs(x(i)-y(j)) .lt. mindiff) then
                        mindiff = abs(x(i)-y(j))
                        posx = i
                        posy = j
                     endif
                  endif
               enddo
            endif
         enddo
         
C    --- EXIT LOOP IF NO VALUES ARE LEFT ---

         if ((posx .eq. 0) .or. (posy .eq. 0)) goto 100

C    --- CHECK WHETHER ASSIGNMENT VIOLATES ORDER OF VALUES ---

         violation = .false.
         do i = 1,posx-1
            if (x2y(i) .ge. posy) violation = .true.
         enddo
         do j = 1,posy
            if (y2x(j) .ge. posx) violation = .true.
         enddo

C    --- SAVE MATCH OR REMOVE (X,Y)-PAIR FROM LIST ---

         if (.not. violation) then
            x2y(posx) = posy
            y2x(posy) = posx
         else
            x2y(posx) = -1
            y2x(posy) = -1
         endif
      enddo

C --- RESET ASSIGNMENTS OF VALUES NOT MATCHED ---

 100  continue
      do i = 1,dimx
         if (x2y(i) .lt. 0) x2y(i) = 0
      enddo
      do j = 1,dimy
         if (y2x(j) .lt. 0) y2x(j) = 0
      enddo

C --- RETURN TO CALLING ROUTINE ---

      return
      end

C ----------------------------------------------------------------------
C                             SUBROUTINE ZMATCH
C
C PURPOSE:
C     Finds optimal matches between two sorted lists of complex numbers.
C     The matches are with respect to the number's real parts.
C INPUT PARAMETERS:
C     x,y:       lists (in ascending or descending order) to be matched
C     dimx,dimy: length of x and y
C OUTPUT PARAMETERS:
C     x2y:       integer array of length dimx; x2y(i)=j means that the
C                ith element of x matches best with element j of y;
C                0 means that no match has been found
C     y2x:       integer array of length dimy; y2x(j)=i means that the
C                jth element of y matches best with element i of x;
C                0 means that no match has been found
C ----------------------------------------------------------------------

      subroutine zmatch (x, y, dimx, dimy, x2y, y2x)

      implicit none

      integer    dimx, dimy, x2y(1:dimx), y2x(1:dimy)
      complex*16 x(1:dimx), y(1:dimy)

      logical violation
      integer posx, posy, i, j, k
      real*8  mindiff

C --- INITIALISE VARIABLES ---

      do i = 1,dimx
         x2y(i) = 0
      enddo
      do j = 1,dimy
         y2x(j) = 0
      enddo

C --- ASSIGN ALL VALUES IN X TO VALUES IN Y ---
      
      do k = 1,dimx

C    --- FIND SMALLEST DIFFERENCE BETWEEN VALUES NOT MATCHED SO FAR ---

         mindiff = abs(dble(x(1))-dble(y(dimy)))
     +            +abs(dble(x(dimx))-dble(y(1)))
         posx = 0
         posy = 0
         do i = 1,dimx
            if (x2y(i) .eq. 0) then
               do j = 1,dimy
                  if (y2x(j) .eq. 0) then
                     if (abs(dble(x(i))-dble(y(j))) .lt. mindiff) then
                        mindiff = abs(dble(x(i))-dble(y(j)))
                        posx = i
                        posy = j
                     endif
                  endif
               enddo
            endif
         enddo
         
C    --- EXIT LOOP IF NO VALUES ARE LEFT ---

         if ((posx .eq. 0) .or. (posy .eq. 0)) goto 100

C    --- CHECK WHETHER ASSIGNMENT VIOLATES ORDER OF VALUES ---

         violation = .false.
         do i = 1,posx-1
            if (x2y(i) .ge. posy) violation = .true.
         enddo
         do j = 1,posy
            if (y2x(j) .ge. posx) violation = .true.
         enddo

C    --- SAVE MATCH OR REMOVE (X,Y)-PAIR FROM LIST ---

         if (.not. violation) then
            x2y(posx) = posy
            y2x(posy) = posx
         else
            x2y(posx) = -1
            y2x(posy) = -1
         endif
      enddo

C --- RESET ASSIGNMENTS OF VALUES NOT MATCHED ---

 100  continue
      do i = 1,dimx
         if (x2y(i) .lt. 0) x2y(i) = 0
      enddo
      do j = 1,dimy
         if (y2x(j) .lt. 0) y2x(j) = 0
      enddo

C --- RETURN TO CALLING ROUTINE ---

      return
      end

C#######################################################################
C                     QQFIT  FQFIT
C
C PURPOSE:
C     Performs a double square fit.
C     For a point x2<x<x3 one performs a parabola fit on x1,x2,x3 and
c     one on x2,x3,x4 and averrages over both fits.
C     The data in x(*) must be spaced equidistantly.
C
C     nn = number of largest grid-point below the current x-value xx.
C     xd = (xx-x(n))/dx , i.e. reduced distance of xx form grid-point.
C     x(nn-1)=x1,...,x(nn+2)=x4 ; y(nn-1)=y1,...,y(nn+2)=y4
C
C ----------------------------------------------------------------------

      real*8 function fqfit (x,y1,y2,y3)

C  Square fit.  x_input = (x-x2)/(x2-x1) 

      implicit none
      real*8    x,y1,y2,y3

      fqfit = y2 + 0.5d0*(y3-y1)*x + 0.5d0*(y3-2.d0*y2+y1)*x*x

      return
      end

C ----------------------------------------------------------------------

      subroutine qqfit(x,y,xx,v,f1,f2,n,ierr)

      implicit none

      integer   n,nn,ierr
      real*8    x(n),y(n),v,xx,dx,xd,f1,f2,fqfit

      dx = x(2)-x(1)
      nn = int( (xx-x(1)+dx)/dx )
      xd = (xx-x(1)+dx)/dx - nn

      if( nn .le. 0 )  then
         write(6,'(a,2e12.3)') 'x is outside range. x, x(1):',xx,x(1)
         ierr=1
      elseif( nn .ge. n ) then 
         write(6,'(a,2e12.3)') 'x is outside range. x, x(n):',xx,x(n)
         ierr=2
      elseif( nn .eq. 1 ) then 
         v = fqfit(xd-1.d0,y(1),y(2),y(3))
      elseif( nn .eq. n-1 ) then 
         v = fqfit(xd,y(n-2),y(n-1),y(n))
      else
         f1 = fqfit(xd,y(nn-1),y(nn),y(nn+1))
         f2 = fqfit(xd-1.d0,y(nn),y(nn+1),y(nn+2))
         v  = f1 + xd*(f2-f1)
      endif

      return
      end

C ----------------------------------------------------------------------

C#######################################################################
C Sort a real array rfunc in increasing order
C#######################################################################

      subroutine sortupr(rfunc,indx,dim)

      implicit none
      integer dim,i,j,m,indx(dim)
      real*8  rfunc(dim),rfc

      do i=1,dim-1
         do j=i+1,dim
            if(rfunc(i) .gt. rfunc(j)) then
              rfc=rfunc(j)
              rfunc(j)=rfunc(i)
              rfunc(i)=rfc
              m=indx(j)
              indx(j)=indx(i)
              indx(i)=m
            endif
         enddo
      enddo

      return
      end


C#######################################################################
C Sort an integer array ifunc in decreasing order
C#######################################################################

      subroutine sortdowni(ifunc,indx,dim)

      implicit none
      integer dim,i,j,m,indx(dim),ifunc(dim),ifc

      do i=1,dim-1
         do j=i+1,dim
            if(ifunc(i) .lt. ifunc(j)) then
              ifc=ifunc(j)
              ifunc(j)=ifunc(i)
              ifunc(i)=ifc
              m=indx(j)
              indx(j)=indx(i)
              indx(i)=m
            endif
         enddo
      enddo

      return
      end

C#######################################################################
C Generate permutations of order n. (n=input)
C nfak = n!, prmttn(j,k) = jth entry of kth permutation. wk = work
C#######################################################################

      subroutine permutation(prmttn,wk,nfak,n)

      implicit none
      integer i,j,k,ii,jj,n,nfak,nfak1,swap
      integer prmttn(n+1,*),wk(n,*)

      
      prmttn(1,1) = 1
      prmttn(2,1) = 2
      prmttn(3,1) = +1
      prmttn(1,2) = 2
      prmttn(2,2) = 1
      prmttn(3,2) = -1
      nfak=2

      if(n.le.2) return

C --- The last column, i+1, carries the signum (odd,even transpositions)

      do i = 3, n
         nfak1=nfak 
         nfak=0
C ---    Copy prmttn to wk
         do ii = 1,nfak1
            do jj = 1,i
              wk(jj,ii) = prmttn(jj,ii)
            enddo
         enddo
         do j = 1, nfak1
            nfak=nfak+1
C ---       Append a column
            do jj = 1,i-1
               prmttn(jj,nfak) = wk(jj,j)
            enddo
            prmttn(i,nfak) = i
            prmttn(i+1,nfak) = wk(i,j)   ! signum (odd/even)
C ---       transputation
            do k = i,2,-1
               nfak=nfak+1
               do jj = 1,i
                  prmttn(jj,nfak) = prmttn(jj,nfak-1)
               enddo
               swap=prmttn(k,nfak)
               prmttn(k,nfak) = prmttn(k-1,nfak)
               prmttn(k-1,nfak) = swap
               prmttn(i+1,nfak) = -prmttn(i+1,nfak-1)  ! signum
            enddo  !  k
         enddo     !  j
      enddo        !  i

c      write(6,'(a)') '---------------------------------------------'
c      do j = 1, nfak
c         write(6,'(i7,i3,a,9i3)') j, prmttn(n+1,j),
c     +                         ' |',(prmttn(i,j),i=1,n)
c      enddo
c      write(6,'(a)') '---------------------------------------------'

      return
      end

C#######################################################################
C ILOOKUP: scans a list of integers for an element (what) and
C returns the index (idx) of this element if found, 0 otherwise.
C#######################################################################

      subroutine ilookup(idx,list,n,what)

      implicit none
      integer  idx,n,i
      integer  what
      integer  list(n)

      do i=1,n
         if (list(i).eq.what) then
            idx=i
            goto 10
         endif
      enddo
      idx=0

  10  continue
      return
      end ! of ilookup


C#######################################################################

       function numder1(func,x,nvar,n,work,macheps,h)

       implicit none

      integer n,nvar,h
      real*8 numder1,x(nvar),func
      real*8 work(nvar)

      integer f
      real*8 temp,temp1,macheps,dq,vp,vm

C-----------------------------------------------------------------------
C get step value for derivatives
C-----------------------------------------------------------------------
      temp = 1/3.0d0
      temp1 = macheps**temp
      temp = temp1*x(n)
      dq = max(temp,temp1)

C-----------------------------------------------------------------------
C force step value to be accurately represented
C-----------------------------------------------------------------------
      temp = x(n) + dq
      dq = temp - x(n)

      do f=1,nvar
         work(f)  = x(f) 
      enddo
      work(n)  = x(n) + dq
      vp=func(work,h)

      work(n)  = x(n) - dq
      vm=func(work,h)

      numder1 = (vp - vm)/(2.0d0*dq)

      return
      end

C#######################################################################
C
C NUMDER1R
C
C obtain numerical first derivative of function using Ridders 
C extrapolation
C
C func: function routine
C x(nvar): coordinate 
C n: Derivative with respect to n th coordiate is formed
C h: Initial step size
C err: error
C work(nvar,2): work array
C 
C#######################################################################

      function numder1r(func,x,nvar,n,hx,err,work,h)

      integer ntab,n,nvar,h
      real*8 numder1r,err,hx,x(nvar),func,con,con2,big,safe
      real*8 work(nvar,2)
      parameter (con=1.4,con2=con*con,big=1.e30,ntab=10,safe=2.)

      integer i,j
      real*8 errt,fac,hh,a(ntab,ntab)

      numder1r=0.0d0
      if (hx .lt. 1.d-8) then
         hx = 0.001
      endif
      do i=1,nvar
         work(i,1)=x(i)
         work(i,2)=x(i)
      enddo

      hh=hx
      work(n,1)=x(n)+hh
      work(n,2)=x(n)-hh
      a(1,1)=(func(work(1,1),h)-func(work(1,2),h))
      a(1,1)=a(1,1)/(2.0*hh)
      err=big
      do i=2,ntab
         hh=hh/con
         work(n,1)=x(n)+hh
         work(n,2)=x(n)-hh
         a(1,i)=(func(work(1,1),h)-func(work(1,2),h))
         a(1,i)=a(1,i)/(2.0*hh)
         fac=con2
         do j=2,i
            a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.0)
            fac=con2*fac
            errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
            if (errt .le. err) then
               err=errt
               numder1r=a(j,i)
            endif
         enddo
         if(abs(a(i,i)-a(i-1,i-1)) .ge. safe*err) return
      enddo

      return
      end

C#######################################################################
C
C NUMDERR
C
C obtain numerical first derivative of function using Ridders 
C extrapolation
C
C func: function routine
C x: coordinate 
C hx: Initial step size
C err: error
C 
C#######################################################################

      function numderr(func,x,hx,err)

      integer ntab
      real*8 numderr,err,hx,x,func,con,con2,big,safe
      real*8 work(2)
      parameter (con=1.4,con2=con*con,big=1.e30,ntab=10,safe=2.)

      integer i,j
      real*8 errt,fac,hh,a(ntab,ntab)

      numderr=0.0d0
      if (hx .lt. 1.d-8) then
         hx = 0.001
      endif
!      do i=1,nvar
!         work(1)=x
!         work(2)=x
!      enddo

      hh=hx
      work(1)=x+hh
      work(2)=x-hh
      a(1,1)=(func(work(1))-func(work(2)))
      a(1,1)=a(1,1)/(2.0*hh)
      err=big
      do i=2,ntab
         hh=hh/con
         work(1)=x+hh
         work(2)=x-hh
         a(1,i)=(func(work(1))-func(work(2)))
         a(1,i)=a(1,i)/(2.0*hh)
         fac=con2
         do j=2,i
            a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.0)
            fac=con2*fac
            errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
            if (errt .le. err) then
               err=errt
               numderr=a(j,i)
            endif
         enddo
         if(abs(a(i,i)-a(i-1,i-1)) .ge. safe*err) return
      enddo

      return
      end

C ----------------------------------------------------------------------
C                             FUNCTION DSINC
C
C PURPOSE:
C     Computes the value of the function dsinc(x) := dsin(x)/x.
C     Special care is taken for the singularity at x = 0. All 
C     calculations are performed with double precision.
C INPUT PARAMETERS:
C     x: parameter as stated above
C FUNCTION VALUE:
C     value of the "dsinc" function as defined above
C ----------------------------------------------------------------------

      real*8 function dsinc (x)

      implicit none

      real*8 inv6, inv20, inv42, tiny
      parameter (inv6  = 1.0d0/ 6.0d0, inv20 = 1.0d0/20.0d0,
     +           inv42 = 1.0d0/42.0d0, tiny = 4.154252d-2)

C "Tiny" is chosen such that tiny**8 / 8! = eps where "eps" is the
C machine precision (here eps = 2.2 10^(-16)).

      real*8 x,x2

C --- COMPUTE FUNCTION VALUE ---

      if (dabs(x) .ge. tiny) then
         dsinc = dsin(x)/x
      else
         x2 = x**2
         dsinc = 1.0d0-inv6*x2*(1.0d0-inv20*x2*(1.0d0-inv42*x2))
      endif

      return
      end

C#######################################################################
C     find a real*8 tuple in an unsorted list of real*8 tuples
C     returns i,
C
C     input
C     N,M    first and second dimension of 'list'
C            first and only dimension of 'tuple'
C     list   array of length N of tuples of length M
C     tuple  The tuple to find in 'list'
C     tol    element wise tolerence of matching
C
C     output:
C     i <  0: not found
C     i >  0: is the i'th tuple in the list.
C#######################################################################

      subroutine finddtuple(list, tuple, M,N,tol,i)
      implicit none
      integer M,N,i
      real*8  list(M,N), tuple(M), tol

      integer k,l
      logical lmatch

      i=-1
      do k=1,N
         lmatch = .true.
         do l=1,M
            lmatch = lmatch .and. (abs(tuple(l)-list(l,k)) .le. tol)
            if (.not. lmatch) goto 221
         end do
         i = k
         goto 222
 221     continue
      end do

 222  continue

      return
      end

