portata=84 B=30 ifondo=0.02 ks=40 m=2 k=((84/30)^2/9.8)^(1/3) hu=(portata/B/ifondo^0.5/ks)^(3/5) Evalle=1.5*k+m hvalle=seq(0.01,hu,0.01) Ehvalle=hvalle+portata^2/2/9.81/B^2/hvalle^2 diffEhvalle=Ehvalle-Evalle indmin=which.min(abs(diffEhvalle)) hvalle=hvalle[indmin] h=seq(hvalle,0.65,0.01) E=h+portata^2/(B*h)^2/2/9.8 hm=c(h+0.005) J=(portata/(B*ks*hm^(5/3)))^2 DE=diff(E) Ds=DE/(ifondo-J[1:(length(J)-1)]) s=cumsum(Ds) s=c(0,s,200) h=c(h,hu) plot(s,h,type="l") Emonte=1.5*k+m hmonte=seq(k,6,0.01) Ehmonte=hmonte+portata^2/2/9.81/B^2/hmonte^2 diffEhmonte=Ehmonte-Emonte indmin=which.min(abs(diffEhmonte)) hmonte=hmonte[indmin] h1=seq(hmonte,k,-0.01) E1=h1+portata^2/(B*h1)^2/2/9.8 hm1=c(h1-0.005) J1=(portata/(B*ks*hm1^(5/3)))^2 DE1=diff(E1) Ds1=DE1/(ifondo-J1[1:(length(J1)-1)]) s1=cumsum(Ds1) s1=c(0,s1) plot(s1,h1,type="l") mins1=min(s1) maxs=max(s) plot((s1+200),(h1+4-0.02*200-s1*0.02),type="l",ylim=c(0,(max(h1)+2)),xlim=c(0,200),lwd=3,col="blue") lines(s,h+4-s*0.02,lwd=3,col="blue") stot=c((s1+200),s) hfondo=4-stot*0.02 lines(stot,hfondo,col="brown",lwd=2) lines(stot,hfondo+k,col="red",lty="dashed") lines(stot,hfondo+hu,col="green",lty="dashed") spintam=0.5*9810*h1^2+1000*portata^2/B/h1/B spintav=0.5*9810*h^2+1000*portata^2/B/h/B plot((s1+200),spintam,type="l",ylim=c(0,max(c(max(spintam),max(spintav)))),xlim=c(0,200),lwd=3,col="black") lines(s,spintav,lwd=3,col="red")