;+
; NAME:
;        PROB2EXP
;
; PURPOSE:
;        This function determines the probability 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 = FWT2EXP( Wtmax, 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, (1=Default).
;
;        SUM:      Set this keyword to return the total probability
;                  function, see OUTPUTS below, (0=Default).
;
; OUTPUTS:
;        This function returns the probabilities 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 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
;
;        To calculate the probability, one must take into account that events
;        can fall anywhere within the time resolution, BIN.
;
; EXAMPLE:
;        Given data with 8 microsec resolution, with p1 = 500 Hz, p2 = 200 Hz,
;        t1 = 20 microsec, and t2 = 1 ms, to determine the Probabilities from
;        0 to 1000 microsec, we do the following:
;
;        Bin  = 8e-6
;        p1   = 500.
;        p2   = 200.
;        t1   = 20e-6.
;        t2   = 1e-3
;        wtmax= 1000.e-6/Bin
;        nbin = wtmax/Bin
;        wt   = Bin*FINDGEN(nbin)/(nbin-1.)
;        psd  = PROB2EXP( wtmax, 1/p1, 1/p2, t1, t2, BIN=Bin )
;
; MODIFICATION HISTORY:
;        Written by:    Han Wen, September 1996.
;        30-SEP-1996    Misc. bug fixes. Accept array for Wt parameter.
;-
function POI_, Ncount, Mu

         lnP  = 0d0*Mu
         h    = WHERE(Mu gt 0,ngt0)
         if (ngt0 gt 0) then $
              lnP(h)  = Ncount(h)*ALOG(DOUBLE(Mu(h)))-Mu(h)-LNGAMMA(Ncount(h)+1)
         return, EXP(lnP)
end

function PROB2EXP, 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)

         b    = Bin
         t1   = Offset1
         t2   = Offset2
         p1   = 1./Mean1
         p2   = 1./Mean2
         p    = p1 + p2
         mu   = p*b

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

;   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 Prob = Prob # REPLICATE(1.,nq)

         K1   = LONG(t1/b)
         if (nWT gt 1) then begin
              h1   = WHERE( Wt eq K1*b )
              h1p1 = WHERE( Wt eq (K1+1)*b )
         endif else begin
              h1   = K1
              h1p1 = K1 + 1
         endelse

;   Probabilities of waiting times around the source 1 offset

         if (nmax lt K1) then goto, DONE
         if (h1(0) ne -1) then begin
              mu_k      = p*(K1*b     -(t1 + 0*t2))
              p0q_k     = (p1/(p*mu))*(poi_(0,mu_k+mu) + mu + mu_k - 1)
              Prob(h1,0)= p0q_k
         endif

         if (nmax lt (K1+1)) then goto, DONE
         if (h1p1(0) ne -1) then begin
              mu_kp1    = p*((K1+1)*b - (t1 + 0*t2))
              p0q_kp1   = (p1/(p*mu)) $
                          *(poi_(0,mu_kp1 + mu) -2*poi_(0,mu_kp1) $
                                                 + mu - mu_kp1 + 1)
              Prob(h1p1,0)= p0q_kp1
         endif

         if (nmax lt (K1+2)) then goto, DONE
         nbin = nmax - (K1+1)
         n0   = K1+2
         n    = (LINDGEN(nbin) + n0)

;   Probabilities of waiting times from source 1 and q events from source 2

         for  q=0L,nq-1 do begin
              K2   = LONG((t1 + q*t2)/b)
              nmin = K2 + 2
              if (nWT gt 1) then begin
                   h    = WHERE( (Wt ge nmin*b) and (Wt le nmax*b), nbin )
                   nm   = LONG(Wt(h)/b) # REPLICATE(1,q+1)
              endif else begin
                   nbin = nmax - nmin + 1
                   nm   = n(nmin-n0:nmax-n0) # REPLICATE(1,q+1)
              endelse

              z    = REPLICATE(1,nbin) # LINDGEN(q+1)
              mu_n = p*(nm*b -(t1 + q*t2))
              pNq  = (((p1/p)*(p2/p)^q)/mu) $
                       *((z + 1)*(  poi_(q-z,mu_n-mu)  $
                                   +poi_(q-z,mu_n+mu)  $
                                 -2*poi_(q-z,mu_n   )))
              if (q gt 0) then pNq = TOTAL(pNq,2)

              i    = q*(1 - KEYWORD_SET(Sum))
              if (nWT gt 1) then Prob(h,i) = Prob(h,i) + pNq $
              else Prob(nmin:nmax) = Prob(nmin:nmax)   + pNq

              if (MAX(pNq) lt eps) then goto, DONE
         endfor
DONE:    return, Prob
end