;HODOGRAPH.PRO
;AUTHOR:	JONATHAN GARNER
;CREATED:	SPRING 2005
;PURPOSE:	BRING IN A WIND PROFILE, PLOT A HODOGRAPH, AND
;			COMPUTE/DISPLAY VERTICAL WIND SHEAR PARAMETERS.

;FUNCTION TO PLOT CIRCLES ON THE HODOGRAPH
FUNCTION CIRCLE, xcenter, ycenter, radius
points = (2*!PI/99.0)*FINDGEN(100)
x=xcenter+radius*COS(points)
y=ycenter+radius*SIN(points)
RETURN, TRANSPOSE([[x],[y]])
END

;BRING IN WIND DATA
;FIRST PUT THE WIND DATA IN AN ASCII FILE-->'wind.dat'
;MAY HAVE TO MESS WITH THE TEXT LABELS 'level' IN ORDER TO
;FIT THE INCOMING DATA.

OPENR, lun, 'wind.dat', /GET_LUN
header=STRARR(3)
READF, lun, header

junk=""
rows=0
READS, header(1), junk, rows, FORMAT='(A17, x, I0)'

x=FLTARR(10000)
y=FLTARR(10000)
x0=0.0
y0=0.0
count=0

WHILE (NOT EOF(lun)) DO BEGIN
READF, lun, x0, y0
x(count)=x0
y(count)=y0
count=count+1
ENDWHILE

x=x(0:count-1)
y=y(0:count-1)

;FIND STORM MOTION

b=0
FOR i=0, 12 DO BEGIN
avg=(x[i] + b)
b=avg
avgx=avg/(i+1)
ENDFOR

b=0
FOR i=0, 12 DO BEGIN
avg=(y[i] + b)
b=avg
avgy=avg/(i+1)
ENDFOR

ushr=((x[12]+x[11])/2-(x[1]+x[0])/2)
vshr=((y[12]+y[11])/2-(y[1]+y[0])/2)
shear=sqrt(ushr^2+vshr^2)
umotion=[(avgx+(7.5/(shear))*vshr)]
vmotion=[(avgy-(7.5/(shear))*ushr)]
motion=[(umotion),(vmotion)]
meanwind=[(avgx),(avgy)]

;PLOT HODOGRAPH AND STORM MOTION

window, 0, xsize=900, ysize=750

level=['SFC','0.5 km','1 km','1.5 km','2 km','2.5 km','3 km','3.5 km','4 km','4.5 km','5 km', $
'5.5 km','6 km','6.5 km','7 km','7.5 km','8 km','8.5 km','9 km','9.5 km','10 km']
plot, x, y, xrange=[-20,40], yrange=[-19,39], ticklen=0, xstyle=4, ystyle=4, /isotropic
xyouts, x, y, level, color=200, charthick=1.1, align=1.2
plots, x, y, psym=6, symsize=.5, color=200
plots, motion, psym=2, symsize=1.5, color=150
srmotion='Vrm'
xyouts, umotion, vmotion, srmotion, color=200, charthick=1.1, align=-.5

;PLOT X AND Y AXES AND ALSO

xaxis=[[0,-30],[0,50]]
yaxis=[[-30,0],[60,0]]
plots, xaxis
plots, yaxis

plots, 10, 0, psym=1, symsize=1, color=100
plots, 20, 0, psym=1, symsize=1, color=100
plots, 30, 0, psym=1, symsize=1, color=100
plots, 40, 0, psym=1, symsize=1, color=100
plots, -10, 0, psym=1, symsize=1, color=100
plots, -20, 0, psym=1, symsize=1, color=100
plots, -30, 0, psym=1, symsize=1, color=100
plots, 0, 10, psym=1, symsize=1, color=100
plots, 0, 20, psym=1, symsize=1, color=100
plots, 0, 30, psym=1, symsize=1, color=100
plots, 0, 40, psym=1, symsize=1, color=100
plots, 0, -10, psym=1, symsize=1, color=100
plots, 0, -20, psym=1, symsize=1, color=100
plots, 0, -30, psym=1, symsize=1, color=100

;plots, CIRCLE(0, 0, 5), p=0, color=100
;plots, CIRCLE(0, 0, 10), p=0, color=100
;plots, CIRCLE(0, 0, 15), p=0, color=100
;plots, CIRCLE(0, 0, 20), p=0, color=100
;plots, CIRCLE(0, 0, 25), p=0, color=100
;plots, CIRCLE(0, 0, 30), p=0, color=100
;plots, CIRCLE(0, 0, 35), p=0, color=100
;plots, CIRCLE(0, 0, 40), p=0, color=100
;plots, CIRCLE(0, 0, 45), p=0, color=100
;plots, CIRCLE(0, 0, 50), p=0, color=100

radx10='10 m/s'
xyouts, 9.8, -1.2, radx10, color=100, charthick=1.1, align=1
radx20='20 m/s'
xyouts, 19.8, -1.2, radx20, color=100, charthick=1.1, align=1
radx30='30 m/s'
xyouts, 29.8, -1.2, radx30, color=100, charthick=1.1, align=1
radx40='40 m/s'
xyouts, 39.8, -1.2, radx40, color=100, charthick=1.1, align=1
rady10='10 m/s'
xyouts, 4, 10.5, rady10, color=100, charthick=1.1, align=1
rady20='20 m/s'
xyouts, 4, 20.5, rady20, color=100, charthick=1.1, align=1
rady30='30 m/s'
xyouts, 4, 30.5, rady30, color=100, charthick=1.1, align=1
rady40='40 m/s'
xyouts, 4, 40.5, rady40, color=100, charthick=1.1, align=1

;SHEAR HEADERS

shrtxt='BULK SHEAR:'
xyouts, -13.5, 40, shrtxt, color=200, charthick=1.25, align=2
srhtxt='STORM-RELATIVE HELICITY (SRH):'
xyouts, 3.95, 34, srhtxt, color=200, charthick=1.25, align=2
srtxt='STORM-RELATIVE WINDS:'
xyouts, -3.77, 28, srtxt, color=200, charthick=1.25, align=2
svrtxt='STORM EVALUATION:'
xyouts, -7.87, 20, svrtxt, color=200, charthick=1.25, align=2

;COMPUTE AND DISPLAY 0-6 KM SHEAR

deepshr=sqrt([x[12]-x[0]]^2+[y[12]-y[0]]^2)
xyouts, -3, 36, deepshr, color=200, charthick=1.25, align=2
shrtxt6='0-6 km Shear:'
xyouts, -11.6, 36, shrtxt6, color=200, charthick=1.25, align=2
unitstxt6='m/s'
xyouts, -4.75, 36, unitstxt6, color=200, charthick=1.25, align=2

;COMPUTE AND DISPLAY 0-1 KM SHEAR

shr1=sqrt([x[2]-x[0]]^2+[y[2]-y[0]]^2)
xyouts, -3, 38, shr1, color=200, charthick=1.25, align=2
shrtxt1='0-1 km Shear:'
xyouts, -11.7, 38, shrtxt1, color=200, charthick=1.25, align=2
unitstxt1='m/s'
xyouts, -4.75, 38, unitstxt1, color=200, charthick=1.25, align=2

;CALCULATE, PLOT, AND DISPLAY LOW-LEVEL STORM-RELATIVE WINDS

sra=[[x[0],y[0]],[motion]]
srb=[[x[1],y[1]],[motion]]
src=[[x[2],y[2]],[motion]]
srd=[[x[3],y[3]],[motion]]
sre=[[x[4],y[4]],[motion]]
srf=[[x[5],y[5]],[motion]]
srg=[[x[6],y[6]],[motion]]

plots, sra, p=0, color=200
plots, srb, p=0, color=200
plots, src, p=0, color=200
plots, srd, p=0, color=200
plots, sre, p=0, color=200
plots, srf, p=0, color=200
plots, srg, p=0, color=200

sr1km=sqrt((x[2]-umotion)^2+(y[2]-vmotion)^2)
sr5km=sqrt((x[10]-umotion)^2+(y[10]-vmotion)^2)
sr10km=sqrt((x[20]-umotion)^2+(y[20]-vmotion)^2)

xyouts, -2.9, 26, sr1km, color=200, charthick=1.25, align=2
srtxt1='1 km Storm-Relative Wind:'
xyouts, -1, 26, srtxt1, color=200, charthick=1.25, align=2
xyouts, -4, 26, unitstxt1, color=200, charthick=1.25, align=2

xyouts, -2.9, 24, sr5km, color=200, charthick=1.25, align=2
srtxt5='5 km Storm-Relative Wind:'
xyouts, -1, 24, srtxt5, color=200, charthick=1.25, align=2
xyouts, -4, 24, unitstxt1, color=200, charthick=1.25, align=2

xyouts, -2.3, 22, sr10km, color=200, charthick=1.25, align=2
srtxt10='10 km Storm-Relative Wind:'
xyouts, 0, 22, srtxt10, color=200, charthick=1.25, align=2
xyouts, -4, 22, unitstxt1, color=200, charthick=1.25, align=2

;polyfill, [x[0],x[1],umotion], [y[0],y[1],vmotion], line_fill=6, color=100

;CALCULATE AND DISPLAY 0-1 KM AND 0-3 KM SRH

;COMPUTE 0-1 KM SRH
b=0
FOR i=0, 1 DO BEGIN
srh=[(x[i+1]-umotion)*(y[i]-vmotion)-(x[i]-umotion)*(y[i+1]-vmotion)] + b
b=srh
srhlow=b
ENDFOR

;COMPUTE 0-3 KM SRH
b=0
FOR i=0, 4 DO BEGIN
srh=[(x[i+1]-umotion)*(y[i]-vmotion)-(x[i]-umotion)*(y[i+1]-vmotion)] + b
b=srh
srh3=b
ENDFOR

xyouts, -3, 32, srhlow, color=200, charthick=1.25, align=2
srhlowtxt='0-1 km SRH:'
xyouts, -13.1, 32, srhlowtxt, color=200, charthick=1.25, align=2
unitstxtsrh='m2 s-2'
xyouts, -1, 32, unitstxtsrh, color=200, charthick=1.25, align=2
xyouts, -3, 30, srh3, color=200, charthick=1.25, align=2
srhtxt='0-3 km SRH:'
xyouts, -13.1, 30, srhtxt, color=200, charthick=1.25, align=2
xyouts, -1, 30, unitstxtsrh, color=200, charthick=1.25, align=2

;CALCULATE BULK RICHARDSON NUMBER SHEAR

upwmw1=(x[0]*1.52+x[1]*1.43+x[2]*1.34+x[3]*1.25+x[4]*1.18+x[5]*1.10+x[6]*1.04+x[7]*.97 + $
		x[8]*.91+x[9]*.86+x[10]*.8+x[11]*.75+x[12]*.71)/13
vpwmw1=(y[0]*1.52+y[1]*1.43+y[2]*1.34+y[3]*1.25+y[4]*1.18+y[5]*1.10+y[6]*1.04+y[7]*.97 + $
		y[8]*.91+y[9]*.86+y[10]*.8+y[11]*.75+y[12]*.71)/13
upwmw2=(x[0]*1.03+x[1]*.97)/2
vpwmw2=(y[0]*1.03+y[1]*.97)/2
brnshr=0.5*((upwmw1-upwmw2)^2+(vpwmw1-vpwmw2)^2)
print, brnshr

;STORM-TYPE FORECAST

;DERIVE 0-3 KM TOTAL SHEAR
b=0
FOR i=0, 5 DO BEGIN
tshr=sqrt((x[i+1]-x[i])^2+(y[i+1]-y[i])^2) + b
b=tshr
tshr3=b
ENDFOR

;DERIVE 0-6 KM TOTAL SHEAR
b=0
FOR i=0, 9 DO BEGIN
tshr=sqrt((x[i+1]-x[i])^2+(y[i+1]-y[i])^2) + b
b=tshr
tshr6=b
ENDFOR

;The junk below is optional...

stormtxt1='Storm Type:'
xyouts, -14.4, 18, stormtxt1, color=200, charthick=1.25, align=2

tshr6by3=(tshr3/tshr6)

IF (shr1 gt 14 AND deepshr gt 10 AND tshr3 gt 35 AND tshr6 lt 35) THEN BEGIN
	stormtype='/ BOW ECHO'
	xyouts, .75, 18, stormtype, color=200, charthick=1.25, align=2
ENDIF

SUP=((shr1/3)+(deepshr/15)+(srh3/150)+(tshr3/20)+(tshr6/20))

IF (SUP gt 6) THEN BEGIN
	stormtype='SUPERCELL'
	xyouts, -6.35, 18, stormtype, color=200, charthick=1.25, align=2
ENDIF

IF (SUP lt 6) THEN BEGIN
	stormtype='MULTICELL'
	xyouts, -4, 18, stormtype, color=200, charthick=1.25, align=2
ENDIF

;TORNADO FORECAST

tortxt1='Tornado Potential:'
xyouts, -9, 16, tortxt1, color=200, charthick=1.25, align=2

TOR=((shr1/4)+(deepshr/18)+(srhlow/50)+(tshr3/19))

IF (TOR gt 5) THEN BEGIN
	torpotential='TORNADO POSSIBLE'
	xyouts, 3, 16, torpotential, color=200, charthick=1.25, align=2
ENDIF

IF (TOR lt 5) THEN BEGIN
	notorpotential='TORNADOES UNLIKELY'
	xyouts, 4.65, 16, notorpotential, color=200, charthick=1.25, align=2
ENDIF

IF (TOR gt 7) THEN BEGIN
	sigtorpotential='"SIG"'
	xyouts, -11.9, 16, sigtorpotential, color=200, charthick=1.25, align=2
ENDIF

END