Документ взят из кэша поисковой машины. Адрес оригинального документа : http://xmm.vilspa.esa.es/calibration/ept/flux/fluxer.20090326.pro Дата изменения: Mon Nov 3 17:14:31 2008 Дата индексирования: Sun Apr 10 14:12:01 2016 Кодировка: Поисковые слова: mercury program

;
; PRO fluxer
;
; PURPOSE: Create a plot with the statistica distirbution of fluxes
;
; INPUT:
; - "${version}/${obs_id}.txt"
;
; OUTPUT:
; - "${version}.ps"
;
; HISTORY:
;
; 1.0 - <31 Oct 2008 - Original version (MH)
; 1.1 - 31 Oct 2008 - Remove printout
; Improved comments
; Changed format of mean and st. deviation output (MG)
PRO fluxer
;
; Identify SAS versions to be analyzed for:
;
; Chandra
spawn, 'ls | egrep "^SAS"> versions1.dat'
readcol, 'versions1.dat', versions1, format='(a)', /silent
spawn, 'rm versions1.dat'
; XMM-Newton
spawn, 'ls | egrep "^xmmsas"> versions2.dat'
readcol, 'versions2.dat', versions2, format='(a)', /silent
spawn, 'rm versions2.dat'
if (n_elements(versions1) EQ 0) then begin
print,"----------XMM-Newton----------";
versions=[versions2]
endif else begin
print,"----------Chandra-------------";
versions=[versions1]
endelse
;
; Identify files with the flux ratios
;
ve=n_elements(versions)
if (n_elements(versions1) NE 0) then begin
index_chandra=where(versions EQ versions1)
endif
;
; Loop on all available softwarev ersions
;
for vi=0,ve-1 do begin
versionstring_odf=versions(vi)+"/"+"*.txt"
if (n_elements(versions1) NE 0) then begin
check_chandra=where(versions1 EQ versions(vi))
endif else begin
check_chandra=1
endelse
spawn, 'ls -1 '+versionstring_odf+' > odf.dat'
readcol, 'odf.dat', odf, format='(a)', /silent
version=versions(vi)
spawn, 'rm odf.dat'
e=n_elements(odf)
b_num=numlines(odf(0))-1
;
; Definition of instrument based vectors
;
aciss_heg_m1 =fltarr(b_num,e)
aciss_heg_p1 =fltarr(b_num,e)
aciss_leg_m1 =fltarr(b_num,e)
aciss_leg_p1 =fltarr(b_num,e)
aciss_meg_m1 =fltarr(b_num,e)
aciss_meg_p1 =fltarr(b_num,e)
combined =fltarr(b_num,e)
hrci_leg_m123 =fltarr(b_num,e)
hrci_leg_p123 =fltarr(b_num,e)
hrcs_leg_m123 =fltarr(b_num,e)
hrcs_leg_p123 =fltarr(b_num,e)
m1 =fltarr(b_num,e)
m2 =fltarr(b_num,e)
pn =fltarr(b_num,e)
r1 =fltarr(b_num,e)
r2 =fltarr(b_num,e)
;
; Read fluxes
;
readcol, odf(0), name, format='(A,)', /silent
x=findgen(n_elements(odf))
for i=0,e-1 do begin
if (check_chandra EQ 0) then begin
readcol, odf(i), dummy, $
aciss_heg_m1i,$
aciss_heg_p1i,$
aciss_leg_m1i,$
aciss_leg_p1i,$
aciss_meg_m1i,$
aciss_meg_p1i,$
combinedi,$
hrci_leg_m123i,$
hrci_leg_p123i,$
hrcs_leg_m123i,$
hrcs_leg_p123i,$
m1i,$
m2i,$
pni,$
r1i,$
r2i,$
format='(A,f,f, f)', /silent
endif else begin
readcol, odf(i), dummy, $
combinedi,$
pni,$
m1i,$
m2i,$
r1i,$
r2i,$
format='(A,f,f, f)', /silent
aciss_heg_m1i=fltarr(n_elements(pni))*0
aciss_heg_p1i=fltarr(n_elements(pni))*0
aciss_leg_m1i=fltarr(n_elements(pni))*0
aciss_leg_p1i=fltarr(n_elements(pni))*0
aciss_meg_m1i=fltarr(n_elements(pni))*0
aciss_meg_p1i=fltarr(n_elements(pni))*0
hrci_leg_m123i=fltarr(n_elements(pni))*0
hrci_leg_p123i=fltarr(n_elements(pni))*0
hrcs_leg_m123i=fltarr(n_elements(pni))*0
hrcs_leg_p123i=fltarr(n_elements(pni))*0
endelse
;
; Normalize fluxes to the value of the flux in the combined fit
;
for k=0,b_num-1 do begin
combined(k, i)=combinedi(k)
if (combinedi(k) NE 0) then aciss_heg_m1(k, i)= aciss_heg_m1i(k)/combinedi(k)
if (combinedi(k) NE 0) then aciss_heg_p1(k, i)= aciss_heg_p1i(k)/combinedi(k)
if (combinedi(k) NE 0) then aciss_leg_m1(k, i)= aciss_leg_m1i(k)/combinedi(k)
if (combinedi(k) NE 0) then aciss_leg_p1(k, i)= aciss_leg_p1i(k)/combinedi(k)
if (combinedi(k) NE 0) then aciss_meg_m1(k, i)= aciss_meg_m1i(k)/combinedi(k)
if (combinedi(k) NE 0) then aciss_meg_p1(k, i)= aciss_meg_p1i(k)/combinedi(k)
if (combinedi(k) NE 0) then combined(k, i)= combinedi(k)
if (combinedi(k) NE 0) then hrci_leg_m123(k, i)= hrci_leg_m123i(k)/combinedi(k)
if (combinedi(k) NE 0) then hrci_leg_p123(k, i)= hrci_leg_p123i(k)/combinedi(k)
if (combinedi(k) NE 0) then hrcs_leg_m123(k, i)= hrcs_leg_m123i(k)/combinedi(k)
if (combinedi(k) NE 0) then hrcs_leg_p123(k, i)= hrcs_leg_p123i(k)/combinedi(k)
if (combinedi(k) NE 0) then m1(k, i)= m1i(k)/combinedi(k)
if (combinedi(k) NE 0) then m2(k, i)= m2i(k)/combinedi(k)
if (combinedi(k) NE 0) then pn(k, i)= pni(k)/combinedi(k)
if (combinedi(k) NE 0) then r1(k, i)= r1i(k)/combinedi(k)
if (combinedi(k) NE 0) then r2(k, i)= r2i(k)/combinedi(k)
endfor
endfor
;
; Plot output file
;
filename="plot"+versions(vi)+".ps"
set_plot,'ps'
device, filename=filename,/portrait,xsize=18, ysize=28,yoffset=0,/times, /color
noticks=replicate(' ',8)
symsize=1.5
thick=10
loadct,13
;
; Loop on the energy ranges ("b" stands for "bands"?)
;
for k=1,b_num do begin
;
; Determine position of each histogram plot
;
bin=0.02
min=0.0000000001
par_num=b_num
pydynmin=0.93*(k-1)/(par_num)+0.08
pydynmax=0.93*(k)/(par_num)+0.05
pos=[0.06,pydynmin,0.82,pydynmax]
;
; Calculate the histogram
;
if (total(pn(k-1, *)) NE 0) then yp=HISTOGRAM (pn(k-1,*), BINSIZE=bin)
if (total(m1(k-1, *)) NE 0) then ym1=HISTOGRAM (m1(k-1,*), BINSIZE=bin)
if (total(m2(k-1, *)) NE 0) then ym2=HISTOGRAM (m2(k-1,*), BINSIZE=bin)
if (total(r1(k-1, *)) NE 0) then yr1=HISTOGRAM (r1(k-1,*), BINSIZE=bin)
if (total(r2(k-1, *)) NE 0) then yr2=HISTOGRAM (r2(k-1,*), BINSIZE=bin)
if (total(aciss_heg_m1(k-1, *)) NE 0) then yaciss_heg_m1 =HISTOGRAM (aciss_heg_m1(k-1,*) , BINSIZE=bin) *1.0
if (total(aciss_heg_p1(k-1, *)) NE 0) then yaciss_heg_p1 =HISTOGRAM (aciss_heg_p1(k-1,*) , BINSIZE=bin) *1.0
if (total(aciss_leg_m1(k-1, *)) NE 0) then yaciss_leg_m1 =HISTOGRAM (aciss_leg_m1(k-1,*) , BINSIZE=bin) *1.0
if (total(aciss_leg_p1(k-1, *)) NE 0) then yaciss_leg_p1 =HISTOGRAM (aciss_leg_p1(k-1,*) , BINSIZE=bin) *1.0
if (total(aciss_meg_m1(k-1, *)) NE 0) then yaciss_meg_m1 =HISTOGRAM (aciss_meg_m1(k-1,*) , BINSIZE=bin) *1.0
if (total(aciss_meg_p1(k-1, *)) NE 0) then yaciss_meg_p1 =HISTOGRAM (aciss_meg_p1(k-1,*) , BINSIZE=bin) *1.0
if (total(hrci_leg_m123(k-1, *)) NE 0) then yhrci_leg_m123=HISTOGRAM (hrci_leg_m123(k-1,*), BINSIZE=bin) *1.0
if (total(hrci_leg_p123(k-1, *)) NE 0) then yhrci_leg_p123=HISTOGRAM (hrci_leg_p123(k-1,*), BINSIZE=bin) *1.0
if (total(hrcs_leg_m123(k-1, *)) NE 0) then yhrcs_leg_m123=HISTOGRAM (hrcs_leg_m123(k-1,*), BINSIZE=bin) *1.0
if (total(hrcs_leg_p123(k-1, *)) NE 0) then yhrcs_leg_p123=HISTOGRAM (hrcs_leg_p123(k-1,*), BINSIZE=bin) *1.0
;
; Calculate the X-axis for the histogram
;
if (total(pn(k-1, *)) NE 0) then xp=min(pn(k-1,*))+dindgen(n_elements(yp))*bin
if (total(m1(k-1, *)) NE 0) then xm1=min(m1(k-1,*))+dindgen(n_elements(ym1))*bin
if (total(m2(k-1, *)) NE 0) then xm2=min(m2(k-1,*))+dindgen(n_elements(ym2))*bin
if (total(r1(k-1, *)) NE 0) then xr1=min(r1(k-1,*))+dindgen(n_elements(yr1))*bin
if (total(r2(k-1, *)) NE 0) then xr2=min(r2(k-1,*))+dindgen(n_elements(yr2))*bin
if (total(aciss_heg_m1(k-1, *)) NE 0) then xyaciss_heg_m1 =min(aciss_heg_m1(k-1,*) )+dindgen(n_elements(yaciss_heg_m1 ))*bin
if (total(aciss_heg_p1(k-1, *)) NE 0) then xyaciss_heg_p1 =min(aciss_heg_p1(k-1,*) )+dindgen(n_elements(yaciss_heg_p1 ))*bin
if (total(aciss_leg_m1(k-1, *)) NE 0) then xyaciss_leg_m1 =min(aciss_leg_m1(k-1,*) )+dindgen(n_elements(yaciss_leg_m1 ))*bin
if (total(aciss_leg_p1(k-1, *)) NE 0) then xyaciss_leg_p1 =min(aciss_leg_p1(k-1,*) )+dindgen(n_elements(yaciss_leg_p1 ))*bin
if (total(aciss_meg_m1(k-1, *)) NE 0) then xyaciss_meg_m1 =min(aciss_meg_m1(k-1,*) )+dindgen(n_elements(yaciss_meg_m1 ))*bin
if (total(aciss_meg_p1(k-1, *)) NE 0) then xyaciss_meg_p1 =min(aciss_meg_p1(k-1,*) )+dindgen(n_elements(yaciss_meg_p1 ))*bin
if (total(hrci_leg_m123(k-1, *)) NE 0) then xyhrci_leg_m123=min(hrci_leg_m123(k-1,*))+dindgen(n_elements(yhrci_leg_m123))*bin
if (total(hrci_leg_p123(k-1, *)) NE 0) then xyhrci_leg_p123=min(hrci_leg_p123(k-1,*))+dindgen(n_elements(yhrci_leg_p123))*bin
if (total(hrcs_leg_m123(k-1, *)) NE 0) then xyhrcs_leg_m123=min(hrcs_leg_m123(k-1,*))+dindgen(n_elements(yhrcs_leg_m123))*bin
if (total(hrcs_leg_p123(k-1, *)) NE 0) then xyhrcs_leg_p123=min(hrcs_leg_p123(k-1,*))+dindgen(n_elements(yhrcs_leg_p123))*bin
;
; Plot the histogram
;
plot, xp, yp, position=pos, symsize=aymaize, thick=3, psym=10, $
xtitle="flux ratio",font=12, xstyle=4, ystyle=1,charsize=7./par_num,$
ytitle=name(k)+" keV", /noerase, /nodata, xrange=[0.7,1.3], yrange=[min(yp), max(yp)*1.1]
if (total(pn(k-1, *)) NE 0) then oplot, xp, yp, color=0, psym=10, thick=4
if (total(m1(k-1, *)) NE 0) then oplot, xm1, ym1, color=240, psym=10, thick=4
if (total(m2(k-1, *)) NE 0) then oplot, xm2, ym2, color=160, psym=10, thick=4
if (total(r1(k-1, *)) NE 0) then oplot, xr1, yr1, color=60, psym=10, thick=4
if (total(r2(k-1, *)) NE 0) then oplot, xr2, yr2, color=110, psym=10, thick=4
if (total(aciss_heg_m1(k-1, *)) NE 0) then oplot,xyaciss_heg_m1 ,yaciss_heg_m1 , psym=10, thick=6, linestyle=1,color=190
if (total(aciss_heg_p1(k-1, *)) NE 0) then oplot,xyaciss_heg_p1 ,yaciss_heg_p1 , psym=10, thick=6, linestyle=1,color=230
if (total(aciss_leg_m1(k-1, *)) NE 0) then oplot,xyaciss_leg_m1 ,yaciss_leg_m1 ,psym=10, thick=6, linestyle=1 ,color=30
if (total(aciss_leg_p1(k-1, *)) NE 0) then oplot,xyaciss_leg_p1 ,yaciss_leg_p1 ,psym=10, thick=6, linestyle=1 ,color=40
if (total(aciss_meg_m1(k-1, *)) NE 0) then oplot,xyaciss_meg_m1 ,yaciss_meg_m1 , psym=10, thick=6, linestyle=1,color=80
if (total(aciss_meg_p1(k-1, *)) NE 0) then oplot,xyaciss_meg_p1 ,yaciss_meg_p1 , psym=10, thick=6, linestyle=1,color=130
if (total(hrci_leg_m123(k-1, *)) NE 0) then oplot,xyhrci_leg_m123 ,yhrci_leg_m123 , psym=10, thick=6, linestyle=1,color=205
if (total(hrci_leg_p123(k-1, *)) NE 0) then oplot,xyhrci_leg_p123 ,yhrci_leg_p123 , psym=10, thick=6, linestyle=1,color=210
if (total(hrcs_leg_m123(k-1, *)) NE 0) then oplot,xyhrcs_leg_m123 ,yhrcs_leg_m123 , psym=10, thick=6, linestyle=1,color=215
if (total(hrcs_leg_p123(k-1, *)) NE 0) then oplot,xyhrcs_leg_p123 ,yhrcs_leg_p123 , psym=10, thick=6, linestyle=1,color=220
;
; Plot the axis according to the position of the panel in the plot
;
axis, xaxis=1, xstyle=1, xtickname=strarr(30)+' '
if (k EQ 1) then axis, xaxis=0, xstyle=1, font=12
if (k NE 1) then axis,xaxis=0, xstyle=1, font=12
;
; Plot labels outside the plot
;
outstring="statistical basis: "+string(e, format='(I4)')+" observations"
xyouts, 0.5,0.985, outstring, /normal, charsize=1., charthick=1, font=12, ALIGNMENT=0.5
xyouts, 0.5,1.0, " relative flux in different bands [Joint fit = 1]", /normal, charsize=1.5, charthick=1, font=12,ALIGNMENT=0.5
xyouts, 0.85,0.95, "pn ", /normal, charsize=1.3, charthick=1.5, font=12, color=0
xyouts, 0.85,0.925, "MOS1 ", /normal, charsize=1.3, charthick=1.5, font=12, color=240
xyouts, 0.85,0.90, "MOS2 ", /normal, charsize=1.3, charthick=1.5, font=12, color=160
xyouts, 0.85,0.875, "RGS1 ", /normal, charsize=1.3, charthick=1.5, font=12, color=60
xyouts, 0.85,0.85, "RGS2 ", /normal, charsize=1.3, charthick=1.5, font=12, color=110
;
; Select values of the flux different from 0 for the calculation of
; the statistical moments
;
filterp=where (xp NE 0)
xp=xp(filterp)
yp=yp(filterp)
filterm1=where (xm1 NE 0)
xm1=xm1(filterm1)
ym1=ym1(filterm1)
filterm2=where (xm2 NE 0)
xm2=xm2(filterm2)
ym2=ym2(filterm2)
filterr1=where (xr1 NE 0)
xr1=xr1(filterr1)
yr1=yr1(filterr1)
filterr2=where (xr2 NE 0)
xr2=xr2(filterr2)
yr2=yr2(filterr2)
;
; Calculate the statistical moments
;
if (total(pn(k-1, *)) NE 0) then meanp=total(xp*yp)/total(yp)
if (total(pn(k-1, *)) NE 0) then sdevp=sqrt(total((xp-meanp)^2*yp^2)/(total(yp)-1))
if (total(pn(k-1, *)) NE 0) then pn_momlab='!9m!X'+'='+string(format='(1(f5.3))',meanp)+' ('+'!9s!X'+'!D!9m!X!N'+'='+string(format='(1(f5.3))',sdevp)+')'
if (total(m1(k-1, *)) NE 0) then meanm1=total(xm1*ym1)/total(ym1)
if (total(m1(k-1, *)) NE 0) then sdevm1=sqrt(total((xm1-meanm1)^2*ym1^2)/(total(ym1)-1))
if (total(m1(k-1, *)) NE 0) then mos1_momlab='!9m!X'+'='+string(format='(1(f5.3))',meanm1)+' ('+'!9s!X'+'!D!9m!X!N'+'='+string(format='(1(f5.3))',sdevm1)+')'
if (total(m2(k-1, *)) NE 0) then meanm2=total(xm2*ym2)/total(ym2)
if (total(m2(k-1, *)) NE 0) then sdevm2=sqrt(total((xm2-meanm2)^2*ym2^2)/(total(ym2)-1))
if (total(m2(k-1, *)) NE 0) then mos2_momlab='!9m!X'+'='+string(format='(1(f5.3))',meanm2)+' ('+'!9s!X'+'!D!9m!X!N'+'='+string(format='(1(f5.3))',sdevm2)+')'
if (total(r1(k-1, *)) NE 0) then meanr1=total(xr1*yr1)/total(yr1)
if (total(r1(k-1, *)) NE 0) then sdevr1=sqrt(total((xr1-meanr1)^2*yr1^2)/(total(yr1)-1))
if (total(r1(k-1, *)) NE 0) then rgs1_momlab='!9m!X'+'='+string(format='(1(f5.3))',meanr1)+' ('+'!9s!X'+'!D!9m!X!N'+'='+string(format='(1(f5.3))',sdevr1)+')'
if (total(r2(k-1, *)) NE 0) then meanr2=total(xr2*yr2)/total(yr2)
if (total(r2(k-1, *)) NE 0) then sdevr2=sqrt(total((xr2-meanr2)^2*yr2^2)/(total(yr2)-1))
if (total(r2(k-1, *)) NE 0) then rgs2_momlab='!9m!X'+'='+string(format='(1(f5.3))',meanr2)+' ('+'!9s!X'+'!D!9m!X!N'+'='+string(format='(1(f5.3))',sdevr2)+')'
;
; Print the statistical moments as labels
;
if (total(pn(k-1, *)) NE 0) then xyouts, 0.625,0.15+0.1165*(k-1), pn_momlab, /normal, charsize=1., charthick=1.3, font=12, color=0
if (total(m1(k-1, *)) NE 0) then xyouts, 0.625,0.135+0.1165*(k-1), mos1_momlab, /normal, charsize=1., charthick=1.3, font=12, color=240
if (total(m2(k-1, *)) NE 0) then xyouts, 0.625,0.12+0.1165*(k-1), mos2_momlab, /normal, charsize=1., charthick=1.3, font=12, color=160
if (total(r1(k-1, *)) NE 0) then xyouts, 0.625,0.105+0.1165*(k-1), rgs1_momlab, /normal, charsize=1., charthick=1.3, font=12, color=60
if (total(r2(k-1, *)) NE 0) then xyouts, 0.625,0.09+0.1165*(k-1), rgs2_momlab, /normal, charsize=1., charthick=1.3, font=12, color=110
;
; Plot specific Chandra labels
;
if (total(aciss_heg_m1(k-1, *)) NE 0) then xyouts, 0.85,0.70, "...aciss_heg_m1", /normal, charsize=1.3, charthick=1.5 ,color=190
if (total(aciss_heg_p1(k-1, *)) NE 0) then xyouts, 0.85,0.65, "...aciss_heg_p1", /normal, charsize=1.3, charthick=1.5 ,color=230
if (total(aciss_leg_m1(k-1, *)) NE 0) then xyouts, 0.85,0.60, "...aciss_leg_m1", /normal, charsize=1.3, charthick=1.5 ,color=30
if (total(aciss_leg_p1(k-1, *)) NE 0) then xyouts, 0.85,0.55, "...aciss_leg_p1", /normal, charsize=1.3, charthick=1.5 ,color=40
if (total(aciss_meg_m1(k-1, *)) NE 0) then xyouts, 0.85,0.50, "...aciss_meg_m1", /normal, charsize=1.3, charthick=1.5 ,color=80
if (total(aciss_meg_p1(k-1, *)) NE 0) then xyouts, 0.85,0.45, "...aciss_meg_p1", /normal, charsize=1.3, charthick=1.5 ,color=130
if (total(hrci_leg_m123(k-1, *)) NE 0) then xyouts, 0.85,0.40, "...hrci_leg_m123", /normal, charsize=1.3, charthick=1.5,color=205
if (total(hrci_leg_p123(k-1, *)) NE 0) then xyouts, 0.85,0.35, "...hrci_leg_p123", /normal, charsize=1.3, charthick=1.5,color=210
if (total(hrcs_leg_m123(k-1, *)) NE 0) then xyouts, 0.85,0.30, "...hrcs_leg_m123", /normal, charsize=1.3, charthick=1.5,color=215
if (total(hrcs_leg_p123(k-1, *)) NE 0) then xyouts, 0.85,0.25, "...hrcs_leg_p123", /normal, charsize=1.3, charthick=1.5,color=220
;
; Plot footer labels
;
xyouts, 0.5,0.045, "Flux ratio normalised to the flux of the joint fit", /normal, charsize=1., charthick=1., font=12, ALIGNMENT=0.5
xyouts, 0.5,0.025, versions(vi), /normal, charsize=1., charthick=1., font=12, ALIGNMENT=0.5
endfor
device,/close
set_plot,'x'
endfor
;
; That's all Folks!
;
end