pro binomerr,ntotal,ndetect

; pro to calculate binomial upper and lower error bars
;	based on the formalism in burgasser et al. 2003
;	apj, 586, 512, appendix A.
;
; inputs:
;	ntotal = number of total observed sources
;	ndetect = number of positive detections
;
; outputs:
;	prints to screen upper and lower error bars that
;	together encompass 68% probability range (ie
;	equivalent to 1sigma gaussian errors).
;
; 15 jun 2006 // det
;

; can adjust this number if you like. 
; 100 steps gives precision of 0.01.
nsteps = 300

intermediatesum = make_array(ndetect+1)
finalsum = make_array(nsteps)

udiffkeep = 1.0
ldiffkeep = 1.0

; set up x array
x = findgen(nsteps)/float(nsteps)

; at each x ...
for j = 0,n_elements(x)-1 do begin

; ... evaluate equation A3 from burgasser at each i
   for i = 0,ndetect do begin

      factor1 = ( ( factorial(ntotal+1) ) / ( (factorial(i)) * (factorial(ntotal+1-i)) ) )

      factor2 = float(x[j])^float(i)

      factor3 = (1.0 - x[j])^(float(ntotal+1-i))

      intermediatesum[i] = factor1 * factor2 * factor3

   endfor

; sum over i to find total sum for each x

   finalsum[j] = total(intermediatesum)

; if this finalsum is close to the desired 1sigma values,
; save the indices

   ldiff = abs(finalsum[j] - 0.84)
   udiff = abs(finalsum[j] - 0.16)

   if (udiff lt udiffkeep) then begin
      udiffkeep = udiff
      jup = j
   endif
   if (ldiff lt ldiffkeep) then begin
      ldiffkeep = ldiff
      lup = j
   endif

endfor

; now print the best solutions

print,'upper limit = ',x[jup],' where prob = ',finalsum[jup]
print,' best value = ',float(ndetect)/float(ntotal)
print,'lower limit = ',x[lup],' where prob = ',finalsum[lup]

end