;+ ; 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