% programa 5.2, W.L.Luyben, pag. 118
% "gravity-flow tank simulation" 
% modif. E.Quintero.M 18nov'98
%----------------------------------------------------
% (time is in seconds)
% two cases are run, at time equal zero the feed flow
% rate is increased to 100% of the design flow rate
% (15,700 gpm)
% case no.1 starts from 67% of design flow rate
% case no.2 starts from 50% of design flow rate 
clear
format short e
v1=3.40; v2=2.50; h1=2.05; h2=1.2; xx=[]; yy=[];
for n=1:2
   if(n==1)
      v=v1;
      h=h1;
   else
      v=v2;
      h=h2;
   end
   time=0;
   delta=1;
   tprint=0;
   while(time<=800)
      % evaluate derivatives      
      vdot=0.0107*h-0.00205*v^2;
      hdot=0.311-0.0624*v;
      if(time>=tprint)
         if(n==1) 
            xx=[xx;time,v,h];
         else 
            yy=[yy;time,v,h]; 
         end 
         tprint=tprint+20;
      end
      % integrate using euler
      v=v+vdot*delta;
      h=h+hdot*delta;
      time=time+delta;     
   end
end
xx,yy
plot(xx(:,1),xx(:,3)),axis([0 800 -2 10]),grid
xlabel('process time, s'),ylabel('h, ft')
hold on
plot(yy(:,1),yy(:,3))
