pro xspecplot,x,y1,y2,y1err=y1err,y2err=y2err,xerr=xerr,$
	y1over=y1over,y2over=y2over,$
	fraclow=fraclow,verbose=verbose,outfile=outfile,$
	xtitle=xtitle,title=title,$
	Ulstyl=Ulstyl,Upsym=Upsym,Ucol=Ucol,Uylog=Uylog,Uytitle=Uytitle,Uyr=Uyr,$
	Llstyl=Llstyl,Lpsym=Lpsym,Lcol=Lcol,Lylog=Lylog,Lytitle=Lytitle,Lyr=Lyr,$
	oUlstyl=oUlstyl,oUpsym=oUpsym,oUcol=oUcol,oLlstyl=oLlstyl,oLpsym=oLpsym,oLcol=oLcol,$
	linestyle=linestyle,psym=psym,color=color,ylog=ylog,ytitle=ytitle,yrange=yrange,$
	_extra=e
;+
;procedure	xspecplot
;	plot two curves in two panels the way XSPEC does spectra and
;	residuals
;
;syntax
;	xspecplot,x,y1,y2,y1err=y1err,y2err=y2err,xerr=xerr,$
;	y1over=y1over,y2over=y2over,$
;	fraclow=fraclow,verbose=verbose,outfile=outfile,$
;	Ulstyl=Ulstyl,Upsym=Upsym,Ucol=Ucol,Uylog=Uylog,Uytitle=Uytitle,Uyr=Uyr,$
;	Llstyl=Llstyl,Lpsym=Lpsym,Lcol=Lcol,Lylog=Lylog,Lytitle=Lytitle,Lyr=Lyr,$
;	oUlstyl=oUlstyl,oUpsym=oUpsym,oUcol=oUcol,oLlstyl=oLlstyl,oLpsym=oLpsym,oLcol=oLcol,$
;	[PLOT keywords]
;	
;
;parameters
;	x	[INPUT; required] abscissa array
;	y1	[INPUT; required] array that goes in the top panel
;	y2	[INPUT; required] array that goes in the bottom panel
;		* sizes of all arrays _must_ match
;
;keywords
;	y1err	[INPUT] error bar on Y1
;	y2err	[INPUT] error bar on Y2
;	xerr	[INPUT] error bar on Y2
;		* if scalar and +ve, assumed to be a fractional error
;		  if between 0 and 1, and %error if greater than 1
;		* if scalar and -ve, abs value is assumed to be a constant
;		  error
;		* if array of same size as X, assumed to be symmetric
;		* if array of size [2,NX] or [NX,2], assumed to be upper
;		  and lower bounds
;		* ignored if array size doesn't match X
;	y1over	[INPUT] additional arrays to overplot in upper panel
;	y2over	[INPUT] additional arrays to overplot in lower panel
;		* size must be [m,NX] or [NX,m]
;		* if not set, both assumed to be 0*NX
;	fraclow	[INPUT] what fraction of the plot should the lower panel
;		take up?
;		* default is 0.35
;		* hard minimum of 0.05 and hard maximum of 0.95
;	verbose	[INPUT] controls chatter
;	outfile	[INPUT] if set, makes a postscript plot
;		* if set to a string, assumes that that is the name of
;		  the file, otherwise writes to ./xspecplot.ps
;	Ulstyl	[INPUT] line style for upper panel
;	Llstyl	[INPUT] line style for lower panel
;	linestyle	[INPUT] caught and discarded
;	Upsym	[INPUT] PSYM for upper panel
;	Lpsym	[INPUT] PSYM for lower panel
;	psym	[INPUT] caught and discarded
;	Ucol	[INPUT] COLOR for upper panel
;	Lcol	[INPUT] COLOR for lower panel
;	color	[INPUT] caught and discarded
;	Uylog	[INPUT] YLOG for upper panel
;	Lylog	[INPUT] YLOG for lower panel
;	ylog	[INPUT] caught and discarded
;	Uytitle	[INPUT] YTITLE for upper panel
;	Lytitle	[INPUT] YTITLE for lower panel
;	ytitle	[INPUT] caught and discarded
;	Uyr	[INPUT] YRANGE for upper panel
;	Lyr	[INPUT] YRANGE for lower panel
;	oUlstyl	[INPUT] line styles of Y1OVER (must be array of size M)
;	oLlstyl	[INPUT] line styles of Y2OVER (must be array of size M)
;	oUpsym	[INPUT] PSYM for Y1OVER (must be array of size M)
;	oLpsym	[INPUT] PSYM for Y2OVER (must be array of size M)
;	oUcol	[INPUT] COLOR for Y1OVER (must be array of size M)
;	oLcol	[INPUT] COLOR for Y2OVER (must be array of size M)
;	yrange	[INPUT] caught and discarded
;	
;	_extra	[INPUT ONLY] pass defined keywords to PLOT
;
;history
;	vinay kashyap (2012-jun)
;-

;	usage
ok='ok' & np=n_params() & nx=n_elements(x) & n1=n_elements(y1) & n2=n_elements(y2)
if np lt 3 then ok='Insufficient parameters' else $
 if nx eq 0 then ok='X is not defined' else $
  if n1 eq 0 then ok='Y1 is not defined' else $
   if n2 eq 0 then ok='Y2 is not defined' else $
    if nx ne n1 then ok='X and Y1 do not match' else $
     if nx ne n2 then ok='X and Y2 do not match'
if ok ne 'ok' then begin
  print,'Usage: xspecplot,x,y1,y2,y1err=y1err,y2err=y2err,xerr=xerr,$'
  print,'       y1over=y1over,y2over=y2over,$'
  print,'       fraclow=fraclow,verbose=verbose,outfile=outfile,$'
  print,'       Ulstyl=Ulstyl,Upsym=Upsym,Ucol=Ucol,Uylog=Uylog,Uytitle=Uytitle,Uyr=Uyr,$'
  print,'       Llstyl=Llstyl,Lpsym=Lpsym,Lcol=Lcol,Lylog=Lylog,Lytitle=Lytitle,Lyr=Lyr,$'
  print,'       oUlstyl=oUlstyl,oUpsym=oUpsym,oUcol=oUcol,oLlstyl=oLlstyl,oLpsym=oLpsym,oLcol=oLcol,$'
  print,'       [PLOT keywords]'
  print,'  make XSPEC-like 2-panel plots'
  if np ne 0 then message,ok,/informational
  return
endif

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

if keyword_set(outfile) then begin
  dname=!d.NAME
  if size(outfile,/type) eq 7 then psfile=strtrim(outfile[0],2) else $
  	psfile='xspecplot.ps'
  if vv gt 1 then message,'writing to '+psfile,/informational
  set_plot,'ps' & device,file=psfile,/inches,ysize=10,yoffset=1,xsize=8,xoffset=1,/color
endif
;
n1e=n_elements(y1err) & n2e=n_elements(y2err) & nxe=n_elements(xerr)
y1u=0.*y1 & y1l=y1u & y2u=y1u & y2l=y1u & xu=y1u & xl=y1u
;
if n1e eq 1 then begin	;(scalar
  y1e=y1err[0]
  if y1e lt 0 then begin
    y1u[*]=abs(y1e) & y1l[*]=abs(y1e)
  endif
  if y1e gt 0 and y1e lt 1 then begin
    y1u=abs(y1e*y1) & y1l=abs(y1e*y1)
  endif
endif else begin	;scalar)(array
  if n1e eq nx then begin	;(size matches NX
    y1u=y1err & y1l=y1err
  endif else begin		;N=NX)(does not match NX
    if n1e eq 2*nx then begin	;(N=2*NX
      sz=size(y1err) & m1=sz[1] & m2=sz[2]
      if m1 eq 2 then begin & y1u=reform(y1err[0,*]) & y1l=reform(y1err[1,*]) & endif
      if m2 eq 2 then begin & y1u=reform(y1err[*,0]) & y1l=reform(y1err[*,1]) & endif
    endif			;N=NX*2)
  endelse			;N.ne.NX)
endelse			;array)
;
if n2e eq 1 then begin	;(scalar
  y2e=y2err[0]
  if y2e lt 0 then begin
    y2u[*]=abs(y2e) & y2l[*]=abs(y2e)
  endif
  if y2e gt 0 and y2e lt 1 then begin
    y2u=abs(y2e*y2) & y2l=abs(y2e*y2)
  endif
endif else begin	;scalar)(array
  if n2e eq nx then begin	;(size matches NX
    y2u=y2err & y2l=y2err
  endif else begin		;N=NX)(does not match NX
    if n2e eq 2*nx then begin	;(N=2*NX
      sz=size(y2err) & m1=sz[1] & m2=sz[2]
      if m1 eq 2 then begin & y2u=reform(y2err[0,*]) & y2l=reform(y2err[1,*]) & endif
      if m2 eq 2 then begin & y2u=reform(y2err[*,0]) & y2l=reform(y2err[*,1]) & endif
    endif			;N=NX*2)
  endelse			;N.ne.NX)
endelse			;array)
;
if nxe eq 1 then begin	;(scalar
  xe=xerr[0]
  if xe lt 0 then begin
    xu[*]=abs(xe) & xl[*]=abs(xe)
  endif
  if xe gt 0 and xe lt 1 then begin
    xu=abs(xe*x) & xl=abs(xe*x)
  endif
endif else begin	;scalar)(array
  if nxe eq nx then begin	;(size matches NX
    xu=xerr & xl=xerr
  endif else begin		;N=NX)(does not match NX
    if nxe eq 2*nx then begin	;(N=2*NX
      sz=size(xerr) & m1=sz[1] & m2=sz[2]
      if m1 eq 2 then begin & xu=reform(xerr[0,*]) & xl=reform(xerr[1,*]) & endif
      if m2 eq 2 then begin & xu=reform(xerr[*,0]) & xl=reform(xerr[*,1]) & endif
    endif			;N=NX*2)
  endelse			;N.ne.NX)
endelse			;array)
;
flo=0.35 & if keyword_set(fraclow) then flo=(abs(fraclow[0])>0.05)<0.95
;
if not keyword_set(Upsym) then Upsym=10 & if not keyword_set(Lpsym) then Lpsym=4
if not keyword_set(Ucol) then Ucol=255-!p.background & if not keyword_set(Lcol) then Lcol=255-!p.background
if not keyword_set(Uylog) then Uylog=0 & if not keyword_set(Lylog) then Lylog=0
if not keyword_set(Uytitle) then Uytitle='Y1' & if not keyword_set(Lytitle) then Lytitle='Y2'
if not keyword_set(Uyr) then Uyr=minmax(y1) & if not keyword_set(Lyr) then Lyr=minmax(y2)
if n_elements(Uyr) ne 2 then begin
  if Uyr[0] gt min(y1) then Uyr=[min(y1),Uyr[0]] else Uyr=[Uyr[0],max(y1)]
endif
if n_elements(Lyr) ne 2 then begin
  if Lyr[0] gt min(y2) then Lyr=[min(y2),Lyr[0]] else Lyr=[Lyr[0],max(y2)]
endif
;
if not keyword_set(xtitle) then xtitle='X'
if not keyword_set(title) then title=''
;
y1o=0*y1 & y2o=0*y2 & n1o=n_elements(y1over) & n2o=n_elements(y2over)
if n1o gt 0 then begin
  if n1o eq 1 then y1o[*]=y1over[0] else begin
    if n1o eq nx then y1o=y1over else begin
      sz=size(y1over)
      if sz[0] eq 2 then begin
        if sz[1] eq nx or sz[2] eq nx then begin
	  if sz[1] eq nx then y1o=transpose(y1over) else y1o=y1over
	endif else begin
          if vv gt 1 then message,'cannot understand format of Y1OVER; ignoring',/informational
	endelse
      endif else begin
        if vv gt 1 then message,'cannot understand format of Y1OVER; ignoring',/informational
      endelse
    endelse
  endelse
endif
if n2o gt 0 then begin
  if n2o eq 1 then y2o[*]=y2over[0] else begin
    if n2o eq nx then y2o=y2over else begin
      sz=size(y2over)
      if sz[0] eq 2 then begin
        if sz[1] eq nx or sz[2] eq nx then begin
	  if sz[1] eq nx then y2o=transpose(y2over) else y2o=y2over
	endif else begin
          if vv gt 1 then message,'cannot understand format of Y2OVER; ignoring',/informational
	endelse
      endif else begin
        if vv gt 1 then message,'cannot understand format of Y2OVER; ignoring',/informational
      endelse
    endelse
  endelse
endif
sz=size(y1o) & if sz[0] eq 2 then mn1=sz[1] else mn1=1
sz=size(y2o) & if sz[0] eq 2 then mn2=sz[1] else mn2=1
;
zoUlstyl=intarr(mn1) & zoLlstyl=intarr(mn2)
if n_elements(oUlstyl) eq 1 then zoUlstyl[*]=oUlstyl[0] else $
 if n_elements(oUlstyl) eq mn1 then zoUlstyl=oUlstyl
if n_elements(oLlstyl) eq 1 then zoLlstyl[*]=oLlstyl[0] else $
 if n_elements(oLlstyl) eq mn1 then zoLlstyl=oLlstyl
;
zoUpsym=intarr(mn1)-3 & zoLpsym=intarr(mn2)-3
if n_elements(oUpsym) eq 1 then zoUpsym[*]=oUpsym[0] else $
 if n_elements(oUpsym) eq mn1 then zoUpsym=oUpsym
if n_elements(oLpsym) eq 1 then zoLpsym[*]=oLpsym[0] else $
 if n_elements(oLpsym) eq mn1 then zoLpsym=oLpsym
;
zoUcol=intarr(mn1)+(255-!P.BACKGROUND) & zoLcol=intarr(mn2)+(255-!P.BACKGROUND)
if n_elements(oUcol) eq 1 then zoUcol[*]=oUcol[0] else $
 if n_elements(oUcol) eq mn1 then zoUcol=oUcol
if n_elements(oLcol) eq 1 then zoLcol[*]=oLcol[0] else $
 if n_elements(oLcol) eq mn1 then zoLcol=oLcol

;	initialize
pmulti=!p.multi & pposition=!p.position

!p.multi=[0,1,2]

;!p.position=[0.12,flo*0.84,0.97,0.94]
plot,/nodata,x,y1,ylog=Uylog,ytitle=Uytitle,yrange=Uyr,title=title,xticks=1,xtickname=replicate(' ',2),$
	position=[0.12,flo*0.84,0.97,0.94], _extra=e
oplot,x,y1,psym=Upsym,color=Ucol
for i=0,nx-1 do if y1u[i]+y1l[i] gt 0 then oplot,x[i]*[1,1],y1[i]+[-y1l[i],y1u[i]],color=Ucol
;
if mn1 eq 1 then oplot,x,y1o,linestyle=zoUlstyl[0],psym=zoUpsym[0],color=zoUcol[0], _extra=e else $
  for i=0,mn1-1 do oplot,x,y1o[i,*],linestyle=zoUlstyl[i],psym=zoUpsym[i],color=zoUcol[i], _extra=e

;!p.position=[0.12,0.1,0.97,flo*0.84]
plot,/nodata,x,y2,ylog=Lylog,ytitle=Lytitle,yrange=Lyr,xtitle=xtitle,$
	position=[0.12,0.1,0.97,flo*0.84], _extra=e
for i=0,nx-1 do if y2u[i]+y2l[i] gt 0 then oplot,x[i]*[1,1],y2[i]+[-y2l[i],y2u[i]],color=Lcol
;
oplot,x,y2,psym=Lpsym,color=Lcol
if mn2 eq 1 then oplot,x,y2o,linestyle=zoLlstyl[0],psym=zoLpsym[0],color=zoLcol[0], _extra=e else $
  for i=0,mn2-1 do oplot,x,y2o[i,*],linestyle=zoLlstyl[i],psym=zoLpsym[i],color=zoLcol[i], _extra=e

if vv gt 1000 then stop,'halting; type .CON to continue'

if keyword_set(outfile) then begin
  device,/close & set_plot,dname
  if vv gt 0 then spawn,'ls -l '+psfile
endif

!p.multi=pmulti & !p.position=pposition

return
end

peasecolr & loadct,3 & peasecolr

nx=361 & x=findgen(nx)*!pi/180. & y0=sin(x) & ysig=randomn(seed,nx)*0.1
y1=y0+ysig & y2=ysig
xspecplot,x,y1,y2,y1err=0.1,y2err=y2err,xerr=xerr,$
	y1over=y0,y2over=0,$
	fraclow=0.35,verbose=101,outfile=1,$
	xtitle=xtitle,title='X vs Y1,Y2',$
	Ulstyl=Ulstyl,Upsym=Upsym,Ucol=Ucol,Uylog=Uylog,Uytitle=Uytitle,Uyr=Uyr,$
	Llstyl=Llstyl,Lpsym=Lpsym,Lcol=Lcol,Lylog=Lylog,Lytitle=Lytitle,Lyr=Lyr,$
	oUlstyl=oUlstyl,oUpsym=oUpsym,oUcol=1,oLlstyl=1,oLpsym=oLpsym,oLcol=2,$
	/xs,/ys,xthick=2,ythick=2,thick=2,charthick=2,charsize=1.2

end