set xrange[0:pi/2]
set yrange[0:2.3]
set size square

set xlabel "{/Symbol q}  [radian]"
set ylabel "phase velocity"
set label "v_A > a" at 0.2,0.2 font "Arial,15"

set xtics ('0' 0, '{/Symbol p}/6' pi/6, '{/Symbol p}/3' pi/3, '{/Symbol p}/2' pi/2)
set ytics ('a' 1, 'v_A' sqrt(2), '(v_A^2+a^2)^{1/2}' sqrt(3))


a=1
VA=sqrt(2)*a

plot 1 with lines dt (10,20) lc rgb "gray" notitle, sqrt(2) with lines dt (10,20) lc rgb "gray" notitle, sqrt(3) with lines dt (10,20) lc rgb "gray" notitle

replot VA*cos(x) title "v_p (shear Alfven wave)" lc rgb "red" linewidth 5
replot sqrt(0.5*(VA**2+a**2+sqrt((VA**2+a**2)**2-4*VA**2*a**2*cos(x)**2))) title "v_{pf} (fast mode)" lc rgb "dark-green" linewidth 5
replot sqrt(0.5*(VA**2+a**2-sqrt((VA**2+a**2)**2-4*VA**2*a**2*cos(x)**2))) title "v_{ps} (slow mode)" lc rgb "blue" linewidth 5

set terminal png
set output "PV1.png"
replot