function prob_gammadist,x,agamma=agamma,bgamma=bgamma,verbose=verbose, _extra=e
;+
;function	prob_gammadist
;	compute the gamma distribution, which is a continuous distribution
;	on the real line with density function
;		f(x) = b^a x^(a-1) exp(-bx) / Gamma(a)
;	(see eqn 4 of van Dyk et al., 2001, ApJ 548, 224)
;	this distribution has mean (a/b) and variance (a/b^2) for a,b>0
;
;	This is equivalent to a Poisson distribution, and is often
;	used as a prior function, with a and b representing our prior
;	expectation that a-1 counts should be observed in an area
;	b*src_area; a non-informative prior therefore has a=1 and b=0,
;	corresponding to 0 counts in 0 area.
;
;syntax
;	g=prob_gammadist(x,agamma=agamma,bgamma=bgamma,verbose=verbose)
;
;parameters
;	x	[INPUT; required] array of abscissa values at which to
;		compute the Gamma distribution
;		* for x < 0, the function will return NaNs
;
;keywords
;	agamma	[INPUT; default=1] "alpha" parameter for the function
;		* AGAMMA must be > 0, or else the function will return Infs
;	bgamma	[INPUT; default=0] "beta" parameter for the function
;		* BGAMMA should be .GE. 0, or else who knows what will happen
;	verbose	[INPUT; default=0] controls chatter
;	_extra	[JUNK] here only to prevent crashing the program
;
;history
;	vinay kashyap (Aug01)
;-

;	usage
ok='ok' & np=n_params() & nx=n_elements(x)
if np eq 0 then ok='Insufficiant parameters' else $
 if nx eq 0 then ok='X is not defined'
if ok ne 'ok' then begin
  print,'Usage: g=prob_gammadist(x,agamma=agamma,bgamma=bgamma,verbose=verbose)'
  print,'  compute the gamma distribution,'
  print,'  g(x) = bgamma^agamma x^(agamma-1) exp(-bgamma x) / Gamma(agamma)'
  if np ne 0 then message,ok,/info
  return,-1L
endif

;	keywords
aa=1. & bb=0.
if n_elements(agamma) ne 0 then aa=agamma[0]+0. else message,'Setting A=1',/info
if n_elements(bgamma) ne 0 then bb=bgamma[0]+0. else message,'Setting B=0',/info
vv=0 & if keyword_set(verbose) then vv=long(verbose[0]) > 1

;	reformat input
xx=[x[*]]
oo=where(xx lt 0,moo) & ok=where(xx gt 0,mok) & o0=where(xx eq 0,mo0)

;	initialize output
if moo ne 0 and vv gt 0 then message,strtrim(moo,2)+$
	' points will be set to NaN',/info
g=fltarr(nx) & if moo ne 0 then g[oo]=!values.F_NAN

;	compute function
if bb ne 0 then tmp1=aa*alog(bb) else tmp1=0.
tmp2=fltarr(nx) & if mok gt 0 then tmp2[ok]=(aa-1.)*alog(xx[ok])
tmp3=fltarr(nx) & if mok gt 0 then tmp3[ok]=-bb*xx[ok]
if mok gt 0 then g[ok]=tmp1+tmp2[ok]+tmp3[ok]-lngamma(aa)
if mo0 gt 0 then g[o0]=tmp1+tmp2[o0]+tmp3[o0]-lngamma(aa)
if mok gt 0 then g[ok]=exp(g[ok])
if mo0 gt 0 then g[o0]=exp(g[o0])
if aa ne 1 and mo0 gt 0 then g[o0]=0.

if vv gt 20 and !d.NAME eq 'X' then plot,x,g,xtitle='!3X',$
	ytitle='!3G(X;!4a,b!3)',title='!4a!3='+strtrim(aa,2)+' ; !4b!3='+$
	strtrim(bb,2)
if vv gt 100 then stop,'Halting.  type .CON to continue'

return,g
end