function ppd_src,ns,nb,sgrid,smap=smap,scred=scred,clev=clev,ppdstr=ppdstr,$
	funit=funit,asrc=asrc,abkg=abkg,agamma=agamma,bgamma=bgamma,$
	priorbg=priorbg,nsgrid=nsgrid,srcmin=srcmin,srcmax=srcmax,$
	verbose=verbose, _extra=e
;+
;function	ppd_src
;	returns the posterior probability density for the source
;	intensity, marginalized over background intensity, assuming
;	gamma-function priors (see van Dyk et al., 2001, ApJ 548, 224)
;
;	D:ns,nb ; I:asrc,abkg,f ; M:s,b ; p(M|DI)=p(M|I)*p(D|MI)/p(D|I)
;
;	r = abkg/asrc & A = nb+alfa_B & B = f*(beta_B+r)
;
;	p(b|I) = ((f*beta_B)^alfa_B/Gamma(alfa_B)) * b^(alfa_B-1) * exp(-b*beta_B*f)
;	p(nb|b,I) = ((r*f*b)^nb/Gamma(nb+1)) * exp(-b*f*r)
;	p(b|nb,I) = (B^A/Gamma(A)) * b^(A-1) * exp(-b*B)
;
;	I_k = b^(ns+A-k-1) * s^(k+alfa_S-1) * exp(-s*f*(1+beta_S) -b*(f+B)) /
;		(Gamma(k+1)*Gamma(ns-k+1))
;	I_kb = (Gamma(ns+A-k)/Gamma(k+1)/Gamma(ns-k+1)) * (f+B)^(-(ns+A-k)) *
;		s^(k+alfa_S-1) * exp(-s*f*(1+beta_S))
;	I_kbs = (Gamma(ns+A-k)*Gamma(k+alfa_S)/Gamma(k+1)/Gamma(ns-k+1)) *
;		(f+B)^(-(ns+A-k)) * (f*(1+beta_S))^(k+alfa_S)
;
;	p(s|I) = ((f*beta_S)^alfa_S/Gamma(alfa_S)) * s^(alfa_S-1) * exp(-s*beta_S*f)
;	p(ns|M,I) = (f*(s+b))^ns/Gamma(ns+1) * exp(-(s+b)*f)
;	p(sb|D,I) = p(b|nb,I)*p(s|I)*p(ns|M,I) / p(ns|I)
;		  = (Sum_{k=0}^{ns} I_k) / (Sum_{k=0}^{ns} I_kbs)
;
;	p(ns|I) = f^ns * (B^A/Gamma(A)) * ((f*beta_S)^alfa_S/Gamma(alfa_S)) *
;		Sum_{k=0}^{ns} I_kbs
;	p(s|D,I) = (Sum_{k=0}^{ns} I_kb) / (Sum_{k=0}^{ns} I_kbs)
;
;syntax
;	ppd=ppd_src(ns,nb,sgrid,smap=smap,scred=scred,clev=clev,ppdstr=ppdstr,$
;	funit=funit,asrc=asrc,abkg=abkg,agamma=agamma,bgamma=bgamma,$
;	priorbg=priorbg,nsgrid=nsgrid,srcmin=srcmin,srcmax=srcmax,verbose=verbose)
;
;parameters
;	ns	[INPUT; required] observed counts in source region
;	nb	[INPUT; required] observed counts in background region
;		* NS and NB must be scalar -- if vector, adds up all the elements
;	sgrid	[OUTPUT] the output grid over which the PPD is calculated
;		* a regular grid that uses NSGRID, SRCMIN, and SRCMAX
;
;keywords
;	smap	[OUTPUT] the mode of the distribution, corresponding to the
;		maximum a posteriori value of the PPD
;	scred	[OUTPUT] a two-tailed credible range on the source intensity,
;		defined to include CLEV/2 of the area under the PPD above
;		and below SMAP
;		* NOTE: what this means is, if the distribution is skewed,
;		  then the total area under PPD as defined by SCRED is not
;		  necessarily equal to CLEV
;		* The one-tailed credible regions, that is, the range that
;		  encloses an area CLEV under the PPD, are returned in the
;		  structure PPDSTR.  these regions are determined such that
;		  the mode is always _included_ within the region, but is not
;		  necessarily at the center
;	clev	[INPUT] a nominal level at which to determine bounds and
;		get credible regions
;		* default is 0.68
;		* if < 0, abs(CLEV) is used
;		* if > 1 and < 100, then assumed to be given as a percentage
;		* if > 100, then 1-1/CLEV is used
;	ppdstr	[OUTPUT] an anonymous structure containing the fields
;		{
;		  GRID : same as SGRID
;		  MODE : mode, same as SMAP
;		  MEAN : mean of the distribution
;		  SIG : s.d. on the MEAN
;		  MEDIAN : the median of the distribution
;		  CLEV : the probability level defined by the keyword CLEV
;		  HPD1 : the highest posterior density "1-sigma" credible range,
;			 the interval that encloses 68% of the area and also
;			 the highest values of the probability density
;		  HPD2 : as HPD1, for "2-sigma", 95% of the area
;		  HPD3 : as HPD1, for "3-sigma", 99.7% of the area
;		  HPD0 : as HPD1, but corresponding to a fraction CLEV
;		  CRED1 : the "1-sigma" credible range, corresponding to 68%
;		  	  of the area under the PPD that _includes_ the MODE
;		  CRED2 : as CRED1, for "2-sigma", 95% of the area
;		  CRED3 : as CRED1, for "3-sigma", 99.7% of the area
;		  CRED0 : as CRED1, but corresponding to a fraction CLEV of the
;			  area under the PPD that includes the MODE
;		  RNG1 : the two-tailed "1-sigma" range, corresponding to a
;		  	 fractional area under the PPD of 0.34 above and below
;		  	 the MODE
;		  RNG2 : as for RNG1, but the "2-sigma" range, 0.475 of the area
;		  RNG3 : as for RNG1, but the "3-sigma" range, 0.4985 of the area
;		  RNG0 : as for RNG1, but CLEV/2 of the area, same as SCRED
;		  lnpD : the log of the predictive probability, ALOG(p(ns|I))
;		}
;	funit	[INPUT] a factor to _divide_ the counts to convert to different units
;		e.g., ph/cm^2/s
;		* assumed to be the same for both source and background counts
;		  (any differences in exposure times can be encoded in ASRC and ABKG)
;		* default is 1.0
;	asrc	[INPUT] area (or exposure time) for the background contaminated
;		source observation, NS
;		* default is 1.0
;	abkg	[INPUT] area (or exposure time) for the background observation, NB
;		* default is 1.0
;	agamma	[INPUT] alpha parameter for the source intensity prior, assumed to
;		be a gamma distribution (see PROB_GAMMADIST)
;		* default is 1.0 (also the hardcoded minimum)
;		* generally taken to be the number of *expected* source counts+1
;	bgamma	[INPUT] beta parameter for the source intensity prior
;		* default is 0.0 (also the hardcoded minimum)
;		* generally the ratio of that area in which the AGAMMA-1 counts are
;		  obtained in, to ASRC
;		* note that the default values of AGAMMA and BGAMMA make for a
;		  very non-informative, but improper, prior
;	priorbg	[INPUT] the expected number of counts due to the background,
;		a number that is known from prior information and NOT derived
;		from NB itself.
;		* a Gamma function prior is assumed on the background strength, with
;		  alpha=priorbg+1; beta=abkg/asrc, but if priorbg=0, then beta=0
;	nsgrid	[INPUT] number of points in output grid
;		* default=101L
;	srcmin	[INPUT] minimum in source intensity value to consider
;	srcmax	[INPUT] maximum in source intensity value to consider
;		* beware that setting [SRCMIN,SRCMAX] to too small a range will cause
;		  it to act as a prior over and above the gamma prior.  this may not be
;		  desirable.
;		* the default range is the widest of:
;		  -- NS +- 10*sqrt(NS)
;		  -- (NS-(ASRC/ABKG)*NB) +- 10*sqrt(NS+(ASRC/ABKG)^2*NB) (if ABKG > 0)
;		  -- (agamma/bgamma) +- 10*sqrt(agamma/bgamma^2) (if BGAMMA > 0)
;		  as long as SRCMIN isn't -ve or SRCMAX is smaller than 20
;	verbose	[INPUT] controls chatter
;	_extra	[JUNK] here only to prevent crashing the program
;
;subroutines
;	PROB_GAMMADIST
;
;history
;	vinay kashyap (May02)
;	added keyword FUNIT (VK; Oct03)
;	added HPD1,HPD2,HPD3,HPD0 to the output (VK; Mar06)
;-

;	usage
ok='ok' & np=n_params() & nns=n_elements(ns) & nnb=n_elements(nb)
if np lt 2 then ok='Insufficient parameters' else $
 if nns eq 0 then ok='NS: undefined' else $
  if nnb eq 0 then ok='NB: undefined'
if ok ne 'ok' then begin
  print,'Usage: ppd=ppd_src(ns,nb,sgrid,smap=smap,scred=scred,clev=clev,ppdstr=ppdstr,$'
  print,'       funit=funit,asrc=asrc,abkg=abkg,agamma=agamma,bgamma=bgamma,$'
  print,'       priorbg=priorbg,nsgrid=nsgrid,srcmin=srcmin,srcmax=srcmax,verbose=verbose)'
  print,'  return posterior probability density of source intensity, marginalized'
  print,'  over background intensities'
  if np ne 0 then message,ok,/info
  return,-1L
endif

;	check input parameters
if nns gt 1 then Ds=total(ns) else Ds=ns[0]
if nnb gt 1 then Db=total(nb) else Db=nb[0]
if Ds lt 0 or Db lt 0 then begin
  message,'cannot handle -ve observed counts; quitting',/info & return,-1L
endif

;	check input keywords
vv=0L & if keyword_set(verbose) then vv=long(verbose[0]) > 1
;
crlev=0.68 & if keyword_set(clev) then crlev=0.0+clev[0]
if crlev lt 0 then crlev=abs(crlev)
if crlev ge 1 and crlev lt 100 then crlev=crlev/100.
if crlev ge 100 then crlev = 1.0D - 1.0D/crlev
;
facu=1.0 & if keyword_set(funit) then facu=1.0*funit[0]
;
srcar=1.0 & if keyword_set(asrc) then srcar=0.0+asrc[0]
if srcar le 0 then begin
  message,'Cannot handle -ve ASRC; quitting',/info & return,-1L
endif
bkgar=1.0 & if keyword_set(abkg) then bkgar=0.0+abkg[0]
if bkgar lt 0 then begin
  message,'Cannot handle -ve ABKG; quitting',/info & return,-1L
endif
;
alfa_S=1. & if keyword_set(agamma) then alfa_S=0.0+agamma[0]
if alfa_S lt 1 then begin
  if vv gt 0 then message,'AGAMMA cannot be < 1; resetting',/info
  alfa_S = 1.
endif
;
beta_S=0. & if n_elements(bgamma) gt 0 then beta_S=0.0+bgamma[0]
if beta_S lt 0 then begin
  if vv gt 0 then message,'BGAMMA cannot be < 0; resetting',/info
  beta_S = 0.
endif
;
alfa_B=1. & beta_B=0.
if keyword_set(priorbg) then begin
  alfa_B=priorbg[0]+1.>1.
  beta_B=bkgar/srcar
endif
;
msgrid=101L & if keyword_set(nsgrid) then msgrid=long(abs(nsgrid[0])) > 2
;
s_ML=Ds/facu & s_MLe=(sqrt(s_ML+0.75)+1.)/facu
smin1=0 & smax1=20/facu & smin=smin1 & smax=smax1
smin2=s_ML-10*s_MLe > smin1
smax2=s_ML+10*s_MLe > smax1
smin=smin2 & smax=smax2
if bkgar gt 0 then begin
  s_ML=Ds/facu-(srcar/bkgar)*Db/facu
  s_MLe=(sqrt((sqrt(Ds+0.75)+1.)^2+(srcar/bkgar)^2*(sqrt(Db+0.75)+1.)^2))/facu
  smin3=s_ML-10*s_MLe > smin1
  smax3=s_ML+10*s_MLe > smax1
  if smin3 lt smin2 then smin=smin3
  if smax3 gt smax2 then smax=smax3
endif
if beta_S gt 0 then begin
  smin4=(alfa_S/beta_S) - 10*sqrt(alfa_S/beta_S^2) > smin1
  smax4=(alfa_S/beta_S) + 10*sqrt(alfa_S/beta_S^2) > smax1
  if smin4 lt smin2 then smin=smin4
  if smax4 gt smax2 then smax=smax4
endif
if keyword_set(srcmin) then begin
  if srcmin[0] gt smin then begin
    if vv gt 0 then message,'SRCMIN > default; presumably you know what you are doing',/info
  endif
  smin=srcmin[0]
endif
if keyword_set(srcmax) then begin
  if srcmax[0] lt smax then begin
    if vv gt 0 then message,'SRCMAX < default; presumably you know what you are doing',/info
  endif
  smax=srcmax[0]
endif
if smin gt smax then begin
  message,"SRCMIN > SRCMAX?  can't be havin' wi' that",/info
  jnk=smax & smax=smin & smin=jnk
endif
deltS=(smax-smin)/(msgrid-1L) & sgrid=findgen(msgrid)*deltS+smin

;	output
lnppd=fltarr(msgrid) & ppd=dblarr(msgrid)

;	now compute the coefficients
kk=lindgen(Ds+1L)
lnI_kbs = lngamma(Ds+Db+(alfa_B-1)-kk+1L) +$
	lngamma(kk+alfa_S) -$
	lngamma(kk+1L) -$
	lngamma(Ds-kk+1L) -$
	(Ds+Db+(alfa_B-1)-kk+1L)*alog(1.0+(bkgar/srcar)+beta_B) -$
	(Ds+Db+(alfa_B-1)-kk+1L)*alog(facu) -$
	(kk+alfa_S)*alog(1.+beta_S) -$
	(kk+alfa_S)*alog(facu)
normkbs=max(lnI_kbs) & dkbs=lnI_kbs-normkbs
sumkbs=total(exp(dkbs)) & sumkbs=alog(sumkbs)+normkbs
lnI_kb = lngamma(Ds+Db+(alfa_B-1)-kk+1L) -$
	lngamma(kk+1L) -$
	lngamma(Ds-kk+1L) -$
	(Ds+Db+(alfa_B-1)-kk+1L)*alog(1.0+(bkgar/srcar)+beta_B) -$
	(Ds+Db+(alfa_B-1)-kk+1L)*alog(facu)

;	r = abkg/asrc & A = nb+alfa_B & B = f*(beta_B+r)
;	I_kbs = (Gamma(ns+A-k)*Gamma(k+alfa_S)/Gamma(k+1)/Gamma(ns-k+1)) *
;		(f+B)^(-(ns+A-k)) * (f*(1+beta_S))^(k+alfa_S)
;	p(D|I) = f^ns * (B^A/Gamma(A)) * ((f*beta_S)^alfa_S/Gamma(alfa_S)) *
;		Sum_{k=0}^{ns} I_kbs
if beta_S gt 0 then p_D_I=alfa_S*alog(beta_S) else p_D_I=0.
p_D_I = p_D_I + ns*alog(facu) +$
	(Db+alfa_B)*alog(beta_B+(bkgar/srcar)) +$
	(Db+alfa_B)*alog(facu) -$
	lngamma(Db+alfa_B) +$
	alfa_S*alog(facu) -$
	lngamma(alfa_S) +$
	sumkbs
;p_D_I = exp(p_D_I)

;	for each source intensity grid point, compute the ppd
for i=0L,msgrid-1L do begin
  s=sgrid[i]
  if s eq 0 then begin
    tmp=lnI_kb[0] & normkb=tmp & sumkb=1.
    lnppd[i]=lnI_kb[0]-sumkbs ;lnI_kbs[0]
  endif else begin
    tmp=lnI_kb+(kk+alfa_S-1L)*alog(s)
    normkb=max(tmp) & dkb=tmp-normkb & sumkb=total(exp(dkb))
    lnppd[i]=alog(sumkb)+normkb-sumkbs-(1.+beta_S)*s*facu
  endelse
  if vv ge 150 then begin
    print,'ln{p('+strtrim(s,2)+'|'+strtrim(Ds,2)+','+strtrim(Db,2)+',I)}='+strtrim(lnppd[i],2)
    if vv gt 200 then stop,s
  endif
endfor
lnormppd=max(lnppd) & ppd=exp(lnppd-lnormppd)
normppd=total(ppd*deltS) & ratnorm=exp(-lnormppd)/normppd	;just in case
ppd=ppd*exp(lnormppd)
;
if vv gt 90 then print,'normalization (should be v/s is) =',exp(-lnormppd),normppd
if ratnorm gt 1.5 or ratnorm lt 1/1.5 and vv gt 0 then begin
  message,'	grid appears to be non optimal;',/info
  message,'	increase NSGRID, or fiddle with SRCMIN and SRCMAX',/info
endif

;	find the mode
if facu eq 1 then spline_p,sgrid,ppd,ss,pp,interval=deltS/100. else begin
  if facu lt 100. then spline_p,sgrid,ppd,ss,pp else begin
    ss=sgrid & pp=ppd
  endelse
endelse
pp=pp>0
jnk=max(pp,imx) & smap=ss[imx]

;	find the highest posterior density intervals
os=reverse(sort(pp)) & cpp=total(pp[os],/cumul) & cpp=cpp/max(cpp) & xx=ss[os]
om=where(xx le smap,mom) & hpdm=fltarr(4)+smap
op=where(xx ge smap,mop) & hpdp=fltarr(4)+smap
if mom gt 1 then hpdm=interpol(xx[om],cpp[om],[crlev,0.68,0.95,0.997])
if mop gt 1 then hpdp=interpol(xx[op],cpp[op],[crlev,0.68,0.95,0.997])

;	find the mean
;smean=total(ss*pp*deltS/100.) & ssig=sqrt(total(ss^2*pp*deltS/100.)-smean^2)
smean=total(sgrid*ppd*deltS) & ssig=sqrt(total(sgrid^2*ppd*deltS)-smean^2)
if vv gt 90 and finite(ssig) le 0 then stop,smean,ssig

;	make cumulative distribution
sp=findgen(msgrid+1L)*deltS+smin-deltS/2.>0 & cp=fltarr(msgrid+1L)
for i=1L,msgrid do cp[i]=cp[i-1L]+ppd[i-1L] & cp=cp/max(cp)
smedian=interpol(sp,cp,0.5)

;	make cumulative 2-tailed distributions
cmode=interpol(cp,sp,smap)
ucp=cp-cmode & ou=where(ucp ge 0,mou)
if mou gt 1 then begin
  usp=[smap,sp[ou]] & ucp=[0.,ucp[ou]] & ucp=ucp/max(ucp)
endif else begin
  usp=[smap,sp[msgrid]] & ucp=[0.,1.]
endelse
lcp=-(cp-cmode) & ol=where(lcp ge 0,mol)
if mol gt 1 then begin
  lsp=[sp[ol],smap] & lcp=[lcp[ol],0.] & lcp=lcp/max(lcp)
endif else begin
  lsp=[sp[0],smap] & lcp=[1.,0.]
endelse
levs=[crlev,0.68,0.95,0.997]
uu=interpol(usp,ucp,levs) & ll=interpol(lsp,lcp,levs)
scred=[ll[0],uu[0]]

;	find credible regions that have the right area AND
;	include the mode in all cases
ulev=0.5D + levs/2.D - (0.5D - cmode)
llev=0.5D - levs/2.D - (0.5D - cmode)
oo=where(llev lt 0,moo)
if moo gt 0 then begin
  ulev[oo]=ulev[oo]-llev[oo] & llev[oo]=0.
endif
oo=where(ulev gt 1,moo)
if moo gt 0 then begin
  llev[oo]=llev[oo]-(ulev[oo]-1.) & ulev[oo]=1.
endif
ureg=interpol(sp,cp,ulev) < smax
lreg=interpol(sp,cp,llev) > smin

;	make the summary structure
ppdstr=create_struct('GRID',sgrid,'MODE',smap,'MEAN',smean,'SIG',ssig,$
	'MEDIAN',smedian,'CLEV',crlev,$
	'HPD1',[hpdm[1],hpdp[1]],'HPD2',[hpdm[2],hpdp[2]],$
	'HPD3',[hpdm[3],hpdp[3]],'HPD0',[hpdm[0],hpdp[0]],$
	'CRED1',[lreg[1],ureg[1]],'CRED2',[lreg[2],ureg[2]],$
	'CRED3',[lreg[3],ureg[3]],'CRED0',[lreg[0],ureg[0]],$
	'RNG1',[ll[1],uu[1]],'RNG2',[ll[2],uu[2]],$
	'RNG3',[ll[3],uu[3]],'RNG0',[ll[0],uu[0]],$
	'lnpD',p_D_I)

;	report
tt0='S = '+strtrim(smap,2)+$
	' +'+strtrim(scred[1]-smap,2)+',-'+strtrim(smap-scred[0],2)+$
	' @'+strtrim(crlev,2)
tt1='MaxLik='+strtrim(s_ML,2)+'+-'+strtrim(s_MLe,2)+$
	' ; MEAN='+strtrim(smean,2)+'+-'+strtrim(ssig,2)+$
	' ; MEDIAN='+strtrim(smedian,2)
if vv gt 0 then begin
  message,tt0,/info
  if vv gt 1 then message,tt1,/info,/noname
  if vv gt 3 then message,'AGAMMA='+strtrim(alfa_S,2)+' ; '+$
	'BGAMMA='+strtrim(beta_S,2)+' ; '+$
	'ASRC='+strtrim(srcar,2)+' ; '+$
	'ABKG='+strtrim(bkgar,2)+' ; '+$
	'FUNIT='+strtrim(facu,2),/informational,/noname
  if vv gt 20 then begin
    plot,sgrid,ppd,charsize=2,$
	xtitle='!4k!3!uS!n',ytitle='!3p(!4k!3!uS!n|n!dS!n,n!dB!n,I)',$
	title=strtrim(smap,2),subtitle=tt0
    tmp=prob_gammadist(sgrid,agamma=alfa_S,bgamma=beta_S)
    oplot,sgrid,tmp*total(ppd)/total(tmp),line=1,thick=2
  endif
  if vv ge 100 then stop,'HALTing.  type .CON to continue'
endif else begin
  if np eq 2 and not arg_present(smap) and not arg_present(ppdstr) then begin
    message,tt0,/info
    message,tt1,/info,/noname
  endif
endelse

return,ppd
end