portata=84
B=30
ifondo=0.02
ks=40
k=((84/30)^2/9.8)^(1/3)
hu=(portata/B/ifondo^0.5/ks)^(3/5)
h=seq(0.30,0.65,0.01)
E=h+portata^2/(B*h)^2/2/9.8
hm=c(h+0.005,hu)
J=(portata/(B*ks*hm^(5/3)))^2
DE=diff(E)
Ds=DE/(ifondo-J[1:35])
s=cumsum(Ds)
s=c(0,s)
plot(s,h,type="l")