function colbehr,Ssrc,Msrc,Hsrc, Sbkg=Sbkg,Mbkg=Mbkg,Hbkg=Hbkg,$
	Sarea=Sarea,Marea=Marea,Harea=Harea, Seff=Seff,Meff=Meff,Heff=Heff,$
	Sidx=Sidx,Midx=Midx,Hidx=Hidx, Sscl=Sscl,Mscl=Mscl,Hscl=Hscl,$
	Stable=Stable,Mtable=Mtable,Htable=Htable,$
	nsim=nsim,nburnin=nburnin,post=post,level=level,$
	outputF=outputF,BEHRdir=BEHRdir,verbose=verbose, _extra=e
;+
;function	colbehr
;	computes hardness ratios based on 3-band data by invoking
;	BEHR (Bayesian Estimate of Hardness Ratios) twice, and
;	returns the inputs as well as the following quantities
;	and their credible ranges in a structure:
;	  {S, M, H, lS=log10(S), lM=log10(M), lH=log10(H),
;	  SpM=S+M, MpH=M+H, SpH=S+H, T=S+M+H,
;	  SmM=S-M, MmH=M-H, SmH=S-H,
;	  R1=S/M, R2=M/H, R3=S/H,
;	  C1=C_SM, C2=C_MH, C3=C_SH,
;	  HR1=(S-M)/(S+M), HR2=(M-H)/(M+H), HR3=(S-H)/(S+H),
;	  HRA=(S-M)/(S+M+H), HRB=(M-H)/(S+M+H), HRC=(S-H)/(S+M+H)}
;
;	Warning: Because of the necessity of combining two separate
;	runs of BEHR, only the MCMC option is used.  Thus, if the
;	length of the chain is small, the computed values may be
;	subject to computational instability.
;
;	Warning: BEHR assumes that the passbands in question
;	do not overlap, and that the counts input to the program
;	are statistically independent.  It is up to the users to
;	ensure the validity of this assumption.  No checks are
;	made to verify it.
;
;	Reference:
;	"Bayesian Estimation of Hardness Ratios: Modeling and Computations",
;	  Park, T., Kashyap, V.L., Siemiginowska, A., van Dyk, D., Zezas, A.,
;	  Heinke, C., and Wargelin, B., 2006, ApJ, 652, 610
;	http://hea-www.harvard.edu/AstroStat/BEHR/
;
;syntax
;	behr3=colbehr(Ssrc,Msrc,Hsrc,Sbkg=Sbkg,Mbkg=Mbkg,Hbkg=Hbkg,$
;	Sarea=Sarea,Marea=Marea,Harea=Harea,$
;	Seff=Seff,Meff=Meff,Heff=Heff,$
;	Sidx=Sidx,Midx=Midx,Hidx=Hidx,$
;	Sscl=Sscl,Mscl=Mscl,Hscl=Hscl,$
;	Stable=Stable,Mtable=Mtable,Htable=Htable,$
;	/post,level=level,nsim=nsim,nburnin=nburnin,hpd=hpd,$
;	outputf=outputf,BEHRdir=BEHRdir, verbose=verbose)
;
;parameters
;	Ssrc	[INPUT; required] counts in source region in the soft (S) band
;	Msrc	[INPUT; required] counts in source region in the medium (M) band
;	Hsrc	[INPUT; required] counts in source region in the hard (H) band
;		* can be arrays; if so, the array with the most elements
;		  determines the size of the output and the shortfalls in the
;		  others, if any, are made up by replicating the first elements
;
;keywords
;	Sbkg	[INPUT] counts in background region in the S band
;	Mbkg	[INPUT] counts in background region in the M band
;	Hbkg	[INPUT] counts in background region in the H band
;		* if not given, assumed to be 0
;		* if size smaller than Xsrc, first element gets replicated
;	Sarea	[INPUT] background scaling factor in the S band
;	Marea	[INPUT] background scaling factor in the M band
;	Harea	[INPUT] background scaling factor in the H band
;		* (background region area)/(source region area)
;		* if not given, assumed to be 1
;		* can also include differences in exposure time into
;		  the ratio, in the same manner as geometric area
;		* if size smaller than Xsrc, first element gets replicated
;	Seff	[INPUT] effective area in S band
;	Meff	[INPUT] effective area in M band
;	Heff	[INPUT] effective area in H band
;		* if none are set, all are assumed to be 1,
;		  else if one is set, all are assumed to be equal to that one,
;		  else if two are set and unequal, third is assumed to be 1
;		* can also be the effective area relative to some
;		  special point on the detector (e.g., aimpoint)
;		  or even some specific detector (e.g., ACIS-I v/s ACIS-S)
;		* if size smaller than Xsrc, first element gets replicated
;	Sidx	[INPUT] index of prior on S (range = 0+)
;	Midx	[INPUT] index of prior on M (range = 0+)
;	Hidx	[INPUT] index of prior on H (range = 0+)
;		* if none are set, all are assumed to be 0.5,
;		  else if one is set, all are assumed to be equal to that one
;		  else if two are set and are unequal, third is assumed to be 0.5
;		* if size smaller than Xsrc, first element gets replicated
;		* similar to AGAMMA of PPD_SRC()
;	Sscl	[INPUT] scaling index of prior on Ssrc
;	Mscl	[INPUT] scaling index of prior on Msrc
;	Hscl	[INPUT] scaling index of prior on Hsrc
;		* if none are set, all are assumed to be 0
;		  else if one is set, all are assumed to be equal to that one
;		  else if two are set and are unequal, third is assumed to be 0
;		* if size smaller than Xsrc, first element gets replicated
;		* similar to BGAMMA of PPD_SRC()
;	Stable	[INPUT] filename containing a tabulated prior for Ssrc
;	Mtable	[INPUT] filename containing a tabulated prior for Msrc
;	Htable	[INPUT] filename containing a tabulated prior for Hsrc
;		* the table prior must be an ascii file with the following format:
;		  line 1: number of entries, say NLIN
;		  line 2: labels for the columns, ignored
;		  lines 3..NLIN+2: two whitespace separated columns of numbers,
;	            with each row containing the source intensity and the posterior
;	            density, in that order
;		* the default filenames are "./tblprior_{soft|med|hard}.txt"
;		* the default filenames are used iff Stable, Mtable, and Htable are set
;		  but are not found
;		* WARNING: if regex is used in the filename specification, only the
;		  first file from the list will be used.  furthermore, if specified,
;		  the table priors are applied to _all_ SSRC, MSRC, and HSRC
;	post	[INPUT] if set, suggests the values of (Sidx,Sscl), (Midx,Mscl),
;		and (Hidx,Hscl) going forward, i.e., what you should set the
;		priors to in your next calculation for the same source -- the
;		suggested values are stored in the output structure
;	level	[INPUT] percentage confidence level at which to report error
;		(default = 68)
;	details	[INPUT] compute various ratios (true/false)?
;		(default = true)
;	nsim	[INPUT] number of draws if algo=gibbs (default=10000)
;	nburnin	[INPUT] number of burn-in draws if algo=gibbs
;		(default=5000 or NSIM/2, whichever is smaller)
;	outputF	[INPUT] root of filename in which to place output
;		(default = 'none')
;		* output will be placed in the files OUTPUTF.txt and OUTPUTF_draws.txt
;		* NOTE: if OUTPUTF='none', then MC draws will be in BEHR_draws.txt
;	BEHRdir	[INPUT] full path to directory where BEHR executable resides
;		(default = '/data/fubar/kashyap/AstroStat/BEHR')
;	verbose	[INPUT] controls chatter
;	_extra	[JUNK] here only to prevent crashing the program
;
;requirements
;	uses subroutines HIPD_INTERVAL() and MODALPOINT()
;	BEHR executable must be installed in BEHRDIR
;	BEHR should be executable under the shell via SPAWN
;	BEHR output assumed to be compatible with 12-19-2005 version
;
;side-effects
;	potentially creates numerous ascii files in $cwd or `basedir OUTPUTF`
;
;history
;	vinay kashyap (Mar07; based on behr_hug.pro)
;	bug correction with NaNs not being caught in some cases
;	  (VK; Mar07)
;	added keywords Stable,Mtable,Htable (VK; Feb08)
;
;etymology
;	getting color-color diagrams using BEHR
;	(that's my story and I'm sticking to it)
;-

vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1

;	where is it?
if not keyword_set(BEHRdir) then begin
  stdout='BEHR: Command not found.'
  if !version.os_family eq 'unix' then spawn,'which BEHR',stdout
  if strpos(strlowcase(stdout[0]),'command not found',0) ge 0 then $
	BEHRdir='/data/fubar/kashyap/AstroStat/BEHR' else $
	BEHRdir=stdout[0]
endif
BEHRexe=filepath('BEHR',root_dir=BEHRdir)
if vv ge 5 then message,'Using : '+BEHRexe,/informational

;	usage
ok='ok'
np=n_params() & nss=n_elements(Ssrc) & nms=n_elements(Msrc) & nhs=n_elements(Hsrc)
if np lt 3 then ok='Insufficient parameters'
;
spawn,BEHRexe,usage
if strpos(strlowcase(usage[0]),'command not found',0) ge 0 then $
	ok='BEHR: command not found'
;
if nss gt 0 then cmd=BEHRexe+' Ssrc='+strtrim(Ssrc[0],2) else $
  cmd=BEHRexe+' Ssrc=0'
spawn,cmd,help

if ok ne 'ok' then begin
  print,'Usage: behr3=colbehr(Ssrc,Msrc,Hsrc, Sbkg=Sbkg,Mbkg=Mbkg,Hbkg=Hbkg,$
  print,'       Sarea=Sarea,Marea=Marea,Harea=Harea, Seff=Seff,Meff=Meff,Heff=Heff,$
  print,'       Sidx=Sidx,Midx=Midx,Hidx=Hidx, Sscl=Sscl,Mscl=Mscl,Hscl=Hscl,$
  print,'       Stable=Stable,Mtable=Mtable,Htable=Htable, /post,level=level,$'
  print,'       nbin=nbins,outputf=outputf,BEHRdir=BEHRdir,verbose=verbose)'
  print,''
  print,'  an IDL wrapper to Bayesian Estimate of Hardness Ratios (BEHR)'
  print,'  returns a structure containing the relevant outputs neatly'
  print,'  summarized into fields for IDL consumption'
  print,''
  print,'  REFERENCES:'
  print,'  Park, T., Kashyap, V.L., Siemiginowska, A., van Dyk, D., Zezas, A.,
  print,'    Heinke, C., and Wargelin, B., 2006, ApJ, 652, 610'
  print,'    http://hea-www.harvard.edu/AstroStat/BEHR/'
  print,''
  if np ne 0 then message,ok,/informational
  print,''
  if vv ne 0 then begin
    if mp eq 0 then for i=0,n_elements(usage)-1 do print,usage[i] else $
  	for i=0,n_elements(help)-1 do print,help[i]
  endif
  return,-1L
endif

;	how many are input?
ndat=nss > nms > nhs

;	set up the input parameters
;	source counts
if nss eq 0 then xSsrc=0 else xSsrc=Ssrc
if nms eq 0 then xMsrc=0 else xMsrc=Msrc
if nhs eq 0 then xHsrc=0 else xHsrc=Hsrc
ss=lonarr(ndat)+xSsrc[0] & if nss gt 0 then ss[0L:nss-1L]=xSsrc[*]
ms=lonarr(ndat)+xMsrc[0] & if nms gt 0 then ms[0L:nms-1L]=xMsrc[*]
hs=lonarr(ndat)+xHsrc[0] & if nhs gt 0 then hs[0L:nhs-1L]=xHsrc[*]
;
;	background counts
nsb=n_elements(Sbkg) & if nsb eq 0 then xSbkg=0 else xSbkg=Sbkg
nmb=n_elements(Mbkg) & if nmb eq 0 then xMbkg=0 else xMbkg=Mbkg
nhb=n_elements(Hbkg) & if nhb eq 0 then xHbkg=0 else xHbkg=Hbkg
sb=lonarr(ndat)+xSbkg[0] & if nsb gt 0 then sb[0L:nsb-1L]=xSbkg[*]
mb=lonarr(ndat)+xMbkg[0] & if nmb gt 0 then mb[0L:nmb-1L]=xMbkg[*]
hb=lonarr(ndat)+xHbkg[0] & if nhb gt 0 then hb[0L:nhb-1L]=xHbkg[*]
;
;	background area scaling factors
nsa=n_elements(Sarea) & if nsa eq 0 then xSarea=1. else xSarea=Sarea
nma=n_elements(Marea) & if nma eq 0 then xMarea=1. else xMarea=Marea
nha=n_elements(Harea) & if nha eq 0 then xHarea=1. else xHarea=Harea
sa=dblarr(ndat)+xSarea[0] & if nsa gt 0 then sa[0L:nsa-1L]=xSarea[*]
ma=dblarr(ndat)+xMarea[0] & if nma gt 0 then ma[0L:nma-1L]=xMarea[*]
ha=dblarr(ndat)+xHarea[0] & if nha gt 0 then ha[0L:nha-1L]=xHarea[*]
;
;	effective area
nse=n_elements(Seff) & if nse eq 0 then xSeff=1.0 else xSeff=Seff
nme=n_elements(Meff) & if nme eq 0 then xMeff=1.0 else xMeff=Meff
nhe=n_elements(Heff) & if nhe eq 0 then xHeff=1.0 else xHeff=Heff
if nse eq 0 and nme eq 0 and nhe ne 0 then begin & xSeff=xHeff & xMeff=xHeff & endif
if nse eq 0 and nme ne 0 and nhe eq 0 then begin & xSeff=xMeff & xHeff=xMeff & endif
if nse ne 0 and nme eq 0 and nhe eq 0 then begin & xMeff=xSeff & xHeff=xSeff & endif
se=dblarr(ndat)+xSeff[0] & if nse gt 0 then se[0L:nse-1L]=xSeff[*]
me=dblarr(ndat)+xMeff[0] & if nme gt 0 then me[0L:nme-1L]=xMeff[*]
he=dblarr(ndat)+xHeff[0] & if nhe gt 0 then he[0L:nhe-1L]=xHeff[*]
;
;	index of prior
nsi=n_elements(Sidx) & if nsi eq 0 then xSidx=0.5 else xSidx=Sidx
nmi=n_elements(Midx) & if nmi eq 0 then xMidx=0.5 else xMidx=Midx
nhi=n_elements(Hidx) & if nhi eq 0 then xHidx=0.5 else xHidx=Hidx
if nsi eq 0 and nmi eq 0 and nhi ne 0 then begin & xSidx=xHidx & xMidx=xHidx & endif
if nsi eq 0 and nmi ne 0 and nhi eq 0 then begin & xSidx=xMidx & xHidx=xMidx & endif
if nsi ne 0 and nmi eq 0 and nhi eq 0 then begin & xMidx=xSidx & xHidx=xSidx & endif
si=fltarr(ndat)+xSidx[0] & if nsi gt 0 then si[0L:nsi-1L]=xSidx[*]
mi=fltarr(ndat)+xMidx[0] & if nmi gt 0 then mi[0L:nmi-1L]=xMidx[*]
hi=fltarr(ndat)+xHidx[0] & if nhi gt 0 then hi[0L:nhi-1L]=xHidx[*]
;
;	scale index of prior
nsc=n_elements(Sscl) & if nsc eq 0 then xSscl=0. else xSscl=Sscl
nmc=n_elements(Mscl) & if nmc eq 0 then xMscl=0. else xMscl=Mscl
nhc=n_elements(Hscl) & if nhc eq 0 then xHscl=0. else xHscl=Hscl
if nsc eq 0 and nmc eq 0 and nhc ne 0 then begin & xSscl=xHscl & xMscl=xHscl & endif
if nsc eq 0 and nmc ne 0 and nhc eq 0 then begin & xSscl=xMscl & xHscl=xMscl & endif
if nsc ne 0 and nmc eq 0 and nhc eq 0 then begin & xMscl=xSscl & xHscl=xSscl & endif
sc=fltarr(ndat)+xSscl[0] & if nsc gt 0 then sc[0L:nsc-1L]=xSscl[*]
mc=fltarr(ndat)+xMscl[0] & if nmc gt 0 then mc[0L:nmc-1L]=xMscl[*]
hc=fltarr(ndat)+xHscl[0] & if nhc gt 0 then hc[0L:nhc-1L]=xHscl[*]
;
nlev=n_elements(level) & if nlev eq 0 then xlevel=68. else xlevel=level
clev=fltarr(ndat)+xlevel[0] & if nlev gt 0 then clev[0L:nlev-1L]=xlevel[*]
;
msim=n_elements(nsim) & if msim eq 0 then nnsim=10000L else nnsim=nsim
numsim=lonarr(ndat)+nnsim[0] & if msim gt 0 then numsim[0L:msim-1L]=nnsim[*]
;
mburn=n_elements(nburnin) & if mburn eq 0 then nnburnin=(5000L < numsim/2) else nnburnin=nburnin
numburn=lonarr(ndat)+nnburnin[0] & if mburn gt 0 then numburn[0L:mburn-1L]=nnburnin[*]
;
nfmt=strtrim(fix(alog10(ndat)+1),2) & fmt='(i'+nfmt+'.'+nfmt+')'
nout=n_elements(outputf) & if nout eq 0 then outputf='none'
if ndat gt 1 then filroot=outputf[0]+'_'+string(lindgen(ndat)+1L,fmt) else $
	filroot=outputf[0]
if ndat eq 1 then if nout gt 1 then filroot[0L:nout-1L]=outputf[*]
;
mPost=n_elements(post) & if mPost eq 0 then iPost=0 else iPost=Post
outPost=bytarr(ndat)+iPost[0] & if mPost gt 0 then outPost[0L:mPost-1L]=iPost[*]
soutPost=strarr(ndat)+'false' & oo=where(outPost ne 0,moo) & if moo gt 0 then soutPost[oo]='true'
;
sinthpd=strarr(ndat)+'true'
;
;	table priors
if keyword_set(Stable) then begin
  stbl='tblprior_soft.txt'
  sfil=findfile(string(Stable),count=nsfil)
  if nsfil ne 0 then stbl=sfil[0]
endif
if keyword_set(Mtable) then begin
  Mtbl='tblprior_med.txt'
  Mfil=findfile(string(Mtable),count=nMfil)
  if nMfil ne 0 then Mtbl=Mfil[0]
endif
if keyword_set(Htable) then begin
  Htbl='tblprior_hard.txt'
  Hfil=findfile(string(Htable),count=nHfil)
  if nHfil ne 0 then Htbl=Hfil[0]
endif

;	define the output
dat=fltarr(ndat) & nat=dat+!values.F_NAN
Sstr=create_struct('mode',nat,'mean',nat,'median',nat,'lowerbound',nat,'upperbound',nat)
Mstr=Sstr & Hstr=Sstr					;S, M, H,
lSstr=Sstr & lMstr=Sstr & lHstr=Sstr			;logS, logM, logH,
SmMstr=Sstr & MmHstr=Sstr & SmHstr=Sstr			;SmM=S-M, MmH=M-H, SmH=S-H,
SpMstr=Sstr & MpHstr=Sstr & SpHstr=Sstr & Tstr=Sstr	;SpM=S+M, MpH=M+H, SpH=S+H, T=S+M+H
R1str=Sstr & R2str=Sstr & R3str=Sstr			;R1=S/M, R2=M/H, R3=S/H,
C1str=Sstr & C2str=Sstr & C3str=Sstr			;C1=C_SM, C2=C_MH, C3=C_SH,
HR1str=Sstr & HR2str=Sstr & HR3str=Sstr			;HR1=(S-M)/(S+M), HR2=(M-H)/(M+H), HR3=(S-H)/(S+H),
HRAstr=Sstr & HRBstr=Sstr & HRCstr=Sstr			;HRA=(S-M)/(S+M+H), HRB=(M-H)/(S+M+H), HRC=(S-H)/(S+M+H)}
;
postlamSidx=si & postlamMidx=mi & postlamHidx=hi
postlamSscl=dat & postlamMscl=dat & postlamHscl=dat

;	call BEHR for each data point
for i=0L,ndat-1L do begin
  if vv ne 0 and vv le 4 then print,strtrim(i,2),format='($,a10)'
  algo_method='gibbs'

  if ss[i]-sb[i]/sa[i] lt 15 or ms[i]-mb[i]/ma[i] lt 15 or hs[i]-hb[i]/ha[i] lt 15 then begin
    if ss[i]-sb[i]/sa[i] gt 100 or ms[i]-mb[i]/ma[i] gt 100 or hs[i]-hb[i]/ha[i] gt 100 then begin
      if numsim[i] lt 50000 then message,'WARNING: '+$
	'consider raising NSIM; too many counts in one band and too few in the other',/informational
    endif
  endif

  cmd1=BEHRexe+$
	' softsrc='+strtrim(ss[i],2)+$
	' hardsrc='+strtrim(ms[i],2)+$
	' softbkg='+strtrim(sb[i],2)+$
	' hardbkg='+strtrim(mb[i],2)+$
	' softarea='+strtrim(sa[i],2)+$
	' hardarea='+strtrim(ma[i],2)+$
	' softeff='+strtrim(se[i],2)+$
	' hardeff='+strtrim(me[i],2)+$
	' softidx='+strtrim(si[i],2)+$
	' hardidx='+strtrim(mi[i],2)+$
	' softscl='+strtrim(sc[i],2)+$
	' hardscl='+strtrim(mc[i],2)
  if keyword_set(Stable) then cmd1=cmd1+' softtbl='+stbl
  if keyword_set(Mtable) then cmd1=cmd1+' hardtbl='+mtbl
  cmd1=cmd1+$
	' post='+soutPost[i]+$
	' level='+strtrim(clev[i],2)+$
	' algo=gibbs'+$
	' details=true'+$
	' nsim='+strtrim(numsim[i],2)+$
	' nburnin='+strtrim(numburn[i],2)+$
	' HPD=true'+$
	' output='+filroot[i]+$
	' outputMC=true'
  if vv gt 4 then print,cmd1
  ;
  if numburn[i] ge numsim[i] then begin
    message,'NBURNIN must be less than NSIM; skipping this command:',$
    /informational
    print,cmd1
  endif else spawn,cmd1,stdout

  MCdrawsfil='BEHR_draws.txt'
  if strpos(strlowcase(filroot[i]),'none') lt 0 then MCdrawsfil=filroot[i]+'_draws.txt'
  var=fltarr(2,numsim[i]-numburn[i])
  openr,umc,MCdrawsfil,/get_lun
  readf,umc,var & Sdraw=reform(var[0,*]) & Mdraw=reform(var[1,*])
  close,umc & free_lun,umc

  nn=n_elements(stdout) & j=0L & k=-1
  while k lt 0 do begin
    if strpos(stdout[j],'Mode',0) ge 0 and $
       strpos(stdout[j],'Mean',0) ge 0 and $
       strpos(stdout[j],'Median',0) ge 0 and $
       strpos(stdout[j],'Lower Bound',0) ge 0 and $
       strpos(stdout[j],'Upper Bound',0) ge 0 then k=j
    if k gt 0 and vv ge 10 then print,stdout[j]
    j=j+1L
    if j eq nn then message,'BUG?  BEHR returns nothing'
  endwhile
  for j=k,nn-1L do begin
    if strpos(stdout[j],'p(lamS|data)',0) ge 0 then begin
      if vv ge 10 then print,stdout[j]
      i0=strpos(stdout[j],'(',strpos(stdout[j],'p(lamS|data)',0)+strlen('p(lamS|data)'))
      i1=strpos(stdout[j],')',i0+1)
      cc=strmid(stdout[j],i0+1,i1-i0-1)
      ccc=strsplit(cc,',',/extract)
      postlamSidx[i]=float(ccc[0])
      postlamSscl[i]=float(ccc[1])
    endif
    if strpos(stdout[j],'p(lamH|data)',0) ge 0 then begin
      if vv ge 10 then print,stdout[j]
      i0=strpos(stdout[j],'(',strpos(stdout[j],'p(lamH|data)',0)+strlen('p(lamH|data)'))
      i1=strpos(stdout[j],')',i0+1)
      cc=strmid(stdout[j],i0+1,i1-i0-1)
      ccc=strsplit(cc,',',/extract)
      postlamMidx[i]=float(ccc[0])
      postlamMscl[i]=float(ccc[1])
    endif
  endfor

  cmd2=BEHRexe+$
	' softsrc='+strtrim(ss[i],2)+$
	' hardsrc='+strtrim(hs[i],2)+$
	' softbkg='+strtrim(sb[i],2)+$
	' hardbkg='+strtrim(hb[i],2)+$
	' softarea='+strtrim(sa[i],2)+$
	' hardarea='+strtrim(ha[i],2)+$
	' softeff='+strtrim(se[i],2)+$
	' hardeff='+strtrim(he[i],2)+$
	' softidx='+strtrim(si[i],2)+$
	' hardidx='+strtrim(hi[i],2)+$
	' softscl='+strtrim(sc[i],2)+$
	' hardscl='+strtrim(hc[i],2)
  if keyword_set(Stable) then cmd2=cmd2+' softtbl='+mtbl
  if keyword_set(Mtable) then cmd2=cmd2+' hardtbl='+htbl
  cmd2=cmd2+$
	' post='+soutPost[i]+$
	' level='+strtrim(clev[i],2)+$
	' algo=gibbs'+$
	' details=true'+$
	' nsim='+strtrim(numsim[i],2)+$
	' nburnin='+strtrim(numburn[i],2)+$
	' HPD=true'+$
	' output='+filroot[i]+$
	' outputMC=true'
  if vv gt 4 then print,cmd2
  ;
  if numburn[i] ge numsim[i] then begin
    message,'NBURNIN must be less than NSIM; skipping this command:',$
    /informational
    print,cmd2
  endif else spawn,cmd2,stdout

  MCdrawsfil='BEHR_draws.txt'
  if strpos(strlowcase(filroot[i]),'none') lt 0 then MCdrawsfil=filroot[i]+'_draws.txt'
  var=fltarr(2,numsim[i]-numburn[i])
  openr,umc,MCdrawsfil,/get_lun
  readf,umc,var & Sdraw=reform(var[0,*]) & Hdraw=reform(var[1,*])
  close,umc & free_lun,umc

  nn=n_elements(stdout) & j=0L & k=-1
  while k lt 0 do begin
    if strpos(stdout[j],'Mode',0) ge 0 and $
       strpos(stdout[j],'Mean',0) ge 0 and $
       strpos(stdout[j],'Median',0) ge 0 and $
       strpos(stdout[j],'Lower Bound',0) ge 0 and $
       strpos(stdout[j],'Upper Bound',0) ge 0 then k=j
    if k gt 0 and vv ge 10 then print,stdout[j]
    j=j+1L
    if j eq nn then message,'BUG?  BEHR returns nothing'
  endwhile
  for j=k,nn-1L do begin
    if strpos(stdout[j],'p(lamH|data)',0) ge 0 then begin
      if vv ge 10 then print,stdout[j]
      i0=strpos(stdout[j],'(',strpos(stdout[j],'p(lamH|data)',0)+strlen('p(lamH|data)'))
      i1=strpos(stdout[j],')',i0+1)
      cc=strmid(stdout[j],i0+1,i1-i0-1)
      ccc=strsplit(cc,',',/extract)
      postlamHidx[i]=float(ccc[0])
      postlamHscl[i]=float(ccc[1])
    endif
  endfor

  y=Sdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  Sstr.mode[i]=fmode & Sstr.mean[i]=mean(y) & Sstr.median[i]=median(y) & Sstr.lowerbound[i]=z[0] & Sstr.upperbound[i]=z[1]
  y=Mdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  Mstr.mode[i]=fmode & Mstr.mean[i]=mean(y) & Mstr.median[i]=median(y) & Mstr.lowerbound[i]=z[0] & Mstr.upperbound[i]=z[1]
  y=Hdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  Hstr.mode[i]=fmode & Hstr.mean[i]=mean(y) & Hstr.median[i]=median(y) & Hstr.lowerbound[i]=z[0] & Hstr.upperbound[i]=z[1]

  y=alog10(Sdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    lSstr.mode[i]=fmode & lSstr.mean[i]=mean(y[os]) & lSstr.median[i]=median(y[os]) & lSstr.lowerbound[i]=z[0] & lSstr.upperbound[i]=z[1]
  endif
  y=alog10(Mdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    lMstr.mode[i]=fmode & lMstr.mean[i]=mean(y[os]) & lMstr.median[i]=median(y[os]) & lMstr.lowerbound[i]=z[0] & lMstr.upperbound[i]=z[1]
  endif
  y=alog10(Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    lHstr.mode[i]=fmode & lHstr.mean[i]=mean(y[os]) & lHstr.median[i]=median(y[os]) & lHstr.lowerbound[i]=z[0] & lHstr.upperbound[i]=z[1]
  endif

  y=Sdraw-Mdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  SmMstr.mode[i]=fmode & SmMstr.mean[i]=mean(y) & SmMstr.median[i]=median(y) & SmMstr.lowerbound[i]=z[0] & SmMstr.upperbound[i]=z[1]
  y=Mdraw-Hdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  MmHstr.mode[i]=fmode & MmHstr.mean[i]=mean(y) & MmHstr.median[i]=median(y) & MmHstr.lowerbound[i]=z[0] & MmHstr.upperbound[i]=z[1]
  y=Sdraw-Hdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  SmHstr.mode[i]=fmode & SmHstr.mean[i]=mean(y) & SmHstr.median[i]=median(y) & SmHstr.lowerbound[i]=z[0] & SmHstr.upperbound[i]=z[1]

  y=Sdraw+Mdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  SpMstr.mode[i]=fmode & SpMstr.mean[i]=mean(y) & SpMstr.median[i]=median(y) & SpMstr.lowerbound[i]=z[0] & SpMstr.upperbound[i]=z[1]
  y=Mdraw+Hdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  MpHstr.mode[i]=fmode & MpHstr.mean[i]=mean(y) & MpHstr.median[i]=median(y) & MpHstr.lowerbound[i]=z[0] & MpHstr.upperbound[i]=z[1]
  y=Sdraw+Hdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  SpHstr.mode[i]=fmode & SpHstr.mean[i]=mean(y) & SpHstr.median[i]=median(y) & SpHstr.lowerbound[i]=z[0] & SpHstr.upperbound[i]=z[1]
  y=Sdraw+Mdraw+Hdraw & z=hipd_interval(y,/fsample,fmode=fmode,clev=clev[i])
  Tstr.mode[i]=fmode & Tstr.mean[i]=mean(y) & Tstr.median[i]=median(y) & Tstr.lowerbound[i]=z[0] & Tstr.upperbound[i]=z[1]

  y=float(Sdraw)/float(Mdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    R1str.mode[i]=fmode & R1str.mean[i]=mean(y[os]) & R1str.median[i]=median(y[os]) & R1str.lowerbound[i]=z[0] & R1str.upperbound[i]=z[1]
  endif
  y=float(Mdraw)/float(Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    R2str.mode[i]=fmode & R2str.mean[i]=mean(y[os]) & R2str.median[i]=median(y[os]) & R2str.lowerbound[i]=z[0] & R2str.upperbound[i]=z[1]
  endif
  y=float(Sdraw)/float(Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    R3str.mode[i]=fmode & R3str.mean[i]=mean(y[os]) & R3str.median[i]=median(y[os]) & R3str.lowerbound[i]=z[0] & R3str.upperbound[i]=z[1]
  endif

  y=alog10(Sdraw)-alog10(Mdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    C1str.mode[i]=fmode & C1str.mean[i]=mean(y[os]) & C1str.median[i]=median(y[os]) & C1str.lowerbound[i]=z[0] & C1str.upperbound[i]=z[1]
  endif
  y=alog10(Mdraw)-alog10(Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    C2str.mode[i]=fmode & C2str.mean[i]=mean(y[os]) & C2str.median[i]=median(y[os]) & C2str.lowerbound[i]=z[0] & C2str.upperbound[i]=z[1]
  endif
  y=alog10(Sdraw)-alog10(Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    C3str.mode[i]=fmode & C3str.mean[i]=mean(y[os]) & C3str.median[i]=median(y[os]) & C3str.lowerbound[i]=z[0] & C3str.upperbound[i]=z[1]
  endif

  y=(Sdraw-Mdraw)/(Sdraw+Mdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    HR1str.mode[i]=fmode & HR1str.mean[i]=mean(y[os]) & HR1str.median[i]=median(y[os]) & HR1str.lowerbound[i]=z[0] & HR1str.upperbound[i]=z[1]
  endif
  y=(Mdraw-Hdraw)/(Mdraw+Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    HR2str.mode[i]=fmode & HR2str.mean[i]=mean(y[os]) & HR2str.median[i]=median(y[os]) & HR2str.lowerbound[i]=z[0] & HR2str.upperbound[i]=z[1]
  endif
  y=(Sdraw-Hdraw)/(Sdraw+Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    HR3str.mode[i]=fmode & HR3str.mean[i]=mean(y[os]) & HR3str.median[i]=median(y[os]) & HR3str.lowerbound[i]=z[0] & HR3str.upperbound[i]=z[1]
  endif

  y=(Sdraw-Mdraw)/(Sdraw+Mdraw+Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    HRAstr.mode[i]=fmode & HRAstr.mean[i]=mean(y[os]) & HRAstr.median[i]=median(y[os]) & HRAstr.lowerbound[i]=z[0] & HRAstr.upperbound[i]=z[1]
  endif
  y=(Mdraw-Hdraw)/(Sdraw+Mdraw+Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    HRBstr.mode[i]=fmode & HRBstr.mean[i]=mean(y[os]) & HRBstr.median[i]=median(y[os]) & HRBstr.lowerbound[i]=z[0] & HRBstr.upperbound[i]=z[1]
  endif
  y=(Sdraw-Hdraw)/(Sdraw+Mdraw+Hdraw) & os=where(finite(y) ne 0,mos)
  if mos gt 2 then begin
    z=hipd_interval(y[os],/fsample,fmode=fmode,clev=clev[i])
    HRCstr.mode[i]=fmode & HRCstr.mean[i]=mean(y[os]) & HRCstr.median[i]=median(y[os]) & HRCstr.lowerbound[i]=z[0] & HRCstr.upperbound[i]=z[1]
  endif

endfor

;	define the output
BEHR3=create_struct('softsrc',ss,'mediumsrc',ms,'hardsrc',hs,$
	'softbkg',sb,'mediumbkg',mb,'hardbkg',hb,$
	'softarea',sa,'mediumarea',ma,'hardarea',ha,$
	'softeff',se,'mediumeff',me,'hardeff',he,$
	'softidx',si,'mediumidx',mi,'hardidx',hi,$
	'softscl',sc,'mediumscl',mc,'hardscl',hc,$
	'level',clev,'algo','gibbs','nsim',numsim,'nburnin',numburn,$
	'nbins',0,'HPD',sinthpd,'filroot',filroot,$
	'S',Sstr,'M',Mstr,'H',Hstr,'lS',lSstr,'lM',lMstr,'lH',lHstr,$
	'SmM',SmMstr,'MmH',MmHstr,'SmH',SmHstr,'SpM',SpMstr,'MpH',MpHstr,'SpH',SpHstr,'T',Tstr,$
	'R1',R1str,'R2',R2str,'R3',R3str,'C1',C1str,'C2',C2str,'C3',C3str,$
	'HR1',HR1str,'HR2',HR2str,'HR3',HR3str,'HRA',HRAstr,'HRB',HRBstr,'HRC',HRCstr,$
	'postlamS_idx',postlamSidx,'postlamS_scl',postlamSscl,$
	'postlamM_idx',postlamMidx,'postlamM_scl',postlamMscl,$
	'postlamH_idx',postlamHidx,'postlamH_scl',postlamHscl,$
	'help', 'S, M, H ; lS=log10(S), lM=log10(M), lH=log10(H) ; '+$
		'SpM=S+M, MpH=M+H, SpH=S+H, T=S+M+H ; '+$
		'SmM=S-M, MmH=M-H, SmH=S-H ; '+$
		'R1=S/M, R2=M/H, R3=S/H ; C1=C_SM, C2=C_MH, C3=C_SH ; '+$
		'HR1=(S-M)/(S+M), HR2=(M-H)/(M+H), HR3=(S-H)/(S+H) ; '+$
		'HRA=(S-M)/(S+M+H), HRB=(M-H)/(S+M+H), HRC=(S-H)/(S+M+H)')

return,BEHR3
end