; plot the example grid pattern
;
; To run this program, do
;
;	IDL> .r plot_grid.pro
;
; JaeSub Hong, 2003-2005, version 1.5
; Please report any problem or suggestion at jaesub@head.cfa.harvard.edu
; 

; read the file
; choose your file
file_grid = "test_powerlaw.rdb"
file_grid = "pl_835.rdb"
file_grid = "br_835_old.rdb"
file_grid = "br_835.rdb"
file_grid = "gaussian.rdb"
file_grid = "edge.rdb"
file_grid = "lorentz.rdb"
file_grid = "redge.rdb"
file_grid = "rs_835.rdb"
file_grid = "rs_835_.rdb"
file_grid = "rs_835_2.rdb"
file_grid = "rs_2.rdb"
file_grid = "rs_.rdb"
file_grid = "rs.rdb"
file_grid = "test_powerlaw.rdb"
file_grid = "test_bremss.rdb"
n_data = rdrdb(file_grid,data=grid,comment=comment)

; get the energy range and par values from the comment
; this is for cosmetics, not really necessary.
for i=0, n_elements(comment)-1 do begin
	pos=strpos(comment[i], '-range',0)
	if pos ge 0 then begin
		str_range = strtrim(strmid(comment[i],pos+7,20))
		goto,next
	endif
	pos=strpos(comment[i], '-par1',0)
	if pos ge 0 then begin
		str_par1 = strtrim(strmid(comment[i],pos+6,40))
		goto,next
	endif
	pos=strpos(comment[i], '-par2',0)
	if pos ge 0 then begin
		str_par2 = strtrim(strmid(comment[i],pos+6,40))
		goto,next
	endif
	next:
endfor

fields = strsplit(str_range,'[:,]',/extract,/reg)
Elo = float(fields[0])
Ehi = float(fields[1])
Er = Ehi-Elo

fields = strsplit(str_par1,':',/extract)
par1v = fields[1]
par1v = strsplit(par1v,',',/extract)
par1n = fields[0]
par1n = strsplit(par1n,'.',/extract)
par1n = par1n[1]

fields = strsplit(str_par2,':',/extract)
par2v = fields[1]
par2n = fields[0]
par2v = strsplit(par2v,',',/extract)
par2n = strsplit(par2n,'.',/extract)
par2n = par2n[1]

; convert Ex% to Q50 and Q25/Q75
Q25 = (grid.E0-Elo)/Er
Q50 = (grid.E1-Elo)/Er
Q75 = (grid.E2-Elo)/Er

; set the axes
x=alog10(Q50/(1-Q50))
y=3.*Q25/Q75
xr=[min(x),max(x)]
yr=[min(y),max(y)]
xr=xr+(xr[1]-xr[0])*[-0.2,0.2]
yr=yr+(yr[1]-yr[0])*[-0.2,0.2]

; set the window
plot,[0],[0], $
	/nodata, $
	xr=xr,/xs, $
	yr=yr,/ys, $
	xtit='log!d10!n(!s!a m !r!s!b1-m!r!s-!r !s-!r !s-!r -!3!n)', $
	ytit='3 Q!d25!n/Q!d75!n'

; plot grid patterns
gidt = strmid(grid.gridid,0,2)
gidn = fix(strmid(grid.gridid,2,3))

w=where(gidt eq 'gx',cw)
x1 = x[w]
y1 = y[w]
gidn1 = gidn[w]

w=where(gidt eq 'gy',cw)
x2 = x[w]
y2 = y[w]
gidn2 = gidn[w]


xn=-0.05
yn=0.03
plotsym,0,0.5,/fill
for i=0, max(gidn1) do begin
	w=where(gidn1 eq i,cw)
	oplot,x1[w],y1[w],psym=8
	if i eq max(gidn1) then pref=par1n+'=' $
	else pref=''
	xyouts,x1[w[0]]+xn,y1[w[0]]+yn,pref+par1v[i]
endfor

xn=-0.03
yn=-0.05
for i=0, max(gidn2) do begin
	w=where(gidn2 eq i,cw)
	oplot,x2[w],y2[w],col=250,psym=8
	if i eq max(gidn2) then pref='='+par2n $
	else pref=''
	xyouts,x2[w[0]]+xn,y2[w[0]]+yn,par2v[i]+pref,col=250
endfor



end