;+
; NAME:
;        PDF2EXP
;
; PURPOSE:
;        This function determines the probability density function of the
;        waiting times between successive events from source 1 in the
;        presence of possible events from source 2 arriving between source 1
;        events. Both sources are assumed Poisson with a respective
;        non-extended dead time.
;
; CATEGORY:
;        Math.
;
; CALLING SEQUENCE:
;
;        Result = PDF2EXP( Wt, Mean1, Mean2, Offset1, Offset2 )
;
; INPUTS:
;        Wt:       An array of waiting times or the maximum waiting time to
;                  determine the probability density function at.
;
;        Mean1:    The average waiting time for source 1 with no dead time.
;
;        Mean2:    The average waiting time for source 2 with no dead time.
;
;        Offset1:  The dead time of source 1.
;
;        Offset2:  The dead time of source 2.
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
;
;        BIN:      Time resolution of the waiting times.  This keyword must be
;                  set if the Wt parameter is a scalar, (1=Default).
;
;        SUM:      Set this keyword to return the total probability density
;                  function, see OUTPUTS below, (0=Default).
;
; OUTPUTS:
;        This function returns the probability densities at the waiting times
;        specified by the Wt parameter if Wt is an array.  If Wt is a scalar,
;        the waiting times chosen are [0, BIN, 2*BIN, 3*BIN,... N*BIN < Wt].
;        Normally, this function will return a 2-D array, (nWT, nQ) where
;        the first dimension corresponds to the waiting times and the second
;        dimension specifies the number of events which occurred from
;        source 2 in a given waiting time. If the SUM keyword is set, then
;        the total probability density is returned (i.e. TOTAL((nWT, nQ),2)).
;
;        The probability density function follows the formula:
;
;        psd  = SUM(f_q)
;
;        p1   = 1/Mean1
;        p2   = 1/Mean2
;        q    = Number of events from source 2 occurring in the
;               waiting time, T
;        T'   = T - (Offset1 + q*Offset2)
;
;        f(p,T) = p*exp(-p*T),         Waiting time distribution for Poisson source
;        P(n,u) = u^n*exp(-u)/n!,      Prob. of n counts for Poisson source
;                                      with mean, u.
;
;        f_q  = f(p1,T')*P(q,p2*T'),    T > Offset1 + q*Offset2
;
; EXAMPLE:
;        Given data with 8 microsec resolution, with p1 = 500 Hz, p2 = 200 Hz,
;        t1 = 20 microsec, and t2 = 1 ms, to determine the PDF from 0 to 1000
;        microsec, we do the following:
;
;        Bin  = 8e-6
;        p1   = 500.
;        p2   = 200.
;        t1   = 20e-6.
;        t2   = 1e-3
;        nbin = wtmax/Bin
;        wt   = Bin*FINDGEN(nbin)/(nbin-1.)
;        psd  = PDF2EXP( wtmax, 1/p1, 1/p2, t1, t2, BIN=Bin )
;
;        or, equivalently,
;
;        Bin  = 8e-6
;        p1   = 500.*Bin
;        p2   = 200.*Bin
;        t1   = 20e-6./Bin
;        t2   = 1e-3/Bin
;        wtmax= 1000.e-6/Bin
;        wt   = FINDGEN(wtmax+1)
;        psd  = PDF2EXP( wt, 1/p1, 1/p2, t1, t2 )
;
; MODIFICATION HISTORY:
;        Written by:    Han Wen, September 1996.
;-
function PDF2EXP, Wt, Mean1, Mean2, Offset1, Offset2, BIN=Bin, SUM=Sum

;   Check parameters

         NP   = N_PARAMS()
         if (NP ne 5) then message, 'Must be called with 5 parameters: '+ $
              'Wt, Mean1, Mean2, Offset1, Offset2'

         if (N_ELEMENTS(Bin) eq 0) then Bin=1

         nWT  = N_ELEMENTS(Wt)
         Wtmax= MAX(Wt)

         t1   = Offset1
         t2   = Offset2
         p1   = 1./Mean1
         p2   = 1./Mean2

         nmax = LONG(1.*WTmax/Bin)
         if (nmax  lt 1 ) then message, 'BIN keyword set too large.'
         if (nWT eq 1) then Pden = DBLARR( nmax+1 ) else Pden = DBLARR( nWT )
         if (WTmax lt t1) then return, Pden

;   Determine max. number of events from Source 2

         eps  = 1e-25
         qmax = LONG((WTmax - 1.*t1)/t2)
         nq   = qmax + 1
         if (KEYWORD_SET(Sum) eq 0) then Pden = Pden # REPLICATE(1.,nq)

;   Probability density of a waiting time between successive events from source 1
;   and q events from source 2 occurring between the two source 1 events.

         for  q=0L,nq-1 do begin
              WTmin     = t1 + q*t2
              if (nWT gt 1) then begin
                   n    = WHERE( (Wt ge WTmin) and (Wt le WTmax) )
                   s    = Wt(n) - WTmin
              endif else begin
                   nmin = CEIL(WTmin/Bin)
                   nbin = nmax - nmin + 1
                   n    = (DINDGEN(nbin) + nmin)
                   s    = n*Bin - WTmin
              endelse
              fNq       = p1*(p2*s)^q*exp(-(p1 + p2)*s)/GAMMA(q+1)
              i         = q*(1 - KEYWORD_SET(Sum))
              Pden(n,i) = Pden(n,i) + fNq
              if (MAX(fNq) lt eps) then goto, DONE
         endfor
DONE:    return, Pden
end