;+
; NAME:
; PROBHEXP
;
; PURPOSE:
; This function determines the probability function of the binned
; waiting times between successive events for a hyperexponential
; distribution distorted by a non-extended dead time,
;
; f(t) = P1*r1*exp(-r1*(t-tau)) + (1-P1)*r2*exp(-r2*(t-tau))
;
; CATEGORY:
; Math.
;
; CALLING SEQUENCE:
;
; Result = PROBHEXP( Wtmax, Mean1, Mean2, Offset, Prob1 )
;
; 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 the FIRST exponential.
;
; Mean2: The average waiting time for the SECOND exponential.
;
; Offset: The non-extended dead time.
;
; Prob1: Probability of generating events from the FIRST exponential.
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
;
; BIN: Time resolution of the waiting times, (1=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].
; NOTE: To calculate the probability, one must take into account that
; events can fall anywhere within the time resolution, BIN.
;
; EXAMPLE:
; For times binned with a bin time interval of 1, then the possible
; binned waiting times are [0,1,2..] corresponding to REAL waiting
; times of ([0,1], [0,2], [1,3], ...]. If we want to determine the
; probabilities of binned waiting times of [0,10] when the average
; number of counts = 1/bin (without dead time) and the deadtime
; is 2.5, with 15% of the photons being affected by this
; software/hardware bug with a scale factor or 2, then
;
; avg_cts = 1.0
; avg_cts2 = 0.5
; tavg = 1./avg_cts
; tavg2 = 1./avg_cts2
; dead_time = 2.5
; max_waiting_time = 10
; fraction_unaffected = 0.15
; Probs = PROBHEXP( max_waiting_time, tavg, tavg2,
; dead_time, fraction_unaffected )
;
; MODIFICATION HISTORY:
; Written by: Han Wen, October 1996.
; 06-DEC-1996 Renamed from PROBWBUG -> PROBHEXP, Changed Frac ->
; (1 - Prob1).
;-
function PROBHEXP, Wt, Mean1, Mean2, Offset, Prob1, BIN=Bin
NP = N_PARAMS()
if (NP ne 5) then message, $
'Must be called with 5 parameters: Wt, Mean1, Mean2, Offset, Prob1'
p1 = PROBEXP( Wt, Mean1, Offset, BIN=Bin )
p2 = PROBEXP( Wt, Mean2, Offset, BIN=Bin )
return, Prob1*p1 + (1-Prob1)*p2
end