********************************************************************************
*                                                                              *
*  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 ### mutli1 ###
c
      subroutine multi1(i,j,mt,ipc,pe)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension dome(10),pe(10)

      do k=1,impole
         pe(k)=0.d0
      enddo

      if (mt.gt.4) return

      rr=sqrt(vxisq(j)+vetasq(i)-1.d0)		
      costh=0.d0
      if (abs(rr).gt.precis) then
         costh=veta(i)*vxi(j)/rr	 
      endif

      rr=rr*r2
      pe(1)=0.d0

      if (mt.eq.0) then
c        m=0
         dome( 1)=costh
         dome( 2)=( 3.d0*costh*dome(1) - 1.d0)	  / 2.d0
         dome( 3)=( 5.d0*costh*dome(2) - 2.d0*dome(1))/ 3.d0
         dome( 4)=( 7.d0*costh*dome(3) - 3.d0*dome(2))/ 4.d0
         dome( 5)=( 9.d0*costh*dome(4) - 4.d0*dome(3))/ 5.d0
         dome( 6)=(11.d0*costh*dome(5) - 5.d0*dome(4))/ 6.d0
         dome( 7)=(13.d0*costh*dome(6) - 6.d0*dome(5))/ 7.d0
         dome( 8)=(15.d0*costh*dome(7) - 7.d0*dome(6))/ 8.d0
         
         xrn=rr*rr
         pe(1)=        dome(1)*excdi(ipc)/xrn
         xrn=xrn*rr
         pe(2)=pe(1) + dome(2)*excqu(ipc)/xrn
         xrn=xrn*rr
         pe(3)=pe(2) + dome(3)*excoc(ipc)/xrn
         xrn=xrn*rr
         pe(4)=pe(3) + dome(4)*exche(ipc)/xrn
         xrn=xrn*rr
         pe(5)=pe(4) + dome(5)*exc5 (ipc)/xrn
         xrn=xrn*rr
         pe(6)=pe(5) + dome(6)*exc6 (ipc)/xrn
         xrn=xrn*rr
         pe(7)=pe(6) + dome(7)*exc7 (ipc)/xrn
         xrn=xrn*rr
         pe(8)=pe(7) + dome(8)*exc8 (ipc)/xrn
         return	
c        for m not equal zero the associated Legendre polynomials 
c        are multiplied by sqrt((k-|q|)!/(k+|q|)!) factor 
      elseif (mt.eq.1) then
c        m=1
         costh2=costh*costh
         costh4=costh2*costh2
         costh6=costh4*costh2
         
         sini2=abs(1.d0-costh2)
         sini =sqrt(sini2)
         
         dome(1)=-sini/sqrt(2.d0)
         dome(2)=-3.d0*sini*costh/sqrt(6.d0)
         dome(3)=(3.d0/2.d0)*sini*(1.d0-5.d0*costh2)/sqrt(12.d0)
         dome(4)=(5.d0/2.d0)*sini*costh*(3.d0-7.d0*costh2)/sqrt(20.d0)
         dome(5)=(15.d0/8.d0)*sini
     &        *(-1.d0+14.d0*costh2-21.d0*costh4)/sqrt(30.d0)
         dome(6)=(21.d0/8.d0)*sini*costh
     &        *(-5.d0+30.d0*costh2-33.d0*costh4)/sqrt(42.d0)
         dome(7)=(7.d0/16.d0)*sini
     &        *(5.d0-135.d0*costh2+495.d0*costh4
     &        -429.d0*costh6)/sqrt(56.d0)
         dome(8)=(9.d0/16.d0)*sini*costh*(35.d0-385.d0*costh2
     &        +1001.d0*costh4-715.d0*costh6)/sqrt(72.d0)

c        minus sign for odd values of mt
c         pe=-pe

         xrn=rr*rr
         pe(1)=        dome(1)*excdi(ipc)/xrn
         xrn=xrn*rr
         pe(2)=pe(1) + dome(2)*excqu(ipc)/xrn
         xrn=xrn*rr
         pe(3)=pe(2) + dome(3)*excoc(ipc)/xrn
         xrn=xrn*rr
         pe(4)=pe(3) + dome(4)*exche(ipc)/xrn
         xrn=xrn*rr
         pe(5)=pe(4) + dome(5)*exc5 (ipc)/xrn
         xrn=xrn*rr
         pe(6)=pe(5) + dome(6)*exc6 (ipc)/xrn
         xrn=xrn*rr
         pe(7)=pe(6) + dome(7)*exc7 (ipc)/xrn
         xrn=xrn*rr
         pe(8)=pe(7) + dome(8)*exc8 (ipc)/xrn
         return
      elseif (mt.eq.2) then
c        m=2
         costh2=costh*costh
         costh4=costh2*costh2
         costh6=costh4*costh2
         costh8=costh6*costh2
         sini2=abs(1.d0-costh2)

         dome(1)=0.d0
         dome(2)=3.d0*sini2/sqrt(24.d0)
         dome(3)=15.d0*costh*sini2/sqrt(120.d0)
         dome(4)=(15.d0/2.d0)*(-1.d0+8.d0*costh2
     &        -7.d0*costh4)/sqrt(360.d0)
         dome(5)=(105.d0/2.d0)*costh*(-1.d0+4.d0*costh2
     &        -3.d0*costh4)/sqrt(840.d0)
         dome(6)=(105.d0/8.d0)
     &        *(1.d0-19.d0*costh2+51.d0*costh4
     &        -33.d0*costh6)/sqrt(1680.d0)
         dome(7)=(63.d0/8.d0)*costh
     &        *(15.d0-125.d0*costh2+253.d0*costh4
     &        -143.d0*costh6)/sqrt(3024.d0)
         dome(8)=(315.d0/8.d0)*(-1.d0+34.d0*costh2-176.d0*costh4
     &        +286.d0*costh6-143.d0*costh8)/sqrt(5040.d0)

         xrn=rr*rr
         pe(1)=        dome(1)*excdi(ipc)/xrn
         xrn=xrn*rr
         pe(2)=pe(1) + dome(2)*excqu(ipc)/xrn
         xrn=xrn*rr
         pe(3)=pe(2) + dome(3)*excoc(ipc)/xrn
         xrn=xrn*rr
         pe(4)=pe(3) + dome(4)*exche(ipc)/xrn
         xrn=xrn*rr
         pe(5)=pe(4) + dome(5)*exc5 (ipc)/xrn
         xrn=xrn*rr
         pe(6)=pe(5) + dome(6)*exc6 (ipc)/xrn
         xrn=xrn*rr
         pe(7)=pe(6) + dome(7)*exc7 (ipc)/xrn
         xrn=xrn*rr
         pe(8)=pe(7) + dome(8)*exc8 (ipc)/xrn
         return

      elseif (mt.eq.3) then
c        m=3
         costh2=costh*costh
         costh4=costh2*costh2

         sini2=abs(1.d0-costh2)
         sini =(sini2)
         sini3=sini2*sini

         dome(1)=0.d0
         dome(2)=0.d0
         dome(3)=-15.d0*sini3/sqrt(720.d0)
         dome(4)=-105.d0*sini3*costh/sqrt(5040.d0)
         dome(5)=-(105.d0/2.d0)*sini3
     &        *(-1.d0+9.d0*costh2)/sqrt(20160.d0)
         dome(6)=-(315.d0/2.d0)*sini3*costh
     &        *(-3.d0+11.d0*costh2)/sqrt(60480.d0)
         dome(7)=-(315.d0/8.d0)*sini3
     &        *(3.d0-66.d0*costh2+143.d0*costh4)/sqrt(151200.d0)
         dome(8)=-(3465.d0/8.d0)*sini3*costh
     &        *(3.d0-26.d0*costh2+39.d0*costh4)/sqrt(332640.d0)

c        minus sign for odd values of mt
c         pe=-pe

         xrn=rr*rr
         pe(1)=        dome(1)*excdi(ipc)/xrn
         xrn=xrn*rr
         pe(2)=pe(1) + dome(2)*excqu(ipc)/xrn
         xrn=xrn*rr
         pe(3)=pe(2) + dome(3)*excoc(ipc)/xrn
         xrn=xrn*rr
         pe(4)=pe(3) + dome(4)*exche(ipc)/xrn
         xrn=xrn*rr
         pe(5)=pe(4) + dome(5)*exc5 (ipc)/xrn
         xrn=xrn*rr
         pe(6)=pe(5) + dome(6)*exc6 (ipc)/xrn
         xrn=xrn*rr
         pe(7)=pe(6) + dome(7)*exc7 (ipc)/xrn
         xrn=xrn*rr
         pe(8)=pe(7) + dome(8)*exc8 (ipc)/xrn
         return

      elseif (mt.eq.4) then
c        m=4
         costh2=costh*costh
         costh4=costh2*costh2
         costh6=costh4*costh2
         costh8=costh6*costh2

         dome(1)=0.d0
         dome(2)=0.d0
         dome(3)=0.d0
         dome(4)=105.d0*(1.d0-2.d0*costh2+costh4)/sqrt(40320.d0)
         dome(5)=945.d0*costh*(1.d0-2.d0*costh2+costh4)/sqrt(362880.d0)
         dome(6)=(945.d0/2.d0)
     &        *(-1.d0+13*costh2-23.d0*costh4
     &        +11.d0*costh6)/sqrt(1814400.d0)
         dome(7)=(3465.d0/2.d0)*costh
     &        *(-3.d0+19*costh2-29.d0*costh4
     &        +13.d0*costh6)/sqrt(6652800.d0)
         dome(8)=(10395.d0/8.d0)*(1.d0-28.d0*costh2+118.d0*costh4
     &        -156.d0*costh6+65.d0*costh8)/sqrt(19958400.d0)

         xrn=rr*rr
         pe(1)=        dome(1)*excdi(ipc)/xrn
         xrn=xrn*rr
         pe(2)=pe(1) + dome(2)*excqu(ipc)/xrn
         xrn=xrn*rr
         pe(3)=pe(2) + dome(3)*excoc(ipc)/xrn
         xrn=xrn*rr
         pe(4)=pe(3) + dome(4)*exche(ipc)/xrn
         xrn=xrn*rr
         pe(5)=pe(4) + dome(5)*exc5 (ipc)/xrn
         xrn=xrn*rr
         pe(6)=pe(5) + dome(6)*exc6 (ipc)/xrn
         xrn=xrn*rr
         pe(7)=pe(6) + dome(7)*exc7 (ipc)/xrn
         xrn=xrn*rr
         pe(8)=pe(7) + dome(8)*exc8 (ipc)/xrn

      endif
      return
      end