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

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

      costh2=costh*costh
      costh3=costh2*costh
      costh4=costh2*costh2

      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
         return
      elseif (mt.eq.1) then
c        m=1
         sini2=abs(1.d0-costh2)
         sini =sqrt(sini2)
         dome(1)=sini*calp(12)
         dome(2)=costh*sini*calp(5)	  
         dome(3)=sini*(5.d0*costh2-1.d0)*calp(4)
         dome(4)=sini*(7.d0*costh3-3.d0*costh)*calp(3)
         return
      elseif (mt.eq.2) then
c        m=2
         sini2=abs(1.d0-costh2)
         sini =sqrt(sini2)
         dome(2)=sini2*calp(8)	  
         dome(3)=sini2*costh*calp(7)
         dome(4)=sini2*(7.d0*costh2-1.d0)*calp(6)
         return
      elseif (mt.eq.3) then 
c        m=3
         sini2=abs(1.d0-costh2)
         sini =sqrt(sini2)
         sini3=sini2*sini
         dome(3)=sini3*calp(10)
         dome(4)=sini3*costh*calp(9)
         return
      elseif (mt.eq.4) then
c        m=4
         sini2=abs(1.d0-costh2)
         sini =sqrt(sini2)
         sini3=sini2*sini
         sini4=sini2*sini2
         dome(4)=sini4*calp(11)
      endif
      return	
      end