set term svg size 800,600 dynamic enhanced

set termoption dashed

set style line 1 lt 1 lc rgb "black" lw 1
set style line 2 lt 2 lc rgb "black" lw 1
set style line 3 lt 3 lc rgb "black" lw 1
set style line 4 lt 4 lc rgb "black" lw 1
set style line 5 lt 5 lc rgb "black" lw 1
set style line 6 lt 1 lc rgb "black" lw 2

max=1.0
ismax(x) = (x>max)?max:x

# x = xh/zbr
# xw = xbrw/zbr

set xlabel 'Barrier Downwind Distance / Barrier Height'
set ylabel 'effective porosity'
set ylabel 'Fraction Open Field Friction Velocity'

# increase effective porosity with barrier width
# pb = pbr + (1 - exp(-0.5*xw))*0.3*(1-pbr)

# calculate coef. as fn of effective porosity
#pb = 1
#a = 0.008-0.17*pb+0.17*pb**1.05
#b = 1.35*exp(-0.5*pb**0.2)
#c = 10*(1-0.5*pb)
#d = 3 - pb

# calc. frac. of fric. vel.
#fu(x) = 1 - exp(-a*x**2) + b*exp(-0.003*(x+c)**d)

a(pb) = 0.008-0.17*pb+0.17*pb**1.05
b(pb) = 1.35*exp(-0.5*pb**0.2)
c(pb) = 10*(1-0.5*pb)
d(pb) = 3 - pb

# calc. frac. of fric. vel.
fu(x,y) = 1 - exp(-a(y)*x**2) + b(y)*exp(-0.003*(x+c(y))**d(y))

#     Cap fu at 1.0
#      if (fu > 1.0) then
#          fu = 1.0

set xrange [-5.4:35]
#set yrange [0:1]
set key bottom

plot ismax(fu(x,1.0)) ls 1 with lines title 'eff. porosity = 1.0' \
   , ismax(fu(x,0.8)) ls 2 with lines title 'eff. porosity = 0.8' \
   , ismax(fu(x,0.6)) ls 3 with lines title 'eff. porosity = 0.6' \
   , ismax(fu(x,0.4)) ls 4 with lines title 'eff. porosity = 0.4' \
   , ismax(fu(x,0.2)) ls 5 with lines title 'eff. porosity = 0.2' \
   , ismax(fu(x,0.0)) ls 6 with lines title 'eff. porosity = 0.0'

#pause -1