function poidev, xm, SEED = seed
;+
; NAME:
;     POIDEV
; PURPOSE:
;     Generate a Poisson random deviate
; EXPLANATION:
;     Return an integer random deviate drawn from a Poisson distribution with
;     a specified mean.    Adapted from procedure of the same name in 
;     "Numerical Recipes" by Press et al. (1992), Section 7.3
;
;     NOTE: This routine became partially obsolete in V5.0 with the 
;     introduction of the POISSON keyword to the intrinsic functions 
;     RANDOMU and RANDOMN.     However, POIDEV is still useful for adding 
;     Poisson noise to an existing image array, for which the coding is much 
;     simpler than it would be using RANDOMU (see example 1) 
; CALLING SEQUENCE:
;     result = POIDEV( xm, [ SEED = ] )
;
; INPUTS:
;     xm - numeric scalar, vector or array, specifying the mean(s) of the 
;          Poisson distribution
;
; OUTPUT:
;     result - Long integer scalar or vector, same size as xm
;
; OPTIONAL KEYWORD INPUT-OUTPUT:
;     SEED -  Scalar to be used as the seed for the random distribution.  
;             For best results, SEED should be a large (>100) integer.
;             If SEED is undefined, then its value is taken from the system 
;             clock (see RANDOMU).    The value of SEED is always updated 
;             upon output.   This keyword can be used to have POIDEV give 
;             identical results on consecutive runs.     
;
; EXAMPLE:
;     (1) Add Poisson noise to an integral image array, im
;              IDL> imnoise = POIDEV( im)
;
;     (2) Verify the expected mean  and sigma for an input value of 81
;              IDL> p = POIDEV( intarr(10000) + 81)   ;Test for 10,000 points
;              IDL> print,mean(p),sigma(p)
;     Mean and sigma of the 10000 points should be close to 81 and 9
;
; METHOD: 
;     For small values (< 20) independent exponential deviates are generated 
;     until their sum exceeds the specified mean, the number of events 
;     required is returned as the Poisson deviate.   For large (> 20) values,
;     uniform random variates are compared with a Lorentzian distribution 
;     function.
;
; NOTES:
;     Negative values in the input array will be returned as zeros.  
;
;       
; REVISION HISTORY:
;      Version 1               Wayne Landsman        July  1992
;      Added SEED keyword                            September 1992
;      Call intrinsic LNGAMMA function               November 1994
;      Converted to IDL V5.0   W. Landsman   September 1997
;      Use COMPLEMENT keyword to WHERE()        W. Landsman August 2008
;-
  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