FUNCTION g_poiss, mu
; Subroutine to return samples from poisson distribution
; If mu < 11. the correct poisson values are calculated and a random
; value chosen, otherwise a gaussian approximation is used.
; 5/27/92 dd

COMMON rand, seed

; Comment this out to use systime for seed
; if n_elements(seed) eq 0 then seed = 17
if n_elements(mu) eq 0 then mu = 1.0

; pre-calculate the factorials, 0-20
fact = [1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3628800.,$
	39916800.,479001600.,6227020800.,8.71783E10,1.30767E12,$
	2.09228E13,3.55687E14,6.40237E15,1.21645E17,2.43290E18]
probx = fact

x = float(mu)

FOR i=0,n_elements(mu)-1 DO BEGIN

  IF(mu(i) GE 11.0) THEN BEGIN
; Gaussian approximation
    x(i) = FLOAT(mu(i) + sqrt(mu(i))*g_rand())
    x(i) = FLOAT(LONG(x(i)))
  ENDIF ELSE BEGIN
; Poisson calculation
    ; Fill array of cumulative probabilities for 0 to 20 events
    probx(0) = 1.0
    FOR ip=1,20 DO BEGIN
      probx(ip) = probx(ip-1) + (mu(i)^ip)/fact(ip)
    END
    p_norm = EXP(-1.*mu(i))
    probx = p_norm*probx
    r_val = randomu(seed)
    out_val = FLOAT(indgen(23))
    x(i) = INTERPOL(out_val,[0.,probx,1.0],[r_val])
    x(i) = FLOAT(LONG(x(i)))
  ENDELSE

END

RETURN, x
END