      subroutine gauss(ndim, ncomp, npoints, integrand, integral, fail)
      implicit none
      integer n,ndim, npoints,i1,i2,i3,i4,i5,i6,ncomp,fail
      double precision integral(*),x(npoints),w(npoints)
      double precision xx(ndim), ff(ncomp)


      fail=0

      call dsetgauss(0.d0,1.d0,npoints,x,w,fail)

      if(ndim.eq.1)then

         do n=1,ncomp
            integral(n)=0.d0
            do i1=1,npoints
               xx(1)=x(i1)
               call integrand(ndim, xx, ncomp, ff)
               integral(n)=integral(n) + ff(n)*w(i1)
            enddo
         enddo
         return

      else if(ndim.eq.2)then
         
         do n=1,ncomp
            integral(n)=0.d0
            do i1=1,npoints
               xx(1)=x(i1)
               do i2=1,npoints
                  xx(2)=x(i2)
                  call integrand(ndim, xx, ncomp, ff)
                  integral(n)=integral(n) + ff(n)*w(i1)*w(i2)
               enddo
            enddo
         enddo
         return


      else if(ndim.eq.3)then
         
         do n=1,ncomp
            integral(n)=0.d0
            do i1=1,npoints
               xx(1)=x(i1)
               do i2=1,npoints
                  xx(2)=x(i2)
                  do i3=1,npoints
                     xx(3)=x(i3)
                     call integrand(ndim, xx, ncomp, ff)
                     integral(n)=integral(n) + ff(n)*w(i1)*w(i2)*w(i3)
                  enddo
               enddo
            enddo
         enddo
         return

      else if(ndim.eq.4)then
         
         do n=1,ncomp
            integral(n)=0.d0
            do i1=1,npoints
               xx(1)=x(i1)
               do i2=1,npoints
                  xx(2)=x(i2)
                  do i3=1,npoints
                     xx(3)=x(i3)
                     do i4=1,npoints
                        xx(4)=x(i4)
                        call integrand(ndim, xx, ncomp, ff)
                        integral(n)=integral(n) 
     &                       + ff(n)*w(i1)*w(i2)*w(i3)*w(i4)
                     enddo
                  enddo
               enddo
            enddo
         enddo
         return

      else if(ndim.eq.5)then
         
         do n=1,ncomp
            integral(n)=0.d0
            do i1=1,npoints
               xx(1)=x(i1)
               do i2=1,npoints
                  xx(2)=x(i2)
                  do i3=1,npoints
                     xx(3)=x(i3)
                     do i4=1,npoints
                        xx(4)=x(i4)
                        do i5=1,npoints
                           xx(5)=x(i5)
                           call integrand(ndim, xx, ncomp, ff)
                           integral(n)=integral(n) 
     &                          + ff(n)*w(i1)*w(i2)*w(i3)*w(i4)*w(i5)
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         enddo
         return

      else if(ndim.eq.6)then
         
         do n=1,ncomp
            integral(n)=0.d0
            do i1=1,npoints
               xx(1)=x(i1)
               do i2=1,npoints
                  xx(2)=x(i2)
                  do i3=1,npoints
                     xx(3)=x(i3)
                     do i4=1,npoints
                        xx(4)=x(i4)
                        do i5=1,npoints
                           xx(5)=x(i5)
                           do i6=1,npoints
                              xx(6)=x(i6)
                              call integrand(ndim, xx, ncomp, ff)
                              integral(n)=integral(n) 
     &                             + ff(n)*w(i1)*w(i2)*w(i3)*w(i4)
     &                             *w(i5)*w(i6)
                           enddo
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         enddo
         return

      else
         fail=1
         return
      endif


      end



c
c       This subroutine stores the values for
c       the gaussian integration
c       Only works with n=4,6,8,10,16,20,24,32,64,128
c
      subroutine dsetgauss(a,b,n,x,w,fail)
      implicit none
      integer i,n,nm,kiplus,kiminus,fail
      real*8 a,b,gauss4(4),gauss6(6),gauss8(8),work(256)
      real*8 gauss10(10),gauss16(16),gauss20(20)
      real*8 gauss24(24),gauss32(32),gauss64(64),gauss128(128)
      real*8 x(128),w(128),xuni(128),wuni(128)
      data gauss4/0.339981043584856d0,0.652145154862546d0,  
     #0.861136311594053d0,0.347854845137454d0/
      data gauss6/0.238619186083197d0,0.467913934572691,
     #0.661209386466265d0,0.360761573048139,
     #0.932469514203152d0,0.171324492379170d0/
      data gauss8/0.183434642495650d0,0.362683783378362d0,
     #0.525532409916329d0,0.313706645877887d0,
     #0.796666477413627d0,0.222381034453374d0,
     #0.960289856497536d0,0.101228536290376d0/
      data gauss10/0.148874338981631d0,0.295524224714753d0,
     #0.433395394129247d0,0.269266719309996d0,
     #0.679409568299024d0,0.219086362515982d0,
     #0.865063366688985d0,0.149451349150581d0,
     #0.973906528517172d0,0.066671344308688d0/
      data gauss16/0.095012509837637d0,0.189450610455068d0,
     #0.281603550779258d0,0.182603415044923d0,
     #0.458016777657227d0,0.169156519395002d0,
     #0.617876244402643d0,0.149595988816576d0,
     #0.755404408355003d0,0.124628971255533d0,
     #0.865631202387831d0,0.095158511682492d0,
     #0.944575023073232d0,0.062253523938647d0,
     #0.989400934991649d0,0.027152459411754d0/
      data gauss20/0.076526521133497d0,0.152753387130725d0,
     #0.227785851141645d0,0.149172986472603d0,
     #0.373706088715419d0,0.142096109318382d0,
     #0.510867001950827d0,0.131688638449176d0,
     #0.636053680726515d0,0.118194531961518d0,
     #0.746331906460150d0,0.101930119817240d0,
     #0.839116971822218d0,0.083276741576704d0,
     #0.912234428251325d0,0.062672048334109d0,
     #0.963971927277913d0,0.040601429800386d0,
     #0.993128599185094d0,0.017614007139152d0/
      data gauss24/0.064056892862605d0,0.127938195346752d0,
     #0.191118867473616d0,0.125837456346828d0,
     #0.315042679696163d0,0.121670472927803d0,
     #0.433793507626045d0,0.115505668053725d0,
     #0.545421471388839d0,0.107444270115965d0,
     #0.648093651936975d0,0.097618652104113d0,
     #0.740124191578554d0,0.086190161531953d0,
     #0.820001985973902d0,0.073346481411080d0,
     #0.886415527004401d0,0.059298584915436d0,
     #0.938274552002732d0,0.044277438817419d0,
     #0.974728555971309d0,0.028531388628933d0,
     #0.995187219997021d0,0.012341229799987d0/
      data gauss32/0.048307665687738d0,0.096540088514727d0,
     #0.144471961582796d0,0.095638720079274d0,
     #0.239287362252137d0,0.093844399080804d0,
     #0.331868602282127d0,0.091173878695763d0,
     #0.421351276130635d0,0.087652093004403d0,
     #0.506899908932229d0,0.083311924226946d0,
     #0.587715757240762d0,0.078193895787070d0,
     #0.663044266930215d0,0.072345794108848d0,
     #0.732182118740289d0,0.065822222776361d0,
     #0.794483795967942d0,0.058684093478535d0,
     #0.849367613732569d0,0.050998059262376d0,
     #0.896321155766052d0,0.042835898022226d0,
     #0.934906075937739d0,0.034273862913021d0,
     #0.964762255587506d0,0.025392065309262d0,
     #0.985611511545268d0,0.016274394730905d0,
     #0.997263861849481d0,0.007018610009470d0/
      data gauss64/0.02435029266342d0,0.04869095700914d0,
     &0.07299312178780d0,0.04857546744150d0,
     &0.12146281929612d0,0.04834476223477d0,
     &0.16964442042399d0,0.04799938859631d0,
     &0.21742364374001d0,0.04754016571443d0,
     &0.26468716220877d0,0.04696818281529d0,
     &0.31132287199021d0,0.04628479657947d0,
     &0.35722015833767d0,0.04549162792405d0,
     &0.40227015796399d0,0.04459055815801d0,
     &0.44636601725346d0,0.04358372451999d0,
     &0.48940314570705d0,0.04247351510908d0,
     &0.53127946401989d0,0.04126256322055d0,
     &0.57189564620263d0,0.03995374110007d0,
     &0.61115535517239d0,0.03855015317862d0,
     &0.64896547125466d0,0.03705512854024d0,
     &0.68523631305423d0,0.03547221325688d0,
     &0.71988185017161d0,0.03380516183714d0,
     &0.75281990726053d0,0.03205792835485d0,
     &0.78397235894334d0,0.03023465707240d0,
     &0.81326531512280d0,0.02833967261426d0,
     &0.84062929625258d0,0.02637746971505d0,
     &0.86599939815409d0,0.02435270256871d0,
     &0.88931544599511d0,0.02227017380838d0,
     &0.91052213707850d0,0.02013482315353d0,
     &0.92956917213194d0,0.01795171577570d0,
     &0.94641137485840d0,0.01572603047602d0,
     &0.96100879965205d0,0.01346304789672d0,
     &0.97332682778991d0,0.01116813946013d0,
     &0.98333625388463d0,0.00884675982634d0,
     &0.99101337147674d0,0.00650445796881d0,
     &0.99634011677196d0,0.00414703325659d0,
     &0.99930504173577d0,0.00178327963660d0/
      data gauss128/0.01222369896062d0,0.02444618019626d0,
     &0.03666379096873d0,0.02443156909785d0,
     &0.06108196960414d0,0.02440235563385d0,
     &0.08546364050452d0,0.02435855726469d0,
     &0.10979423112764d0,0.02430020016797d0,
     &0.13405919946119d0,0.02422731922281d0,
     &0.15824404271422d0,0.02413995798902d0,
     &0.18233430598534d0,0.02403816868102d0,
     &0.20631559090208d0,0.02392201213669d0,
     &0.23017356422666d0,0.02379155778099d0,
     &0.25389396642269d0,0.02364688358442d0,
     &0.27746262017790d0,0.02348807601650d0,
     &0.30086543887768d0,0.02331522999401d0,
     &0.32408843502441d0,0.02312844882432d0,
     &0.34711772859764d0,0.02292784414359d0,
     &0.36993955534986d0,0.02271353585011d0,
     &0.39254027503327d0,0.02248565203258d0,
     &0.41490637955227d0,0.02224432889359d0,
     &0.43702450103710d0,0.02198971066819d0,
     &0.45888141983355d0,0.02172194953771d0,
     &0.48046407240417d0,0.02144120553878d0,
     &0.50175955913614d0,0.02114764646770d0,
     &0.52275515205118d0,0.02084144778010d0,
     &0.54343830241281d0,0.02052279248617d0,
     &0.56379664822662d0,0.02019187104117d0,
     &0.58381802162876d0,0.01984888123167d0,
     &0.60349045615855d0,0.01949402805730d0,
     &0.62280219391058d0,0.01912752360826d0,
     &0.64174169256231d0,0.01874958693853d0,
     &0.66029763227265d0,0.01836044393492d0,
     &0.67845892244772d0,0.01796032718214d0,
     &0.69621470836951d0,0.01754947582371d0,
     &0.71355437768359d0,0.01712813541906d0,
     &0.73046756674191d0,0.01669655779679d0,
     &0.74694416679706d0,0.01625500090410d0,
     &0.76297433004409d0,0.01580372865266d0,
     &0.77854847550641d0,0.01534301076088d0,
     &0.79365729476219d0,0.01487312259269d0,
     &0.80829175750791d0,0.01439434499294d0,
     &0.82244311695564d0,0.01390696411961d0,
     &0.83610291506091d0,0.01341127127272d0,
     &0.84926298757797d0,0.01290756272029d0,
     &0.86191546893955d0,0.01239613952122d0,
     &0.87405279695803d0,0.01187730734541d0,
     &0.88566771734540d0,0.01135137629108d0,
     &0.89675328804916d0,0.01081866069944d0,
     &0.90730288340176d0,0.01027947896690d0,
     &0.91731019808096d0,0.00973415335479d0,
     &0.92676925087895d0,0.00918300987166d0,
     &0.93567438827792d0,0.00862637779862d0,
     &0.94402028783022d0,0.00806458989049d0,
     &0.95180196134126d0,0.00749798192563d0,
     &0.95901475785370d0,0.00692689256690d0,
     &0.96565436643197d0,0.00635166316171d0,
     &0.97171681874714d0,0.00577263754287d0,
     &0.97719849146391d0,0.00519016183268d0,
     &0.98209610843572d0,0.00460458425670d0,
     &0.98640674272459d0,0.00401625498374d0,
     &0.99012781849173d0,0.00342552604091d0,
     &0.99325711290021d0,0.00283275147146d0,
     &0.99579275853498d0,0.00223828843096d0,
     &0.99773324862551d0,0.00164250301862d0,
     &0.99907745997738d0,0.00104581267832d0,
     &0.99982488794713d0,0.00044938068609d0/


      nm=n/2

      if(n.eq.4)then
         do i=1,n
            work(i)=gauss4(i)
         enddo
      else if(n.eq.6)then
         do i=1,n
            work(i)=gauss6(i)
         enddo
      else if(n.eq.8)then
         do i=1,n
            work(i)=gauss8(i)
         enddo
      else if(n.eq.10)then
         do i=1,n
            work(i)=gauss10(i)
         enddo
      else if(n.eq.16)then
         do i=1,n
            work(i)=gauss16(i)
         enddo
      else if(n.eq.20)then
         do i=1,n
            work(i)=gauss20(i)
         enddo
      else if(n.eq.24)then
         do i=1,n
            work(i)=gauss24(i)
         enddo
      else if(n.eq.32)then
         do i=1,n
            work(i)=gauss32(i)
         enddo
      else if(n.eq.64)then
         do i=1,n
            work(i)=gauss64(i)
         enddo
      else if(n.eq.128)then
         do i=1,n
            work(i)=gauss128(i)
         enddo
      else
         fail=2
         return
      endif

      

      do 10 i=1,nm
      kiplus=nm+i
      kiminus=nm+1-i
      xuni(kiplus)=work(2*i-1)
      xuni(kiminus)=-work(2*i-1)
      wuni(kiplus)=work(2*i)
      wuni(kiminus)=work(2*i)
 10   continue

      do i=1,n
      x(i)=(b-a)/2.d0*xuni(i)+(b+a)/2.
      w(i)=(b-a)/2.d0*wuni(i)
      end do

      return
      end






