c MAIN: test barradmap
      implicit none
c variables
      real * 8 pi
      parameter (pi=3.141592653589793d0)
      integer m
      parameter (m=10)
c common
      real * 8 xmass(m),q0
      integer n
      common/phspc/xmass,q0,n
      real * 8 barredk(4,m),xk(4,m+1),csi,y,phi,jac,vtot(4),
     # mrec2,csimax,avgi,sd,chi2
      integer i,k,l
c functions
      real * 8 random,fun1,fun2
      external random,fun1,fun2

c First test: set up random barred momenta, csi,y,phi,
c call barradmap, check that the n+1 momenta satisfy momentum conservation
      n=7
c do it few times
      do l=1,2

         do i=1,3
            barredk(i,n)=0
         enddo
         do k=1,n-1
c     possible masses
            do i=1,3
               barredk(i,k)=random()
               barredk(i,n)=barredk(i,n)-barredk(i,k)
            enddo
            barredk(4,k)=sqrt(barredk(1,k)**2+barredk(2,k)**2
     #           +barredk(3,k)**2+random())
         enddo
c     massless
         barredk(4,n)=sqrt(barredk(1,n)**2+barredk(2,n)**2
     #        +barredk(3,n)**2)
         q0=0
         do k=1,n
            q0=q0+barredk(4,k)
         enddo
c     mrec2
         mrec2=(q0-barredk(4,n))**2
     #        -barredk(1,n)**2-barredk(2,n)**2-barredk(3,n)**2
         csimax=1-mrec2/q0**2
         csi=random()*csimax
         y=1-random()
c     2*random()
         phi=2*pi*random()
         call barradmap(n,q0,barredk,csi,y,phi,xk,jac)
c     test momentum conservation
         do i=1,4
            vtot(i)=0
         enddo
         do k=1,n+1
            do i=1,4
               vtot(i)=vtot(i)+xk(i,k)
            enddo
         enddo
         write(*,*) (vtot(i),i=1,3), vtot(4)-q0
      enddo

c Second test: integrate a function in the n+1 phase space
c using standard phase space measure, and standard n-phase
c space times the radiation variables times the jacobian.
c Must get the same answer
      n=3
      q0=0
      do l=1,10
         do i=1,n-1
            xmass(i)=random()
            q0=q0+xmass(i)
         enddo
         q0=q0+random()
         xmass(n)=0
         xmass(n+1)=0
         q0=q0+n*random()
         call setexpons
         call setupvegas(3*(n+1)-4)
         call vegas(fun1,avgi,sd,chi2)
         write(*,*) avgi,sd,chi2
         call vegas(fun2,avgi,sd,chi2)
         write(*,*) avgi,sd,chi2
      enddo
      end

      subroutine setupvegas(m)
      implicit none
      integer m
      real * 8 avgi,sd,chi2a
      real * 8 xl,xu,acc
      integer ndim,ncall,itmx,nprn
      common/bveg1/xl(10),xu(10),acc,ndim,ncall,itmx,nprn
      integer j
      do j=1,10
         xl(j)=0
         xu(j)=1
      enddo
      nprn=0
      ndim=m
      itmx=10
      ncall=1000000
      end

      function fun1(x,weight)
      implicit none
      real * 8 fun1,x(10),weight
      integer m
      parameter (m=10)
c common
      real * 8 xmass(m),q0
      integer n
      common/phspc/xmass,q0,n
      real * 8 xk(4,m+1),wt0
c      functions 
      real * 8 fff
      external fff
      call bambo(n+1,q0,xmass,x,xk,wt0)
      fun1=wt0*fff(xk,n+1)
      end

      function fun2(x,weight)
      implicit none
      real * 8 fun2,x(10),weight
      real * 8 pi
      parameter (pi=3.141592653589793d0)
      integer m
      parameter (m=10)
c common
      real * 8 xmass(m),q0
      integer n
      common/phspc/xmass,q0,n
      real * 8 barredk(4,m),xk(4,m+1),wt0,mrec2,csimax,csi,y,phi,jac
c      functions 
      real * 8 fff
      external fff
      call bambo(n,q0,xmass,x(4),barredk,wt0)
c mrec2
      mrec2=(q0-barredk(4,n))**2
     #     -barredk(1,n)**2-barredk(2,n)**2-barredk(3,n)**2
      csimax=1-mrec2/q0**2
      csi=x(1)*csimax
      y=1-2*x(2)
      phi=2*pi*x(3)
      call barradmap(n,q0,barredk,csi,y,phi,xk,jac)
      fun2=wt0*jac*csimax*2*2*pi
      fun2=fun2*fff(xk,n+1)
      end

      subroutine setexpons
      implicit none
      real * 8 pw(4,10)
      common/fffc/pw
      integer i,j
      real * 8 random
      external random
      do i=1,4
         do j=1,10
            pw(i,j)=random()
         enddo
      enddo
      end

      function fff(xk,n)
      implicit none
      integer n
      real * 8 fff,xk(4,n)
      real * 8 pw(4,10)
      common/fffc/pw
      integer i,j
      fff=1
      do i=1,4
         do j=1,n
            fff=fff*abs(xk(i,j))**pw(i,j)
         enddo
      enddo
      end