;+
; NAME:
;        DTPOIDEV
; PURPOSE:
;        Return an integer random deviate drawn from a dead-time-distorted
;        Poisson distribution with a specified mean.  Adapted from POIDEV
;        procedure in "Numerical Recipes" by Press et al. (1986), Section 7.3
;        and IDL Astronomy User's library routine of the same name by
;        Wayne Landsman.
;
; CALLING SEQUENCE:
;        result = DTPOIDEV( Xm, Npts, Tinterval, Deadtime, [ SEED = ] )
;
; INPUTS:
;        Xm:       Mean of the ORIGINAL (undistorted) Poisson distribution.
;
;        Npts:     Number of deviates to generate.
;
;        Tinterval:The length of the measuring time bin.
;
;        Deadtime: The length of the NON-EXTENDED dead time in the SAME units
;                  as Tinterval.
;
; OUTPUT:
;        An array of deviates, [LONARR(Npts)].
;
; 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 dead-time-distorted Poisson noise to an integral image array, im
;             IDL> imnoise = POIDEV( im, 100, tbin, tau)
;
;        (2) Verify the expected mean  and sigma for an input value of 81
;             IDL> p = POIDEV( 81, 10000,1,1e-3) ;Test for 10,000 points
;             IDL> print,avg(p),sigma(p)
;        Average and sigma of the 10000 points should be close to
;        81/(1+81*1e-3)=75 and (1+81*1e-3)^(-3/2)*9=8.
;
; RESTRICTIONS:
;        Negative values in the input array will be returned as zeros.
;        Note: there are very subtle correlations between sequential time bins
;        in the presence of dead time.  The time interval between the beginning
;        of a bin and the first count is correlated to the interval between
;        the last count of the previous bin.  This correlation occurs when
;        the time interval between the last count of the previous bin and the
;        end of the time bin is less than the dead time.
;
; PROCEDURE:
;        For small values (observed counts < 20) independent exponential
;        deviates offset by a dead time are generated until their sum exceeds
;        the specfied mean, the number of events required is returned as the
;        Poisson deviate.   For large (observed counts > 20) values, uniform
;        random variates are compared with a Lorentzian distribution function.
;
; REVISION HISTORY:
;        Adapted from Landsman, POIDEV:     Han Wen, July 1996.
;        26-JUL-1996    Bugfix: Lorentzian distribution function falls below
;                       d-t-d Poisson distribution at low observed counts and
;                       low means. Reapplied technique of independent exponential
;                       deviates (with a dead time offset) for cases when the
;                       observed counts falls below 20.
;-
function DTPOIDEV, Xm, Npts, Tinterval, Deadtime, SEED = seed

         On_error,2

         Npar = N_PARAMS()
         if (Npar ne 4) then message, $
            'Must be called with 4 parameters: Mean, Npts, Tinterval, Deadtime'

         if (Xm le 0) then return, LONARR(Npts)
         if (Deadtime eq 0) then return, POIDEV( LONARR(Npts)+Xm,SEED=Seed )

         tbin      = Tinterval
         tau       = Deadtime

;   Determine the variance of the dead-time-distorted Poisson distribution
;
         x         = xm*tau/tbin
         lambda    = 1/(1+x)
         xobs      = xm*lambda
         sig2      = lambda^3*( xm+0.16667*lambda*x^2*( 6+4*x+x^2 ) )

;   Apply exponential deviates technique for small observed counts (<20)

         if (xobs le 20) then begin

              em        = LONARR(Npts)
              tlen      = Npts * float(tbin)
              tavg      = tbin/float(xm)
              dts       = EXPDEV(Seed, Ndts, MEAN=tavg, OFFSET=tau, SUM=tlen)
              dts(0)    = RANDOMU(Seed)*dts(0)   ; first count can have any
              tcum      = 0d0                    ; interval from left bin edge
              j         = 0L
              for i=0L,Ndts-1 do begin
                   tcum = tcum + dts(i)
                   if (tcum ge tbin) then begin  ; move j to the bin where
                        n0   = FIX(tcum/tbin)    ; the next count lies
                        tcum = tcum MOD tbin
                        j    = j+n0
                   endif
                   em(j)     = em(j)+1           ; add count to current bin
              endfor

              return,em
         endif

;   Calculate all probabilities up to 10 sigma + observed mean
;
         em_max    = CEIL(xm*lambda + 10*SQRT(sig2))
         ems       = LINDGEN(em_max+1)
         Prob      = [POIDEAD(TEMPORARY(ems),xm,tbin,tau),0]

;   Apply Lorentzian distribution comparison technique for higher observed counts (>20)
;
         sq        = SQRT( 2.*sig2 )        ;Sq and g are precomputed
         g         = Prob(long(xobs))

         Ngood     =( Ngood1 = Npts )
         good      = LINDGEN( Ngood)
         good1     = good
         y         = FLTARR(Ngood, /NOZERO ) & em = y

REJECT1: y(good)   = TAN( !PI * RANDOMU( seed, Ngood ) )
         em(good)  = sq*y(good) + xobs
         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)*Prob(fixem<(em_max+1))/g

         good2 = WHERE( RANDOMU (seed, Ngood1) GT T , Ngood)
         if ( Ngood GT 0 ) then begin
                  good1 = good1(good2)
                  good = good1
                  goto, REJECT1
         endif
         output    = long(em)

 return, output
 end