********************************************************************************
*                                                                              * 
*  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 ### zmulti ###	
c
      subroutine zmulti(ixa,pot)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension pot(*),dome(10)

      do j=mxnmu-3,mxnmu
         itt=(j-1)*nni
         do i=1,nni
            kk=i+itt
            rr=sqrt(vxisq(j)+vetasq(i)-1.d0)
            cost=veta(i)*vxi(j)/rr
            rr1=1.d0/(rr*r2)

            dome(1)=cost
            dome(2)=(3.d0*cost*cost-1.d0)*5.d-01
            do n=2,(impole-1)
               dome(n+1)=(dble(2*n+1)*cost*dome(n)-dble(n)*
     &              dome(n-1))/dble(n+1)
            enddo
            
            pe=0.d0
            do m=1,impole
               kxk=ixa+(m-1)*norb
               pe=pe+cmulti(kxk)*dome(m)*(rr1**dble(m+1))
            enddo 
            pot(kk)=r2*vxi(j)*(pe+rr1)	  
         enddo
      enddo
      return
      end

c ### zmultitail ###
c
      subroutine zmultitail(iorb,pot)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension pot(*),dome(10),pottmp(10)

c     potentials are calculated for ni=Pi/2

      j=mxnmu
      itt=(j-1)*nni
      i=(nni-1)/2
      kk=i+itt
      rr=sqrt(vxisq(j)+vetasq(i)-1.d0)
      costh=veta(i)*vxi(j)/rr
      xr=1.d0/(rr*r2)
      
      dome(1)=costh
      dome(2)=(3.d0*costh*costh-1.d0)*0.5d0
      do n=2,(impole-1)
         dome(n+1)=(dble(2*n+1)*costh*dome(n)
     &        -dble(n)*dome(n-1))/dble(n+1)
      enddo
      
      pe=0.d0
      xrr=xr
      do m=1,impole
         xrr=xrr*xr
         kxk=iorb+(m-1)*norb
         pe=pe+cmulti(kxk)*dome(m)*xrr
         pottmp(m)=r2*vxi(j)*(pe+xr)	  
      enddo
      
      write(*,1000) iorn(iorb),bond(iorb),gut(iorb),pottmp(1),
     &     (pottmp(m)-pottmp(m-1),m=2,impole),pottmp(impole)
 1000 format(i4,1x,a8,a1,3x,/4d13.5/5d13.5)
      return
      end