function poidev, xm, SEED = seed
  On_error,2
  compile_opt idl2

 Npts = N_elements( xm)
 
 case NPTS of 
 0: message,'ERROR - Poisson mean vector (first parameter) is undefined'
 1: output = lonarr(1) 
 else: output = make_array( SIZE = size(xm), /NOZERO ) 
 endcase 
 
   index = where( xm LE 20, Nindex, complement=big, Ncomplement=Nbig)

   if Nindex GT 0 then begin

   g = exp( -xm[ index] )           ;To compare with exponential distribution
   em1 = replicate( -1, Nindex )    ;Counts number of events
   t = replicate( 1., Nindex )          ;Counts (log) of total time

  Ngood = Nindex
  good = lindgen( Nindex)                 ;GOOD indexes the original array
  good1 = good                         ;GOOD1 indexes the GOOD vector

 REJECT:  em1[good] = em1[good] + 1      ;Increment event counter
   t = t[good1]*randomu( seed, Ngood )   ;Add exponential deviate, equivalent
                                         ;to multiplying random deviate
   good1 = where( t GT g[good], Ngood1)  ;Has sum of exponential deviates 
                                         ;exceeded specified mean?
   if ( Ngood1 GE 1 ) then begin
           good = good[ good1]
           Ngood = Ngood1
           goto, REJECT
   endif
   output[index] = em1
 endif
     if Nindex EQ Npts then return, output
; ***************************************

    xbig = xm[big]

    sq = sqrt( 2.*xbig )           ;Sq, Alxm, and g are precomputed
    alxm = alog( xbig )
    g = xbig * alxm - lngamma( xbig + 1.)

    Ngood = Nbig  & Ngood1 = Nbig
    good = lindgen( Ngood)
    good1 = good
    y = fltarr(Ngood, /NOZERO ) & em = y


REJECT1:   y[good] = tan( !PI * randomu( seed, Ngood ) )  
   em[good] = sq[good]*y[good] + xbig[good]
   good2 = where( em[good] LT 0. , Ngood )
   if (Ngood GT 0) then begin
            good = good[good2]
            goto, REJECT1
   endif

   fixem = long( em[good1] )
   test = check_math( 0, 1)         ;Don't want overflow messages
   t = 0.9*(1. + y[good1]^2)*exp( fixem*alxm[good1] - $ 
               lngamma( fixem + 1.) - g[good1] )
   good2 = where( randomu (seed, Ngood1) GT T , Ngood)
   if ( Ngood GT 0 ) then begin
            good1 = good1[good2]
            good = good1
            goto, REJECT1
   endif
   output[ big ] = long(em)

 return, output

 end