function livetime, countrate, tau
on_error,2
;+
;
; NAME:  Livetime
;
; PURPOSE: Return livetime fraction
;
; CATEGORY: Data reduction, analysis, time histories, spectroscopy
;       BATSE
; CALLING SEQUENCE: fraction = livetime( countrate [,tau] )
;
; CALLS TO: Interpol
;
; INPUTS:
;	Countrate - Counts/sec: Total counts / realtime second.
;
; OPTIONAL INPUTS:
;	Tau       - Scalar or vector recovery time per event
;		  - default is value for DISCLA, 1.1e-6 seconds
;                 - must be scalar or same number of elements as countrate
;
; OUTPUTS:  livetime fraction
;
; PROCEDURE: Inverts the equation n*tau = r*tau exp(-(r*tau)) where
;		r is the true rate, and n is the observed rate
;		The livetime fraction is defined as n/r.  This algorithm
;		only inverts one side of the equation which has two
;		values of r for each value of n.
;
; MODIFICATION HISTORY: ras, 3 Sep 93
;	fixed at low count rate, 26 Sep 93, ras
;
;-

checkvar, tau, 1.1e-6
taus = tau + countrate*0.0

if n_elements(taus) ne n_elements(countrate) then begin
	print,'Error! Tau must be scalar or same length as Countrate!'
	stop
endif

result = countrate * 0.0
crt = countrate * taus < 1./exp(1.)

w1 = where(crt le .1, nw1) 
;This is the simplest approximation, live time change is linear with counts,
;with a little non-linear correction
if nw1 ge 1 then result(w1) = crt(w1)*(1. + crt(w1) +1.75 * crt(w1)^2)
	
wlt = where(crt le 0.18 and crt gt 0.1, nlt)
;
; For a measured rate less than half the maximum (exp(1)), use a 2nd order fit
; "   " "        "    greater than half, interpolate r*exp(-r*tau)
;
if nlt ge 1 then begin
	coef1 = reform( [1.34932544d-03, 0.93221443,1.7223753],1,3)
	result(wlt) = poly( crt(wlt), coef1)
endif

wgt = where(crt gt 0.18, ngt)
if ngt ge 1 then begin
	r=(findgen(25)+1.0)*.04
	x=   r /exp(r)
	y=r
	result(wgt) = interpol( y,x,crt(wgt) )
endif

result = f_div(crt,result) ;livetime fraction
return, result
end