Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/AstroStat/BEHR/behr_test.pro Дата изменения: Thu May 4 23:25:45 2006 Дата индексирования: Tue Oct 2 01:09:53 2012 Кодировка: Поисковые слова: п п п п п п п п п п р п

;+
;behr_test.pro
; an example program that reads in data from an ascii file containing
; source and background counts and backscale factors, calls the IDL
; wrapper routine to BEHR to generate hardness ratios and plots them.
;
;usage
; infile='behr_test.inp'
; BEHRdir='./' ;REQUIRED full path to location of BEHR executable
; psfile='behr_test.ps' ;if you want the output plots in postscript
; outfile='behr_test.asc' ;if you want the output in an ascii table format
; algo='quad' ;algo='gibbs' for using the Gibbs sampler,
; ;algo='quad' for Quadrature
; ;algo='auto=20' to choose automatically
; .run behr_test
;
; NOTE: additional named input keyword variables for BEHR_HUG()
; may also be defined prior to running this script. Type
; help,behr_hug()
; to obtain a list of accepted keywords. See the documentation
; for BEHR_HUG or the BEHR Tutorial Guide for a description of
; the variables.
;
;restrictions
; it is assumed that the input file contains data as an ASCII table,
; with columns
; -- soft-band counts in source region
; -- hard-band counts in source region
; -- soft-band counts in background region
; -- hard-band counts in background region
; -- ratio of background-to-source area for soft band
; -- ratio of background-to-source area for hard band
;
;requires
; PINTofALE function BEHR_HUG
; IDLAstro routine READCOL (and associated)
;
;history
; vinay kashyap (2005-dec-17)
; added option to print output to ascii file (VK; 2006-may-04)
;-

; check if input file is defined
if n_elements(infile) eq 0 then begin
message,'INFILE is not defined; assuming default of "behr_test.inp"',/informational
infile='behr_test.inp'
endif

; check where the BEHR executable is
if n_elements(BEHRdir) eq 0 then begin
message,'BEHRdir is not defined; please specify directory where executable is'
endif else begin
fil=findfile(BEHRdir[0]+'/BEHR',count=nfil)
if nfil eq 0 then begin
message,BEHRdir[0]+': does not contain BEHR executable'
endif
endelse

; check if input file exists
filename=findfile(infile[0],count=nfil)
if nfil eq 0 then message,infile[0]+': file not found'

; some initializations
if n_elements(verbose) eq 0 then verbose=1

; read from input file
readcol,filename[0],softsrc,hardsrc,softbkg,hardbkg,softarea,hardarea

; call the wrapper
ok='ok'
if n_tags(behr) eq 0 then ok='must run behr_hug' else $
if (tag_names(behr))[0] ne 'SOFTSRC' then ok='structure not in right format' else $
if n_elements(behr.SOFTSRC) ne n_elements(softsrc) then ok='new file read in'

if ok ne 'ok' then behr=behr_hug(softsrc,hardsrc,softbkg,hardbkg,softarea,hardarea,$
softeff=softeff,hardeff=hardeff,softidx=softidx,hardidx=hardidx,$
softscl=softscl,hardscl=hardscl,post=post,$
level=level,algo=algo,details=details,nsim=nsim,nburnin=nburnin,$
nbin=nbins,outputf=outputf,outputR=outputR,outputHR=outputHR,$
outputC=outputC,outputMC=outputMC,outputPr=outputPr,BEHRdir=BEHRdir[0],$
verbose=verbose)

; make plots

; set up display
pmulti=!p.multi & if !d.NAME eq 'X' then !p.multi=[0,2,2]
if keyword_set(psfile) then begin
!p.multi=pmulti
setplot=!d.NAME & set_plot,'ps' & device,file=psfile[0],/landscape,/color
endif else begin
message,'output will be directed to screen, not to PSFILE',/informational
if !d.N_COLORS gt 256 then device,decomposed=0
endelse
loadct,3

; logS v/s logH
plot,behr.S.mode,behr.H.mode,psym=4,symsize=2,charsize=1.4,thick=2,$
xtitle='S',ytitle='H',/xlog,/ylog,xr=[1,3e3],yr=[1,1e3]
for i=0,n_elements(softsrc)-1L do $
oplot,[behr.S.lowerbound[i],behr.S.upperbound[i]]>1e-10,behr.H.mode[i]*[1,1],color=175
for i=0,n_elements(softsrc)-1L do $
oplot,behr.S.mode[i]*[1,1],[behr.H.lowerbound[i],behr.H.upperbound[i]]>1e-10,color=175

; R v/s H+S
plot,behr.R.mode,behr.HpS.mode,psym=4,symsize=2,charsize=1.4,thick=2,$
xtitle='R=S/H',ytitle='H+S',/ylog,xr=[0,10],yr=[1,1e3]
for i=0,n_elements(softsrc)-1L do $
oplot,[behr.R.lowerbound[i],behr.R.upperbound[i]]>1e-10,behr.HpS.mode[i]*[1,1],color=175
for i=0,n_elements(softsrc)-1L do $
oplot,behr.R.mode[i]*[1,1],[behr.HpS.lowerbound[i],behr.HpS.upperbound[i]]>1e-10,color=175

; C v/s H
plot,behr.C.mode,behr.H.mode,psym=4,symsize=2,charsize=1.4,thick=2,$
xtitle='C=log!D10!N(S/H)',ytitle='H',/ylog,xr=[-1.5,1.5],yr=[1,1e3]
for i=0,n_elements(softsrc)-1L do $
oplot,[behr.C.lowerbound[i],behr.C.upperbound[i]],behr.H.mode[i]*[1,1],color=175
for i=0,n_elements(softsrc)-1L do $
oplot,behr.C.mode[i]*[1,1],[behr.H.lowerbound[i],behr.H.upperbound[i]]>1e-10,color=175

; HR v/s H+S
plot,behr.HR.mean,behr.HpS.mode,psym=4,symsize=2,charsize=1.4,thick=2,$
xtitle='HR=(H-S)/(H+S)',ytitle='H+S',/ylog,xr=[-1,1],yr=[1,1e3]
for i=0,n_elements(softsrc)-1L do $
oplot,[behr.HR.lowerbound[i],behr.HR.upperbound[i]],behr.HpS.mode[i]*[1,1],color=175
for i=0,n_elements(softsrc)-1L do $
oplot,behr.HR.mean[i]*[1,1],[behr.HpS.lowerbound[i],behr.HpS.upperbound[i]]>1e-10,color=175

!p.multi=pmulti
if keyword_set(psfile) then begin
device,/close & set_plot,setplot
endif

; print output to ascii table
if keyword_set(outfile) then begin
print,'writing S,H,R,C,HR,H+S to file '+outfile[0]+' ...'
openw,uout,outfile[0],/get_lun

printf,uout,'S S.lowerbound S.upperbound'
for i=0,n_elements(softsrc)-1L do printf,uout,(behr.S.mode)[i],$
' ',(behr.S.lowerbound)[i],' ',(behr.S.upperbound)[i]

printf,uout,''
printf,uout,'H H.lowerbound H.upperbound'
for i=0,n_elements(softsrc)-1L do printf,uout,(behr.H.mode)[i],$
' ',(behr.H.lowerbound)[i],' ',(behr.H.upperbound)[i]

printf,uout,''
printf,uout,'R R.lowerbound R.upperbound'
for i=0,n_elements(softsrc)-1L do printf,uout,(behr.R.mode)[i],$
' ',(behr.R.lowerbound)[i],' ',(behr.R.upperbound)[i]

printf,uout,''
printf,uout,'C C.lowerbound C.upperbound'
for i=0,n_elements(softsrc)-1L do printf,uout,(behr.C.mode)[i],$
' ',(behr.C.lowerbound)[i],' ',(behr.C.upperbound)[i]

printf,uout,''
printf,uout,'HR HR.lowerbound HR.upperbound'
for i=0,n_elements(softsrc)-1L do printf,uout,(behr.HR.mean)[i],$
' ',(behr.HR.lowerbound)[i],' ',(behr.HR.upperbound)[i]

printf,uout,''
printf,uout,'H+S H+S.lowerbound H+S.upperbound'
for i=0,n_elements(softsrc)-1L do printf,uout,(behr.HpS.mean)[i],$
' ',(behr.HpS.lowerbound)[i],' ',(behr.HpS.upperbound)[i]

close,uout & free_lun,uout
print,'... done.'
endif else begin
message,'Set OUTFILE to name of file to contain ascii output',/informational
endelse

save,file='behr_test.save'

end