      subroutine exact_delta ! Bernd Berg Apr 07 2012. 
C Measurement of correlations functions.
      include '../../Libs/Fortran/implicit.sta'
      include '../../Libs/Fortran/constants.par'
      include '../../Libs/Fortran_par/lat.par'
      include '../../Libs/Fortran_par/mc.par'
      include '../../Libs/Fortran_par/common_lat.f'
      include '../../Libs/Scalar_par/common_exact.f'
      dimension ix(nd),dkx(nd),ix1(nd),ix2(nd)
      ltest=.false.
      ltest=.true.
C 
C Initialization of momentum space arrays:
C
      do id=1,nd ! Momenta stepsizes.
        dkx(id)=tpi/nla(id) 
      enddo
      write(iuo,'("nla:",4I4)') nla
      write(iuo,'("dkx:",4G14.6)') dkx
      do is=1,ms 
        call ixcor(ix,nla,is,nd)
        prop(is)=xm**2
        do id=1,nd
          xk(id,is)=ix(id)*dkx(id)
          prop(is)=prop(is)+four*sin(half*xk(id,is))**2
        enddo
        write(iuo,'("ik,xk:",I5,2X,4F9.5)') is,(xk(id,is),id=1,nd)
        prop(is)=one/prop(is)
c       write(iuo,'(I4,6X,5G14.6)') is,(xk(id,is),id=1,nd),prop(is)
      enddo
C
C Calculation of dnm:
C
      do is=1,ms 
        call ixcor(ix,nla,is,nd)
        it=ix(nd)
c       write(iuo,'("is,ix:",I5,2X,4I4)') is,ix
        xkmat(1,is)=zero
        xkinv(1,is)=zero
        ximat(1,is)=zero
        xiinv(1,is)=zero
        do ik=1,ms ! Momentum space sum.
          xdk=zero ! Vector product for exp of numerator.
          do id=1,nd
            xdk=xdk+ix(id)*xk(id,ik)
c           write(iuo,'("is,ik,id,ix,xk,xdk:",4I4,2G15.6)') 
c    &                   is,ik,id,ix(id),xk(id,ik),xdk
          enddo
c         write(iuo,'("is,ik,xdk:",2I4,F9.5)') is,ik,xdk
          xkmat(1,is)=xkmat(1,is)+cos(xdk)
          xiinv(1,is)=xiinv(1,is)+sin(xdk)
        enddo
        xkmat(1,is)=abs(xkmat(1,is)/ms)
        xkinv(1,is)=xkinv(1,is)/ms
      enddo
      write(iuo,'(7F10.6)') (xkmat(1,is2),is2=1,ms)
      if(ltest) stop "exact_delta: ltest."
C
C Building the rest of the matrices:
C
      do is2=1,ms
        call ixcor(ix2,nla,is2,nd)
        do is1=2,ms
          call ixcor(ix1,nla,is1,nd)
          do id=1,nd
            ix(id)=iabs(ix2(id)-ix1(id))
          enddo
          is=isfun(ix,nla,nd)
          if(xkinv(1,is).lt.zero) stop "xkinv.lt.zero"
          xkmat(is1,is2)=xkmat(1,is)
          xkinv(is1,is2)=xkinv(1,is)
          ximat(is1,is2)=ximat(1,is)
          xiinv(is1,is2)=xiinv(1,is)
        enddo
      enddo
      write(iuo,'(/,"Knm matrix Rothe p.30:",/)')
      do is1=1,ms
        write(iuo,'(7F10.6)') (xkmat(is1,is2),is2=1,ms)
      enddo
      write(iuo,'(/,"Inm matrix Rothe p.30:",/)')
      do is1=1,ms
        write(iuo,'(7F10.6)') (ximat(is1,is2),is2=1,ms)
      enddo
      write(iuo,'(/,"Inverse Knm matrix:",/)')
      do is1=1,ms
        write(iuo,'(7F10.6)') (xkinv(is1,is2),is2=1,ms)
      enddo
      write(iuo,'(/,"Inverse Inm matrix:",/)')
      do is1=1,ms
        write(iuo,'(7F10.6)') (xiinv(is1,is2),is2=1,ms)
      enddo
C
C Multiplying matrix with inverse matrix:
C
      do is2=1,ms
        do is1=1,ms
          XI(is1,is2)=zero
          do is=1,ms
            XI(is1,is2)=XI(is1,is2)+xkmat(is1,is)*xkinv(is,is2)
          enddo
          XI(is1,is2)=ABS(XI(is1,is2))
        enddo
      enddo
      write(iuo,'(/,"Test: I=Knm*Knm^{-1} matrix:",/)')
      do is1=1,ms
        write(iuo,'(7F10.6)') (XI(is1,is2),is2=1,ms)
      enddo
C
      write(iuo,'(/,"Inverse Knm matrix:",/)')
      call mat_gau(xkmat,ms,ximat,izero)
      do is1=1,ms
        write(iuo,'(7F10.6)') (xkmat(is1,is2),is2=1,ms)
      enddo
C
      return
      end
