; JaeSub Hong, 2003-2005, version 1.7
; Please report any problem or suggestion at jaesub@head.cfa.harvard.edu
; 
; Calculate quantiles from a distribution
; Refer to 
;	J. Hong, E.M. Schlegel & J.E. Grindlay, 
;	2004 ApJ 614 p508 and references therein
; 
; quantile diagram plotting routine
; now with fits format for grid files
; and the default x axis is log of median
; require rdrdb.pro (included)

;========================================================================

pro qd_setwin, xr=xr, yr=yr, $
	xtit=xtit,ytit=ytit, $
	xty=xty,yty=yty, $
	xtick=xt, range=range, xtv=xtv, $
	noxtick=noxtick, noytick=noytick, $
	noxtoptick=noxtoptick, $
	stu=stu, noerase=noerase, func=func
; setting up window

if not keyword_set(func) then func=['log10g','ratio']

Elo=str(range[0],format='(f10.1)')
Ehi=str(range[1],format='(f10.1)')
if not keyword_set(xtit) then begin
	if func[0] eq 'logit' then begin
		xtit='log!d10!n(!s!a m !r!s!b1-m!r!s-!r !s-!r !s-!r -!3!n)' 
	endif else if func[0] eq 'sqrt' then begin
		xtit='log!d10!n(m)' 
	endif else if func[0] eq 'log10' then begin
		xtit='log!d10!n(10 m + 1-m)' 
	endif else if func[0] eq 'log10b' then begin
		xtit='log!d10!n(100 m + 1-m)/2' 
	endif else if func[0] eq 'log10g' then begin
		xtit='log!d10!n(E!d50!n/'+Elo+')/log!d10!n('+Ehi+'/'+Elo+')' 
	endif
endif
if not keyword_set(ytit) then begin
	ytit='3 (E!d25!n-'+Elo+')/(E!d75!n-'+Elo+')'
endif

if not keyword_set(xt) or keyword_set(noxtoptick) then xs=1 else xs=9

xtickn=''
ytickn=''
if keyword_set(noxtick) then begin
	xtickn=strarr(10)+ ' '
	xtit=' '
endif
if keyword_set(noytick) then begin
	ytickn=strarr(10)+ ' '
	ytit=' '
endif
if not keyword_set(noerase) then noerase=0

plot,[0],[0],/nodata, $
 	xr=xr,xs=xs,yr=yr,/ys,$
	xty=xty,yty=yty, $
	xtickn=xtickn,ytickn=ytickn, $
	xtit=xtit,ytit=ytit, noerase=noerase

; labeled xtick
if keyword_set(xt) and not keyword_set(noxtoptick) then begin
;xt= ['0.6','1','2','3','4','5','6 keV','7'];

	rl = range[1]-range[0]
	if not keyword_set(xtv) then m = (float(xt)-range[0])/rl $
	else m = (xtv-range[0])/rl 

	xv =alog10(m/(1.-m))

	if func[0] eq 'logit' then begin
		xv = alog10(m/(1.-m))
	endif else if func[0] eq 'sqrt' then begin
		xv = sqrt(m)
	endif else if func[0] eq 'log10' then begin
		xv = alog10(10.*m+(1.-m))
	endif else if func[0] eq 'log10b' then begin
		xv = alog10(100.*m+(1.-m))/2.
	endif else if func[0] eq 'log10g' then begin
		rr=range[1]/range[0]
		xv = alog10(rr*m+(1.-m))/alog10(rr)
	endif

	if keyword_set(noxtoptick) then xt=strarr(10)+' '
	axis, xaxis=1,xticks=n_elements(xt)-1, xtickv=xv, xtickn=xt;,xticklen=0.01
; unlabeled small xtick

if keyword_set(stu) then begin
	min_m = (10.^xr[0]/(10.^xr[0]+1.))*rl+range[0]
	max_m = (10.^xr[1]/(10.^xr[1]+1.))*rl+range[0]

	min_m = fix(min_m/stu/10.-1)*stu*10
	max_m = fix(max_m/stu/10.+1)*stu*10
	n_sxt = (max_m-min_m)/stu
	if n_sxt gt 60. then begin
		print,'# of ticks should be less than 60'
		goto, done
	endif
	sxtv = findgen(n_sxt)*stu+min_m

	w=where(sxtv gt range[0] and sxtv lt range[1],n_sxt)
	if n_sxt eq 0 then goto, done
	sxtv=sxtv[w]

	sxtv = (sxtv-range[0])/rl
	if func[0] eq 'logit' then begin
		sxtv = alog10(sxtv/(1.-sxtv))
	endif else if func[0] eq 'sqrt' then begin
		sxtv = sqrt(sxtv)
	endif else if func[0] eq 'log10' then begin
		sxtv = alog10(10.*sxtv+(1.-sxtv))
	endif else if func[0] eq 'log10b' then begin
		sxtv = alog10(100.*sxtv+(1.-sxtv))/2.
	endif else if func[0] eq 'log10g' then begin
		rr=range[1]/range[0]
		sxtv = alog10(rr*sxtv+(1.-sxtv))/alog10(rr)
	endif
	sxt = strarr(n_sxt)+' '
	axis, xaxis=1,xticks=n_sxt-1, xtickv=sxtv, xtickn=sxt, xticklen=0.01
	done:
endif
endif


end

function qd_readgrid, file_grid, comment=comment, n_data=n_data, format=format

; read grid
; if it's multiple, take the average
; assuming the identical format

if not keyword_set(format) then format='fits'

n_file	= n_elements(file_grid)
n_file_=0
n_data=0
for i=0, n_file-1 do begin
;	print,i,file_grid[i]
	if format eq 'fits' then begin
		grid_ = mrdfits(file_grid[i],1,comment,/silent)
		n_data_=n_elements(grid_)
	endif else if format eq 'rdb' then begin
		grid_ = rdrdb2(file_grid[i],n_data=n_data_,comment=comment)
	endif
	if n_data_ eq 0 then begin
		print,'warning no data in ',file_grid[i]
		goto, skip
	endif
	if n_file_ eq 0 then begin	
		ntag 	= n_tags(grid_)
		grid	= grid_
	endif else begin
		for j=3, ntag-1 do grid.(j) = grid.(j)+grid_.(j)
	endelse
	n_file_ = n_file_+1
	if n_data lt n_data_ then n_data=n_data_
	skip:
endfor

if n_file_ gt 0 then begin
	for j=3, ntag-1 do grid.(j) = grid.(j)/n_file_
;	help,grid
	return, grid
endif else return, 0

end

pro qd_grid2qd, grid, comment=comment, $ 
	range=range, cid=cid, $
	qdx=qdx, qdy=qdy, func=func, format=format

if not keyword_set(format) then format='fits'

if not keyword_set(range) then begin
	if format eq 'fits' then begin
		found='no'
		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+9,20))
				strreplace,str_range,"'",""
				found='yes'
				break
			endif
		endfor

	endif else if format eq 'rdb' then begin
		found='no'
		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))
				found='yes'
				break
			endif
		endfor
	endif
	if found eq 'yes' then begin
		fields = strsplit(str_range,'[:,]',/extract,/reg)
		range = [float(fields[0]),float(fields[1])]
	endif 
endif

rl = range[1]-range[0]

; this is safer, but keep compatiblity 
if keyword_set(cid) then begin
	Q25 = (grid.(cid[0])-range[0])/rl
	Q50 = (grid.(cid[1])-range[0])/rl
	Q75 = (grid.(cid[2])-range[0])/rl
endif else begin
	Q25 = (grid.E25-range[0])/rl
	Q50 = (grid.E50-range[0])/rl
	Q75 = (grid.E75-range[0])/rl
endelse

if not keyword_set(func) then func=['logit','ratio']
if func[0] eq 'logit' then begin
	QDx = alog10(Q50/(1.-Q50))
endif else if func[0] eq 'sqrt' then begin
	QDx = sqrt(Q50)
endif else if func[0] eq 'log10' then begin
	QDx = alog10(10.*Q50+(1.-Q50))
endif else if func[0] eq 'log10b' then begin
	QDx = alog10(100.*Q50+(1.-Q50))/2.
endif else if func[0] eq 'log10g' then begin
	rr=range[1]/range[0]
	QDx = alog10(rr*Q50+(1.-Q50))/alog10(rr)
endif

if func[1] eq 'ratio' then begin
	QDy = 3.*Q25/Q75
endif else if func[1] eq 'sqrt_ratio' then begin
	QDy = sqrt(Q25/Q75)
endif
;print,rr, range
;print,min(q50),max(q50)
;print,min(qdx),max(qdx)
;print,min(qdy),max(qdy)
end

pro qd_plotgrid, grid, comment=comment, $ 
	file_grid=file_grid, $ ; read the data
	range=range, $
	line=line, thick=thick, $
	stressx=stressx, stressy=stressy, $
	col=colgrid, cid=cid, full=full, func=func

if not keyword_set(line) 	then line=0
if not keyword_set(colgrid) 	then colgrid=0
if not keyword_set(thick) 	then thick=!p.thick

if n_elements(thick) 	eq 1 	then thick=[thick,thick]
if n_elements(line) 	eq 1 	then line=[line,line]
if n_elements(colgrid) 	eq 1 	then colgrid=[colgrid,colgrid]

if keyword_set(file_grid) then $
	grid=qd_readgrid(file_grid, comment=comment)

qd_grid2qd, grid, comment=comment, $ 
	range=range, cid=cid, $
	qdx=QDx, qdy=QDy, func=func

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

; now the first parameter
w=where(gidt eq 'gx',cw)
gidn1 = gidn[w]
QDx_ = QDx[w]
QDy_ = QDy[w]

if not keyword_set(stressx) 	then stressx=fltarr(max(gidn1)+1)
if keyword_set(full) 	then begin
	for i=0, max(gidn1) do begin
		w=where(gidn1 eq i,cw)
		oplot, QDx_[w], QDy_[w], line=line[0], $
			col=colgrid[0],thick=thick[0],psym=full
	endfor
	return
endif

for i=0, max(gidn1) do begin
	w=where(gidn1 eq i,cw)
	if cw gt 0 then $
	oplot, QDx_[w], QDy_[w], line=line[0], $
		col=colgrid[0],thick=thick[0]+stressx[i];,psym=7
endfor

; next the second parameter
w=where(gidt eq 'gy',cw)
gidn2 = gidn[w]
QDx_ = QDx[w]
QDy_ = QDy[w]

if not keyword_set(stressy) 	then stressy=fltarr(max(gidn2)+1)
for i=0, max(gidn2) do begin
	w=where(gidn2 eq i,cw)
	if cw gt 0 then $
	oplot, QDx_[w], QDy_[w], line=line[1], $
		col=colgrid[1],thick=thick[1]+stressy[i];,psym=7
endfor

end

pro qd_plot_multigrid, file_grid, grid=grid, comment=comment, $
	type=type, col=col,$
	range=range, $
	line=line, thick=thick, cid=cid, $
	stressx=stressx, stressy=stressy, full=full, func=func

if not keyword_set(type) then type='all'
if type eq 'all' then begin
	for i=0, n_elements(file_grid)-1 do begin
		fs = file_search(file_grid[i])
		;print,file_grid[i]
		if fs[0] ne '' then begin
			grid_=qd_readgrid(file_grid[i], comment=comment,n_data=n_data)
			if n_data gt 0 then begin
				qd_plotgrid, grid_, comment=comment, col=col, $
					range=range, line=line, thick=thick, cid=cid, $
					stressx=stressx, stressy=stressy, full=full, func=func
				grid=grid_
			endif
		endif else print,"can't find ", file_grid[i]
	endfor
endif else if type eq 'avg' then begin
	grid=qd_readgrid(file_grid, comment=comment)
	qd_plotgrid, grid, comment=comment, col=col, $
			range=range, line=line, thick=thick, cid=cid, $
			stressx=stressx, stressy=stressy, full=full, func=func
endif

end

pro qd_labelgrid, par, grid, $
	label1=label1, label2=label2, $
	prefix1=prefix1, prefix2=prefix2, $
	suffix1=suffix1, suffix2=suffix2, $
	comment=comment, $ ; when grid data is given directly
	file_grid=file_grid, $ ; read the data
	format=format, $
	range=range, $
	align=align, ori=ori, col=col, xo=xo, yo=yo, cid=cid, func=func

if not keyword_set(format) then format='fits'

if keyword_set(file_grid) then begin
	if format eq 'fits' then begin
		grid = mrdfits(file_grid,1,comment,/silent)
		n_data=n_elements(grid)
	endif else if format eq 'rdb' then begin
		grid = rdrdb2(file_grid,n_data=n_data,comment=comment)
	endif

	if n_data eq 0 then return
endif

if not keyword_set(align) then align=[0.0,0.0]
if not keyword_set(ori) then ori=[0.0,0.0]
if not keyword_set(xo) then xo=[0.0,0.0]
if not keyword_set(yo) then yo=[0.0,0.0]

ncomm=n_elements(comment)
if format eq 'fits' then begin
	bingo=[0,0,0]
	for i=0, ncomm-1 do begin
		pos=strpos(comment[i], 'RANGE',0)
		if pos ge 0 then begin
			str_range = strtrim(strmid(comment[i],pos+9,20))
			strreplace,str_range,"'",""
			bingo[0]=1
		endif
		pos=strpos(comment[i], 'PAR1',0)
		if pos ge 0 then begin
			str_par1 = strtrim(strmid(comment[i],pos+9,80))
			fields = strsplit(str_par1, '/',/extract)
			str_par1=fields[0]
			strreplace,str_par1," ",""
			strreplace,str_par1,"'",""
			bingo[1]=1
		endif
		pos=strpos(comment[i], 'PAR2',0)
		if pos ge 0 then begin
			str_par2 = strtrim(strmid(comment[i],pos+9,80))
			fields = strsplit(str_par2, '/',/extract)
			str_par2=fields[0]
			strreplace,str_par2," ",""
			strreplace,str_par2,"'",""
			bingo[2]=1
		endif
		if total(bingo) eq 3 then break
	endfor
endif else if format eq 'rdb' then begin
	bingo=[0,0,0]
	for i=0, ncomm-1 do begin
		pos=strpos(comment[i], '-range',0)
		if pos ge 0 then begin
			str_range = strtrim(strmid(comment[i],pos+7,20))
			bingo[0]=1
		endif
		pos=strpos(comment[i], '-par1',0)
		if pos ge 0 then begin
			str_par1 = strtrim(strmid(comment[i],pos+6,80))
			bingo[1]=1
		endif
		pos=strpos(comment[i], '-par2',0)
		if pos ge 0 then begin
			str_par2 = strtrim(strmid(comment[i],pos+6,80))
			bingo[2]=1
		endif
		if total(bingo) eq 3 then break
	endfor
endif

if not keyword_set(range) then begin
	fields = strsplit(str_range,'[:,]',/extract,/reg)
	range = [float(fields[0]),float(fields[1])]
endif

if not keyword_set(prefix1) then prefix1=''
if not keyword_set(prefix2) then prefix2=''
if not keyword_set(suffix1) then suffix1=''
if not keyword_set(suffix2) then suffix2=''

fields = strsplit(str_par1,':',/extract)
par1v = fields[1]
par_ = strsplit(par1v,',',/extract)
npar_ = n_elements(par_)
par1v = [prefix1,strarr(npar_)]+par_+[suffix1,strarr(npar_)]

fields = strsplit(str_par2,':',/extract)
par2v = fields[1]
par_ = strsplit(par2v,',',/extract)
npar_ = n_elements(par_)
par2v = [prefix2,strarr(npar_)]+par_+[suffix2,strarr(npar_)]


rl = range[1]-range[0]

; this is safer, but keep compatiblity 
if keyword_set(cid) then begin
	Q25 = (grid.(cid[0])-range[0])/rl
	Q50 = (grid.(cid[1])-range[0])/rl
	Q75 = (grid.(cid[2])-range[0])/rl
endif else begin
	Q25 = (grid.E25-range[0])/rl
	Q50 = (grid.E50-range[0])/rl
	Q75 = (grid.E75-range[0])/rl
endelse

if not keyword_set(func) then func=['logit','ratio']
if func[0] eq 'logit' then begin
	QDx = alog10(Q50/(1.-Q50))
endif else if func[0] eq 'sqrt' then begin
	QDx = sqrt(Q50)
endif else if func[0] eq 'log10' then begin
	QDx = alog10(10.*Q50+(1.-Q50))
endif else if func[0] eq 'log10b' then begin
	QDx = alog10(100.*Q50+(1.-Q50))/2.
endif else if func[0] eq 'log10g' then begin
	rr=range[1]/range[0]
	QDx = alog10(rr*Q50+(1.-Q50))/alog10(rr)
endif

if func[1] eq 'ratio' then begin
	QDy = 3.*Q25/Q75
endif else if func[1] eq 'sqrt_ratio' then begin
	QDy = sqrt(Q25/Q75)
endif
gidt = strmid(grid.gridid,0,2)
gidn = fix(strmid(grid.gridid,2,3))

w=where(gidt eq 'gx',cw)
gidn1 = gidn[w]
grid1_ = grid[w].(1)
grid2_ = grid[w].(2)
QDx_ = QDx[w]
QDy_ = QDy[w]

for i=0, max(gidn1) do begin
	w=where(gidn1 eq i,cw)
	dum=min(abs(grid2_[w] - par[1]),ind)
	if n_elements(label1)-1 ge i then label_ = label1[i] $
	else label_ = par1v[i]
	xyouts, QDx_[w[ind]]+xo[0], QDy_[w[ind]]+yo[0], label_, $
		align=align[0], ori=ori[0],noclip=0,col=col[0]
endfor

w=where(gidt eq 'gy',cw)
gidn2 = gidn[w]
grid1_ = grid[w].(1)
grid2_ = grid[w].(2)
QDx_ = QDx[w]
QDy_ = QDy[w]

for i=0, max(gidn2) do begin
	w=where(gidn2 eq i,cw)
	dum=min(abs(grid1_[w] - par[0]),ind)
	if n_elements(label2)-1 ge i then label_ = label2[i] $
	else label_ = par2v[i]
	xyouts, QDx_[w[ind]]+xo[1], QDy_[w[ind]]+yo[1], label_, $
		align=align[1], ori=ori[1],noclip=0,col=col[1]
endfor


end

pro qd_labelgrid_set, type, grid, comment=comment, $
	col=col, $
	file_grid=file_grid, range=range, cid=cid, func=func, forma=format
; predefined grid label
if type eq '20' then begin
	align =[1.0,0.5]  ; label alignment
	ori   =[-20,0]	  ; orientation
	par   =[0.01,5]	  ; location
	xo    =[0.0,0.0]  ; x offset
	yo    =[0.0,-0.08]; y offset
	prefix1= ['','','','','','','','','','','N!dH!n(x10!u22!n)=','']
	prefix2= ['','','','','','','','!4C!3=']
	suffix1 = ''
	suffix2 = ''
endif else if type eq '21' then begin
	align =[1.0,0.]  ; label alignment
	ori   =[-20,45]	  ; orientation
	par   =[100,0.1]	  ; location
	xo    =[0.0,0.02]  ; x offset
	yo    =[0.0,0.02]; y offset
	prefix1= ['','','','N!dH!n(x10!u22!n)=','','','','','','','']
	prefix2 = ''
	suffix1= ['','','','','',' ','','','','','']
	suffix2 = ['','=kT(keV)','','','','','','','','','']
endif else if type eq 'pl' then begin
	align =[1.0,0.5]  ; label alignment
	ori   =[-20,0]	  ; orientation
	par   =[0.01,4]	  ; location
	xo    =[0.0,0.0]  ; x offset
	yo    =[0.0,-0.05]; y offset
	prefix1= ['','','','','','N!dH22!n=','','','','','']
	prefix2= '!4C!3='
	suffix1 = ''
	suffix2 = ''
endif else if type eq 'pl2' then begin
	align =[1.0,0.5]  ; label alignment
	ori   =[-20,0]	  ; orientation
	par   =[0.01,4]	  ; location
	xo    =[0.0,0.0]  ; x offset
	yo    =[0.0,-0.12]; y offset
	prefix1= ['','','','','','N!dH22!n=','','','','','']
	prefix2= ['','','','','!4C!3=','','']
	suffix1 = ''
	suffix2 = ''
endif else if type eq 'pl3' then begin
	align =[1.0,0.5]  ; label alignment
	ori   =[-20,0]	  ; orientation
	par   =[0.01,4]	  ; location
	xo    =[0.0,0.0]  ; x offset
	yo    =[0.03,-0.10]; y offset
	prefix1= ['','','','','','N!dH22!n=','','','','','']
	prefix2= ['','','','!4C!3=','','']
	suffix1 = ''
	suffix2 = ''
endif else if type eq 'br2' then begin
	align =[1.0,1.0]  ; label alignment
	ori   =[-10,010.]	  ; orientation
	par   =[0.01,1.0]	  ; location
	xo    =[0.02,-0.02]  ; x offset
	yo    =[0.02,-0.02]; y offset
	prefix1= ['','','','','','N!dH22!n=','','','','','']
	prefix2= ['','','','kT(keV)=','','']
	suffix1 = ''
	suffix2 = ''
endif else if type eq 'br3' then begin
	align =[1.0,1.0]  ; label alignment
	ori   =[-10,010.]	  ; orientation
	par   =[0.01,1.0]	  ; location
	xo    =[-0.00,-0.005]  ; x offset
	yo    =[0.02,-0.035]; y offset
	prefix1= ['','','','','','N!dH22!n=','','','','','']
	prefix2= ['','','','kT(keV)=','','']
	suffix1 = ''
	suffix2 = ''
endif else if type eq 'br4' then begin
	align =[1.0,1.0]  ; label alignment
	ori   =[-10,010.]	  ; orientation
	par   =[0.01,3.0]	  ; location
	xo    =[-0.00,-0.005]  ; x offset
	yo    =[0.02,-0.035]; y offset
	prefix1= ['','','','','','N!dH22!n=','','','','','']
	prefix2= ['','','','kT(keV)=','','']
	suffix1 = ''
	suffix2 = ''
endif else if type eq 'pll' then begin
	align =[1.0,0.5]  ; label alignment
	ori   =[-20,0]	  ; orientation
	par   =[0.01,5]	  ; location
	xo    =[0.0,0.0]  ; x offset
	yo    =[0.0,-0.05]; y offset
	prefix1= ['','','','','','','','N!dH!n(x10!u22!n)=','','','','','']
	prefix2= '!4C!3='
	suffix1 = ''
	suffix2 = ''
endif else if type eq 'br5' then begin
	align =[1.0,0.]  ; label alignment
	ori   =[-20,45]	  ; orientation
	par   =[10,0.2]	  ; location
	xo    =[0.0,0.02]  ; x offset
	yo    =[0.0,0.02]; y offset
;	prefix1= ['','','','','N!dH!n(x10!u22!n)=','','','','','','']
	prefix1= ['','','','','N!dH22!n=','','','','','','']
	prefix2 = ''
	suffix1= ['','','','','',' ','','','','','']
	suffix2 = ['','','=kT(keV)','','','','','','','','','']
endif else if type eq 'br6' then begin
	align =[1.0,0.]  ; label alignment
	ori   =[-20,45]	  ; orientation
	par   =[10,0.2]	  ; location
	xo    =[0.0,0.005]  ; x offset
	yo    =[0.0,0.005]; y offset
;	prefix1= ['','','','','N!dH!n(x10!u22!n)=','','','','','','']
	prefix1= ['','','','','N!dH22!n=','','','','','','']
	prefix2 = ''
	suffix1= ['','','','','',' ','','','','','']
	suffix2 = ['','','=kT(keV)','','','','','','','','','']
endif else if type eq 'br' or type eq 'bb' then begin 
	align =[1.0,0.]  ; label alignment
	ori   =[-20,45]	  ; orientation
	par   =[10,0.2]	  ; location
	xo    =[0.0,0.02]  ; x offset
	yo    =[0.0,0.02]; y offset
;	prefix1= ['','','','','N!dH!n(x10!u22!n)=','','','','','','']
	prefix1= ['','','','','N!dH22!n=','','','','','','']
	prefix2 = ''
	suffix1= ['','','','','',' ','','','','','']
	suffix2 = ['','=kT(keV)','','','','','','','','','']
endif else if type eq 'brl' or type eq 'bbl' then begin 
	align =[1.0,0.]  ; label alignment
	ori   =[-20,45]	  ; orientation
	par   =[10,0.2]	  ; location
	xo    =[0.0,0.02]  ; x offset
	yo    =[0.0,0.02]; y offset
	prefix1= ['','','','','N!dH!n(x10!u22!n)=','','','','','','']
	prefix2 = ''
	suffix1= ['','','','','',' ','','','','','']
	suffix2 = ['','=kT(keV)','','','','','','','','','']
endif else if type eq 'fermi_cpl10a_3' then begin
	align =[0.0,0.5]  ; label alignment
	ori   =[20,-30]	  ; orientation
	par   =[1.e6,1.9]	  ; location
	xo    =[0.01,-0.01]  ; x offset
	yo    =[0.01,0.01]; y offset
	prefix1= ['','','','','','','','','','','','','']
	prefix2= [' ',' ','','','','','','','','','','','']
	suffix1 = '  '
	suffix2= ['','','',' ','','','','','','','','','']
	suffix1= ['',' GeV',' GeV',' GeV','','','','','','','','','']
endif else if type eq 'fermi_cpl10a_2' then begin
	align =[0.0,0.5]  ; label alignment
	ori   =[60,-20]	  ; orientation
	par   =[2.e5,0.5]	  ; location
	xo    =[0.01,-0.01]  ; x offset
	yo    =[0.00,0.01]; y offset
	prefix1= ['','','','','','','','','','','','','']
	prefix2= [' ',' ','','','','','','','','','','','']
	suffix1 = '  '
	suffix2= ['','','',' ','','','','','','','','','']
	suffix1= ['','',' GeV',' GeV','','','','','','','','','']
endif else if type eq 'fermi_cpl10a_1' then begin
	align =[1.2,0.5]  ; label alignment
	ori   =[00,00]	  ; orientation
	par   =[5.e5,1.3]	  ; location
	xo    =[0.01,-0.01]  ; x offset
	yo    =[0.00,0.01]; y offset
	prefix1= ['','','','','','','','','','','','','']
	prefix2= [' ',' ','','','','','','','','','','','']
	suffix1 = '  '
	suffix2= ['','','','','','','','','','','','','']
	suffix1= [' GeV','','','','','','','','','','','','']
endif else if type eq 'fermi_cpl10a' then begin
	align =[.5,0.5]  ; label alignment
	ori   =[00,00]	  ; orientation
	par   =[5.e6,0.5]	  ; location
	xo    =[0.00,0.0]  ; x offset
	yo    =[0.03,-0.08]; y offset
	prefix1= ['E!dc!n=','','','','','','','','','','','','']
	prefix2= [' ',' ','','','','','!4C!3=','','','','','','']
	suffix1 = '  '
	suffix2= ['','','','','','','     ','','','','','','']
	suffix1= [' GeV   ','','','','','','','','','','','','']
endif else if type eq 'fermi_cpl' then begin
	align =[.5,0.5]  ; label alignment
	ori   =[00,00]	  ; orientation
	par   =[1.e7,0.5]	  ; location
	par   =[1.e9,0.5]	  ; location
	xo    =[0.00,0.0]  ; x offset
	yo    =[0.03,-0.08]; y offset
	prefix1= ['E!dc!n=','','','','','','','','','','','','']
	prefix2= [' ',' ','','','','','!4C!3=','','','','','','']
	suffix1 = '  '
	suffix2= ['','','','','','','     ','','','','','','']
	suffix1= [' GeV   ','','','','','','','','','','','','']
endif else if type eq 'fermi_cpl2' then begin
	align =[.5,0.5]  ; label alignment
	ori   =[00,00]	  ; orientation
	par   =[1.e7,0.5]	  ; location
	par   =[1.e9,0.5]	  ; location
	xo    =[0.00,0.0]  ; x offset
	yo    =[0.01,-0.04]; y offset
	prefix1= ['E!dc!n=','','','','','','','','','','','','']
	prefix2= [' ',' ','','','!4C!3=','','','','','','']
	suffix1 = '  '
	suffix2= ['','','','','   ','     ','','','','','','']
	suffix1= [' GeV   ','','','','','','','','','','','','']
endif else if type eq 'plc' then begin
	align =[.0,0.5]  ; label alignment
	ori   =[20,0]	  ; orientation
	par   =[0.01,0]	  ; location
	xo    =[0.03,0.0]  ; x offset
	yo    =[-0.05,-0.06]; y offset
	prefix1= ['','','','','','','','','','','','','']
	prefix2= ''
	suffix1 = ''
	suffix2 = ''
endif else if type eq 'brc' then begin 
	align =[1.0,0.5]  ; label alignment
	ori   =[-20,0]	  ; orientation
	par   =[0.01,0.1]	  ; location
	xo    =[0.0,0.02]  ; x offset
	yo    =[0.0,-0.05]; y offset
	prefix1= ['','','','','','','','','','','']
	prefix2 = ''
	suffix1= ['','','','','',' ','','','','','']
	suffix2 = ['','','','','','','','','','','']
endif else if type eq '4' then begin
	align =[0,0.]  ; label alignment
	ori   =[-45,45]	  ; orientation
	par   =[10,100]	  ; location
	xo    =[0.00,0.02]  ; x offset
	yo    =[-0.06,0.02]; y offset
	prefix1= ['','','','','','','','','']
	prefix2 = ''
	suffix1= ['=N!dH!n(x10!u22!n)','','',' ','','','','','']
	suffix2 = ['','=kT(keV)','','','','','','','','','']
endif else if type eq 'tb2' then begin
	align =[1.0,0.]  ; label alignment
	ori   =[-20,45]	  ; orientation
	par   =[1000,0.2]	  ; location
	xo    =[0.0,0.02]  ; x offset
	yo    =[0.0,0.02]; y offset
	prefix1= ['','','','','N!dH!n(x10!u22!n)=','','','','','','']
	prefix2 = ''
	suffix1= ['','','','','',' ','','','','','']
	suffix2 = ['','','=kT(keV)','','','','','','','','','']
endif else begin
	align =[1.0,0.]  ; label alignment
	ori   =[-20,45]	  ; orientation
	par   =[10,0.2]	  ; location
	xo    =[0.0,0.02]  ; x offset
	yo    =[0.0,0.02]; y offset
	prefix1= ['','','','','N!dH!n(x10!u22!n)=','','','','','','']
	prefix2 = ''
	suffix1= ['','','','','',' ','','','','','']
	suffix2 = ['','=kT(keV)','','','','','','','','','']
endelse

;print,'-============ col',col
qd_labelgrid, par, grid, $
	prefix1=prefix1, prefix2=prefix2, $
	suffix1=suffix1, suffix2=suffix2, $
	comment=comment, $ ; when grid data is given directly
	file_grid=file_grid, $ ; read the data
	range=range, $
	align=align, ori=ori, col=col, xo=xo, yo=yo, cid=cid, func=func, format=format

end


pro qd_plot2, range=range, xr=xr, yr=yr, $
	xtick=xtick, col=col, stu=stu, $
	stressx=stressx, stressy=stressy, $
	skip_grid=skip_grid, $
	label=label, multi=grid_multi, $
	line=line, noxtick=noxtick, noytick=noytick, $
	noxtoptick=noxtoptick, xtit=xtit, ytit=ytit, $
	setwin=setwin, xtv=xtv, func=func, $
	file_grid=file_grid, type=grid_type, format=format, grid=grid
	

if not keyword_set(range) 	then range=[0.3,8.0]
if not keyword_set(xr) 		then xr    =[-1.5,0.5]
if not keyword_set(yr) 		then yr    =[0.5,2.5]
if not keyword_set(xtick) 	then xtick =['0.5','1','2','3','4','5','6','7'] ; top x-axis
if not keyword_set(col) 	then col   =[60,250]	  ; color setting
if 0 eq  n_elements(stu) 	then stu   =0.5
if not keyword_set(stressx) 	then stressx   =intarr(20)
if not keyword_set(stressy) 	then stressy   =intarr(20)
if not keyword_set(skip_grid) 	then skip_grid ='no'
if 0 eq n_elements(label) 	then label =1
if not keyword_set(grid_multi) 	then grid_multi ='avg'
if not keyword_set(line) 	then line=0
if not keyword_set(noxtick) 	then noxtick=0
if not keyword_set(noytick) 	then noytick=0
if not keyword_set(noxtoptick) 	then noxtoptick=0


Elo=str(range[0],format='(f10.1)')
Ehi=str(range[1],format='(f10.1)')

if not keyword_set(xtit)	then xtit='log!d10!n(E!d50!n/'+Elo+')/log!d10!n('+Ehi+'/'+Elo+')' 
if not keyword_set(ytit)	then ytit='3 (E!d25!n-'+Elo+')/(E!d75!n-'+Elo+')'

if not keyword_set(setwin) then $
	qd_setwin, xr=xr, yr=yr, range=range, xtick=xtick, stu=stu, $
		xtit=xtit, ytit=ytit, noxtick=noxtick, noytick=noytick, $
		noxtoptick=noxtoptick, func=func, xtv=xtv

if skip_grid eq 'no' then begin
	qd_plot_multigrid, file_grid, grid=grid, comment=comment, $
		type=grid_multi, col=col, thick=!p.thick, cid=cid, $
		stressx=stressx,stressy=stressy, line=line, range=range, func=func
;	help,grid
	if label ne 0 then $
		qd_labelgrid_set, grid_type, grid, comment=comment, $
			col=col, cid=cid, range=range, func=func
endif

setwin=0


end


pro qdplot2
;	print, "qccd plot2: version 0.1"
end