;+
; NAME: poisson_errorbar.pro
;
; PURPOSE: 
; For a given number of counts, determine the statistical uncertainty delta,
; based on Poisson statistics, such that the interval of:
; N - delta -- N + delta has a confidence level of X percent (defailt is 68.3%)
;
; CATEGORY: 
;   statistics, error propagation
;
; CALLING SEQUENCE: 
;   error = poisson_errorbar(counts, [confidence_level = confidence_level])
;
; INPUTS:
;   n_events -- number of detetcted counts
;
; KEYWORD PARAMETERS:
;   confidence_level -- the desired level of confince for the interval
;   double -- set this keyword to force double precision calculations
;
; OUTPUTS:
;   returns the error bar associated with the specified level of confidence
;
; COMMON BLOCKS:
;   poisson_errorbar_common -- contains variables cl, the specified confidence
;                              level and n, the number of counts
;
; RESTRICTIONS:
;   Since the code relies on FX_ROOT to find the root of confidence limits
;   function, the n_events must be a scalar
;
; MODIFICATION HISTORY:
;   2008-11-24  V1. AJS  Initial version
;-

FUNCTION confidence_limits_func, delta_n

COMMON poisson_errorbar_common, cl, n

upper_limit = n + delta_n
lower_limit = n - delta_n

; These are based on Eq. 4 & 5 of Geherls 1986 ApJ 303: 336-346
prob_above = 1.-chisqr_pdf(2.*upper_limit, 2. * n + 2.)
prob_below = chisqr_pdf(2.*lower_limit, 2.*n)

return, prob_above + prob_below + cl - 1.

END


FUNCTION poisson_errorbar, n_events, confindence_limit = confidence_limit, double = double

forward_function conf_limits_func

COMMON poisson_errorbar_common, cl, n

IF n_params() NE 1 THEN BEGIN
   print, '% POISSON_ERRORBAR: errorbar = poisson_errorbar(n_events, confindence_limit = confidence_limit)'
   return, -1
ENDIF

IF n_elements(n_events) NE 1 THEN BEGIN
   print, '% POISSON_ERRORBAR: n_events must be a scalar.'
   print, 'Returning.'
   return, -1
ENDIF

; Default level of confidnce is one sigma (68.3%)
IF NOT keyword_set(confidence_limit) THEN confidence_limit = gauss_pdf(1.) - gauss_pdf(-1.)

; Use double-precision if requested
IF keyword_set(double) THEN n = double(n_events) ELSE n = n_events

cl = confidence_limit

; Use the IDL routine FX_ROOT to solve for the Poissonian error bar
errorbar = fx_root([(sqrt(n)-1.) >.01, sqrt(n) > .1, sqrt(n)+1.] , 'confidence_limits_func')

return, errorbar
END