% This file solves the boiler problem in section 5.2;
%
tt = cputime;
 global m Vl TK Patm;

m = 20; Vl = 100;	% initial conditions
TK = 380; Patm = 8;	% initial guess of T(K) and P(atm)
 
Q = 20*10^3;
lambda = 853;
Cpl = 1;
Tr = 32;
rho = 62.4;
F1 = 35.8;
Kv = 6.03;
P0 = 150;
T1 = 70;

					    
	c0 = [m Vl]';
	[t, c] = ode15s('model51', [0, 0.4], c0);
time = cputime-tt

	N = length(c);

% prepare results for plotting

for I = 1:N,
   Vl = c(I,2);
   m = c(I,1);
   mplot(I) = m;
   x = fsolve('eq51',[TK,Patm]');
   TK = x(1);
   Patm = x(2);
   TF = (x(1)-273.15)*1.8+32;
   TFplot(I) = TF;
   Ppsia = x(2)*14.7;
   Ppsiaplot(I) = Ppsia;
   Vlplot(I) = Vl;
   v(I) = (Q + rho*F1*Cpl*(T1-Tr))/(lambda + Cpl*(TF-Tr));
   if Ppsia < P0
     ve(I) = 0;
    else
     ve(I) = Kv*sqrt(Ppsia*(Ppsia - P0));
    end
end


plot(t, Vlplot, '-', t, Ppsiaplot, '--', t, v, '*');
axis([0 0.4 80 160]);
title('Boiler Dynamic Response (Example 5.2)');
xlabel('Time (hr)');
text(0.24, 150, 'P(psia)');
text(0.24, 115, 'Vl(ft^3)');
text(0.24, 95, 'v(lbs/hr)');
print -deps fig56.eps;
pause;
plot(t, mplot);
title('Vapor Mass Dynamic Response')
xlabel('Time (hr)');
ylabel('Vapor Mass, m (lbs)');
print -deps fig57.eps;
pause;
plot(t, TFplot)
title('Temperature Dynamic Response')
xlabel('Time (hr)');
ylabel('Temperature (F)');
print -deps fig58.eps;
pause;
plot(t, ve);
title('Exit Vapor Rate Dynamic Response');
xlabel('Time (hr)');
ylabel('Exit Vapor Rate (lbs/hr)');
print -deps fig59.eps;
