;------------------------------------------------------------------------------------------
;Subroutine:	wband_momentsplot
;Purpose:	Routine to determine make a 3-panel radar moments plot
;I/O format:	wband_momentsplot,time,height,dbz,dopp,spwd,stn,datestr
;Inputs:	time:	time axis
;		height: height axis
;		dbz:	reflectivity [dBZ]
;		dopp:	Doppler velocity [m/s]
;		spwd:	spectral width [m/s]
;		stn:	signal-to-noise ratio [dBZ]
;		datestr: date string for the plot title
;Keywords	XRA:	Range in time for plotting
;		YRA:	Range in height for plotting
;		SITE: 	Radar site for the title bar.
;		YLOG:	Y axis is logarithmic if set
;		CBASE:  cloud base heights to include in plot if desired
;Output:	PNG plot
;Author:	Matthew Shupe
;Modified:   7/7/2014
;------------------------------------------------------------------------------------------
pro wband_momentsplot,time,height,dbz,dopp,spwd,stn,datestr,XRA=xra,YRA=yra,SITE=site,YLOG=ylog,CBASE=cbase

;colorfile='/home/dataman/IDLcode/shupetools/colorscopy.tbl'
colorfile='c:/data/ACSE/programs/colorscopy.tbl'

;----------------------------------------------------------
;Create the plot window
;Set up all of the margins and individual plot positions.
; (Note:  all "marg" values are in normalized units)
;----------------------------------------------------------
window,0,xsize=650,ysize=650,title='WBAND Visualization'
topmarg=0.042
bottommarg=0.07
lftmarg=0.09
rgtmarg=0.14
wid=(1-(topmarg+bottommarg))/3
pos1=[lftmarg,bottommarg+2*wid,1-rgtmarg,1-topmarg]
pos2=[lftmarg,bottommarg+wid,1-rgtmarg,bottommarg+2*wid]
pos3=[lftmarg,bottommarg,1-rgtmarg,bottommarg+wid]
texttemp='!5blah'
if not(keyword_set(SITE)) then sitest=' ' else sitest=site+', '
;-----------------------------------------------------------
;Load the colortable for plotting radar data
;-----------------------------------------------------------
loadct,44,file=colorfile
;------------------------------------------------------------
;Filter out data where S-to-N is lower than a threshold
;------------------------------------------------------------
pdbz=dbz & pdopp=dopp & pspwd=spwd 
threshold = -15.
iwh=where(stn le threshold)
if iwh[0] ne -1 then begin
	pdbz[iwh]=-99
	pdopp[iwh]=-9
	pspwd[iwh]=-0.1
endif
pdbz[*,n_elements(height)-1]=-999

;-----------------------------------------------------------
;Set up the plot axis ranges.
;-----------------------------------------------------------
if not(keyword_set(yra)) then yra=[0,13]
if not(keyword_set(xra)) then xra=[time[0],time[n_elements(time)-1]]
if xra[1] le xra[0] then xra[1]=xra[1]+24.
if not(keyword_set(ylog)) then ylog=0
;-----------------------------------------------------------
;Plot reflectivity
;	Set up the contour "levels" and "colors"
;-----------------------------------------------------------
!p.color=1
!p.background=0
levels=[-60.,findgen(29)*2.5-50.]
numc=!d.n_colors
if numc gt 256 then numc=256
colcol=fix(indgen(28)*(numc*.78/27))+fix(numc*.196)
colors=[0,colcol,numc-1]
ncblev=29
nlev=n_elements(levels)
contour,pdbz,time,height,levels=levels,c_colors=colors,/cell_fill,/xstyle,/ystyle,yrange=yra,xrange=xra,ylog=ylog,$
	position=pos1,/norm,font=-1,xtickname=strarr(20)+'  ',title=strupcase(sitest)+'WBAND data for '+datestr
if keyword_set(cbase) then oplot,time,cbase,color=1,psym=2,symsize=0.1
;----------------------------------------------------------
;Create the reflectivity colorbar and labels
;----------------------------------------------------------
cben=pos1[3]-.01
cbst=pos1[1]+.05
cbywid=(cbst-cben)/(ncblev-1)
cb_x=[pos1[2]+.02,pos1[2]+.02,pos1[2]+.05,pos1[2]+.05]
cb_y=[cbst-cbywid,cbst,cbst,cbst-cbywid]
cbl_x=pos1[2]+.055
cbstr=strtrim(string(levels,format='(f5.1)'),2)
for i=(nlev-ncblev),nlev-2 do begin
	polyfill,cb_x,cb_y,color=colors[i],/norm
	if i mod 4 eq 1 then xyouts,cbl_x,cb_y[1]+.5*cbywid,cbstr[i],font=-1,/norm
	cb_y=cb_y-cbywid
	if i eq (nlev-2) then xyouts,cbl_x,cb_y[1]+.5*cbywid,cbstr[i+1],font=-1,/norm
endfor
xyouts,pos1[2]+0.04,pos1[1]+.025,/norm,'dBZ',font=-1
xyouts,pos1[0]+0.02,pos1[3]-0.04,/norm,'!6Reflectivity!5',font=-1,charsize=1.8,color=1,charthick=2
;----------------------------------------------------------
;Plot Doppler Velocity
;	Set up the contour "levels" and "colors"
;----------------------------------------------------------
levels=[-10,findgen(25)/(24/5.4)-2.7]
nlev=n_elements(levels)
colors=[0, fix(indgen(nlev-2)*(numc*.80/(nlev-3)))+fix(numc*.196), numc-1]
ncblev=nlev-1
contour, pdopp,time,height,/noerase,levels=levels,c_colors=colors,/cell_fill,position=pos2,/norm,font=-1,$
	/xstyle,/ystyle,yrange=yra,xrange=xra,xtickname=strarr(20)+' ',ytitle='Height AGL (km)',ylog=ylog
if keyword_set(cbase) then oplot,time,cbase,color=1,psym=2,symsize=0.1
;----------------------------------------------------------
;Create the velocity colorbar and label
;----------------------------------------------------------
cben=pos2[3]-.01
cbst=pos2[1]+.05
cbywid=(cbst-cben)/(ncblev-1)
cb_x=[pos2[2]+.02,pos2[2]+.02,pos2[2]+.05,pos2[2]+.05]
cb_y=[cbst-cbywid,cbst,cbst,cbst-cbywid]
cbl_x=pos2[2]+.055
cbstr=strtrim(string(levels,format='(f4.1)'),2)
for i=(nlev-ncblev),nlev-2 do begin
	polyfill,cb_x,cb_y,color=colors[i],/norm
	if (i mod 4)-1 eq 0 then xyouts,cbl_x,cb_y[1]+.5*cbywid,cbstr[i],font=-1,/norm
	cb_y=cb_y-cbywid
	if i eq (nlev-2) then xyouts,cbl_x,cb_y[1]+.5*cbywid,cbstr[i+1],font=-1,/norm
endfor
xyouts,cbl_x+.03,cb_y[1]+.5*cbywid,'down',font=-1,/norm
xyouts,cbl_x+.045,cb_y[1]+(i-.5)*cbywid,'up',font=-1,/norm
xyouts,pos2[2]+0.04,pos2[1]+.025,/norm,'m/s',font=-1
xyouts,pos2[0]+0.02,pos2[3]-0.04,/norm,'!6Doppler Velocity!5',font=-1,charsize=1.8,color=1,charthick=2
;---------------------------------------------------------
;Plot Spectral Width
;	Setup the contour "levels" and "colors"
;----------------------------------------------------------
levels=[-0.01,findgen(31)/20.]
nlev=n_elements(levels)
ncblev=nlev-1
colors=[0,fix(indgen(nlev-2)*(.80*numc/(nlev-3)))+fix(.196*numc), numc-2,numc-1]
contour, pspwd,time,height,/noerase,levels=levels,c_colors=colors,/cell_fill,/xstyle,/ystyle,yrange=yra,xrange=xra,$
	xtickformat='timeticks',position=pos3,/norm,font=-1,xtitle='Time (hours, GMT)',ylog=ylog
if keyword_set(cbase) then oplot,time,cbase,color=1,psym=2,symsize=0.1
;---------------------------------------------------------
;Create the spectral width colorbar and labels
;----------------------------------------------------------
cben=pos3[3]-.01
cbst=pos3[1]+.05
cbywid=(cbst-cben)/(ncblev-1)
cb_x=[pos3[2]+.02,pos3[2]+.02,pos3[2]+.05,pos3[2]+.05]
cb_y=[cbst-cbywid,cbst,cbst,cbst-cbywid]
cbl_x=pos3[2]+.055
cbstr=strtrim(string(levels,format='(f4.2)'),2)
for i=(nlev-ncblev),nlev-2 do begin
	polyfill,cb_x,cb_y,color=colors[i],/norm
	if i mod 5 eq 1 then xyouts,cbl_x,cb_y[1]+.5*cbywid,cbstr[i],font=-1,/norm
	cb_y=cb_y-cbywid
	if i eq (nlev-2) then xyouts,cbl_x,cb_y[1]+.5*cbywid,cbstr[i+1],font=-1,/norm
endfor
xyouts,pos3[2]+0.04,pos3[1]+.025,/norm,'m/s',font=-1
xyouts,pos3[0]+0.02,pos3[3]-0.04,/norm,'!6Spectral Width!5',font=-1,charsize=1.8,color=1,charthick=2
end ;prog
;------------------------------------------------------------------------------------------