      FUNCTION MCPOI(XMU) ! return random integer 0 ...
*     from Poisson distribution with mean XMU
*     (for large XMU slow,  use normal approximation) 
      DATA ARG/0.0/
      SAVE EXPM,ARG
      IF(XMU.NE.ARG) THEN
         ARG =XMU  ! mean value    
         EXPM=EXP(-ARG)
      END IF
      K=0      
      A=1.0
   10 U=RANU()   ! U(0,1) random number
      A=A*U
      IF(A.GE.EXPM) THEN
         K=K+1
         GOTO 10
      END IF 
      MCPOI=K ! K follows Poisson distribution
      END 
