********************************************************************************
*                                                                              *
*  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 ### coulmom ###
c
c     Calculates multipole moments up to order k_max=7.
c
      subroutine coulmom (psi,f4,wgt2,dd0,dd1,dd2,dd3,dd4,
     &                              dd5,dd6,dd7,wk1,wk2)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension psi(*),f4(*),wgt2(*),wk1(*),wk2(*),
     &	dd0(nni,mxnmu),dd1(nni,mxnmu),dd2(nni,mxnmu),
     &	dd3(nni,mxnmu),dd4(nni,mxnmu),
     &  dd5(nni,mxnmu),dd6(nni,mxnmu),dd7(nni,mxnmu)
      data imulmx/8/       

      if (impole.gt.imulmx) then
         write(*,*) 'order of multipole expansion cannot exceed',
     &        imulmx
         stop 'coulmom'
      endif

c     loop over orbitals
c     rr=r  xr2=R/2  costh == xi*eta/rr

      xr2=r/2.d0	
      do iorb=1,norb
         if (itouch(iorb).eq.0) goto 500	    	
         ibeg = i1b (iorb)
         ngrid= i1si(iorb)	  
         
         do mu=1,mxnmu
            do ni=1,nni
               rr=sqrt(vxisq(mu)+vetasq(ni)-1.d0)
c              r==xrr=(R/2)rr see eq.13 (CPC 98 (1996) 346)
               xrr=xr2*rr
               if (abs(rr).lt.precis) then
	          costh=0.d0
               else
	          costh=vxi(mu)*veta(ni)/rr
               endif
c              domei=P_k (Legendre polynomial of order k)
               dome1 = costh
               dome2 =( 3.d0*costh*dome1 - 1.d0      )/ 2.d0
               dome3 =( 5.d0*costh*dome2 - 2.d0*dome1)/ 3.d0
               dome4 =( 7.d0*costh*dome3 - 3.d0*dome2)/ 4.d0
               dome5 =( 9.d0*costh*dome4 - 4.d0*dome3)/ 5.d0
               dome6 =(11.d0*costh*dome5 - 5.d0*dome4)/ 6.d0
               dome7 =(13.d0*costh*dome6 - 6.d0*dome5)/ 7.d0
               dome8 =(15.d0*costh*dome7 - 7.d0*dome6)/ 8.d0
c              dome9 =(17.d0*costh*dome8 - 8.d0*dome7)/ 9.d0
c              dome10=(19.d0*costh*dome9 - 9.d0*dome8)/10.d0

               xw=xrr
               dd0(ni,mu)=dome1*xw
               xw=xw*xrr
               dd1(ni,mu)=dome2*xw
               xw=xw*xrr
               dd2(ni,mu)=dome3*xw
               xw=xw*xrr
               dd3(ni,mu)=dome4*xw
               xw=xw*xrr
               dd4(ni,mu)=dome5*xw
               xw=xw*xrr
               dd5(ni,mu)=dome6*xw
               xw=xw*xrr
               dd6(ni,mu)=dome7*xw
               xw=xw*xrr
               dd7(ni,mu)=dome8*xw
            enddo
         enddo

         if (idbg(100).ne.0) then
            write(*,*) 'momen0: dd0 '
            call pmtx(nni,mxnmu,dd0,1,1,incrni,incrmu)
            write(*,*) 'momen0: dd2 '
            call pmtx(nni,mxnmu,dd2,1,1,incrni,incrmu)
            write(*,*) 'momen0: dd4 '
            call pmtx(nni,mxnmu,dd4,1,1,incrni,incrmu)
         endif

c        Now d0, d1, d2 and d3 is equal to P0, P1*r, P2*r**2 and P3*r**2,
c        respectively. P0, P1,..  are Legendre polynomials of a specified
c        order.

c        Prepare integrands and calculate moments. Multiplication by f4 is 
c        necessary because of the factor 2/(r*cosh(mu)) incorporated in wgt2.

         call prod2 (ngrid,psi(ibeg),psi(ibeg),wk1)
         call prod  (ngrid,f4,wk1)

         if (ihomon.ne.0) then	
            cmulti (iorb	  ) = 0.d0
         else
            call prod2 (ngrid,dd0,wk1,wk2)
            cmulti (iorb	  ) = ddot (ngrid,wgt2,1,wk2,1)
         endif
         
         call prod2 (ngrid,dd1,wk1,wk2)
         cmulti (iorb+norb) = ddot (ngrid,wgt2,1,wk2,1)

         if (ihomon.ne.0) then	
            cmulti (iorb+2*norb) = 0.d0
         else
            call prod2 (ngrid,dd2,wk1,wk2)
            cmulti (iorb+2*norb) = ddot (ngrid,wgt2,1,wk2,1)
         endif

         call prod2 (ngrid,dd3,wk1,wk2)
         cmulti (iorb+3*norb) = ddot (ngrid,wgt2,1,wk2,1)
         
         if (ihomon.ne.0) then	
            cmulti (iorb+4*norb) = 0.d0
         else
            call prod2 (ngrid,dd4,wk1,wk2)
            cmulti (iorb+4*norb) = ddot (ngrid,wgt2,1,wk2,1)
         endif

         call prod2 (ngrid,dd5,wk1,wk2)
         cmulti (iorb+5*norb) = ddot (ngrid,wgt2,1,wk2,1)

         if (ihomon.ne.0) then	
            cmulti (iorb+6*norb) = 0.d0
         else
            call prod2 (ngrid,dd6,wk1,wk2)
            cmulti (iorb+6*norb) = ddot (ngrid,wgt2,1,wk2,1)
         endif

         call prod2 (ngrid,dd7,wk1,wk2)
         cmulti (iorb+7*norb) = ddot (ngrid,wgt2,1,wk2,1)

         if (idbg(101).ne.0) then
            write(*,*)
            write(*,*) 'momen0 - multipole moments for orbital ',iorb
            write(*,1111) cmulti (iorb       ),cmulti (iorb+	norb),
     &           cmulti (iorb+2*norb),cmulti (iorb+3*norb),
     &           cmulti (iorb+4*norb),cmulti (iorb+5*norb),
     &           cmulti (iorb+6*norb),cmulti (iorb+7*norb)

 1111       format(5d13.5)
         endif
 500     continue
      enddo
      return
      end