function dotprod(p,q)
      implicit none
      real * 8 dotprod,p(4),q(4)
      dotprod=p(4)*q(4)-p(1)*q(1)-p(2)*q(2)-p(3)*q(3)
      end

      function random()
      real * 8 random
      call rm48(random,1)
      end

      subroutine rn3vec(vec,r)
c Generates a 3d vector in unit sphere
      implicit none
      real * 8 vec(3),r,r02,norm
      integer j
c generate in unit cube -1<x,y,z<1
 1    call rm48(vec,3)
      do j=1,3
         vec(j)=1-2*vec(j)
      enddo
c generate in unit sphere by hit and miss
      r02=vec(1)**2+vec(2)**2+vec(3)**2
      if(r02.gt.1) goto 1
c normalize
      norm=r/sqrt(r02)
      do j=1,3
         vec(j)=vec(j)*norm
      enddo
      end