;+ ; NAME: ; EXPDEV ; ; PURPOSE: ; Returns an exponentially distributed, positive, random deviate ; using RANDOMU( seed ) as the source of uniform deviates. The ; default unit mean and zero offset of the deviates may be modified ; by setting the appropriate keywords. ; ; CATEGORY: ; Math. ; ; CALLING SEQUENCE: ; ; Result = EXPDEV( Seed [,N] ) ; ; INPUTS: ; Seed: Seed for the random number generator. (It can be undefined) ; ; OPTIONAL INPUTS: ; N: Number of deviates to return, (1=Default). ; ; OPTIONAL INPUT KEYWORD PARAMETERS: ; MEAN: Mean of the deviates, (1=Default). ; ; OFFSET: Offset of the deviates, (0=Default). The resulting ; exponential distribution and its mean is offset by ; this amount. ; ; SUM: Open upper limit to the sum of all the deviates, ; (i.e. TOTAL(EXPDEV(Seed,..)) < SUM), (0=Default). ; Deviates will be generated until their sum is equal to ; or greater than this cutoff. If this keyword is specified, ; then the number of deviates generated is returned in the ; optional argument, N. ; ; CUTOFF: Minimum possible value of the deviates, (1e-37=Default). ; If set, this value MUST be > 0. ; ; MODIFICATION HISTORY: ; Written by: Han Wen, October 1995. ; 02-NOV-1995 Formerly, EXPDT, changed keywords TAU -> OFFSET, ; AVG -> MEAN, added CUTOFF keyword. ; 24-NOV-1995 Return a scalar when called with 1 argument, ; (i.e. EXPDEV(Seed)). ; 03-JAN-1996 Bugfixes: TOTAL(dts)->TOTAL(dts,/DOUBLE), ; N = round(N) -> round(N) > 1 ;- function EXPDEV, Seed, N, MEAN=Mean, OFFSET=Offset, CUTOFF=Cutoff, SUM=Sum NP = N_PARAMS() if (NP lt 2) then N=1 if (N_ELEMENTS(Mean) eq 0) then Mean = 1 $ else if (Mean le 0) then message,'MEAN must be > 0.' if (N_ELEMENTS(Offset) eq 0) then Offset = 0 if (N_ELEMENTS(Cutoff) eq 0) then Cutoff = 1.e-37 $ else if (Cutoff le 0) then message,'CUTOFF must be > 0.' if (N_ELEMENTS(Sum) eq 0) then Sum = 0 $ else if (Sum le Offset) then begin N = 0 dts = -1 return, dts endif if (Sum eq 0) then begin pdf = 1./Mean * RANDOMU(Seed,N) > Cutoff dts = Offset -Mean*alog(pdf*Mean) if (NP lt 2) then dts=dts(0) endif else begin rho = 1./Mean N = rho*Sum/(1.+rho*Offset) ; Estimate of how many N = round(N) > 1 ; deviates to generate nadd = 0L pdf = rho * RANDOMU(Seed,N) > Cutoff dts = Offset -Mean*alog(pdf*Mean) tsum = TOTAL(dts,/DOUBLE) ; If sum of current deviates is LESS than SUM then ; generate additional deviates one at a time until the ; sum of all the deviates is > SUM. if (tsum lt Sum) then begin repeat begin pdf = rho * RANDOMU(Seed) > Cutoff new_dt = Offset -Mean*alog(pdf*Mean) dts = [dts,new_dt] tsum = tsum + new_dt nadd = nadd + 1 endrep until (tsum ge Sum) dts = dts(0:N+nadd-2) if (N eq 0) and (nadd gt 2) then begin N = nadd-2 dts = dts(1:N) endif else N = N + nadd-1 ; If the sum of current deviates is GREATER than or EQUAL to ; SUM then eliminate deviates one at a time until the ; sum of the remaining deviates is < SUM. endif else begin k = N repeat begin k = k-1 tsum = tsum - dts(k) endrep until (tsum lt Sum) N = k if (N gt 0) then dts = dts(0:N-1) $ else dts = -1 endelse if (N le 1) then dts=dts(0) endelse return,dts end