subroutine beamgen
      print *, 'Usage:'
      print *, '  call cohbeamgen(id,nthrow)'
      print *, '  call incohbeamgen(id,nthrow)'
      end

      subroutine cohbeamgen(id,nthrow)
      integer id,nthrow
      include 'cobrems.inc'
      logical hexist
      external hexist
      real xMinimum
      parameter (xMinimum=0.01)
      real freqMaximum
      parameter (freqMaximum=1e-4)
      real phi,theta2,theta
      real freq,f
      real rndm
      external rndm
      if (.not.hexist(id)) then
         print *, 'cohbeamgen error - id',id,' does not exist'
         return
      endif
      do i=1,nthrow
         do itry=1,1000000
            x = xMinimum**rndm(1)
            phi = dpi/4
            freq = dNcdxdp(x,phi)
            f = freq*rndm(2)
            do ip=1,q2points
               if (f.le.q2weight(ip)) then
                  theta2 = q2theta2(ip)
                  goto 5
               endif
            enddo
 5          continue
            freq = x*freq
            if (freq.gt.freqMaximum) then
               print *, 'Warning from cohbeamgen: freq=',freq,
     +                  ' is greater than freqMaximum=',freqMaximum
            endif
            if (freq.ge.freqMaximum*rndm(3)) then
               goto 50
            endif
         enddo
   50    continue
         theta=sqrt(theta2)*me/E
         call hfill(id,x*E,theta,1.0)
      enddo
      end

      subroutine incohbeamgen(id,nthrow)
      integer id,nthrow
      include 'cobrems.inc'
      logical hexist
      external hexist
      real xMinimum
      parameter (xMinimum=0.01)
      real freqMaximum
      parameter (freqMaximum=1e-2)
      real Theta02
      parameter (Theta02=1.8)
      real phi,theta2,theta
      real freq,f
      real rndm
      external rndm
      if (.not.hexist(id)) then
         print *, 'cohbeamgen error - id',id,' does not exist'
         return
      endif
      do i=1,nthrow
         do itry=1,1000000
            x = xMinimum**rndm(1)
            rand2 = rndm(2)
            theta2 = Theta02*rand2/(1-rand2+1e-30)
            freq = dNidxdt2(x,theta2)
            freq = freq*(Theta02+theta2)**2/(2*dpi*Theta02)
            freq = x*freq
            if (freq.gt.freqMaximum) then
               print *, 'Warning from incohbeamgen: freq=',freq,
     +                  ' is greater than freqMaximum=',freqMaximum
            endif
            if (freq.ge.freqMaximum*rndm(3)) then
               goto 50
            endif
         enddo
   50    continue
         theta=sqrt(theta2)*me/E
         call hfill(id,x*E,theta,1.0)
      enddo
      end