// REACTOR RCAI NO ISOTÉRMICO - LINEALIZACIÓN

// parámetros
tau	= 1 // hr , V/v
ko	= 9703*3600 // hr-1
mDH	= 5960 // kcal/kmol
Ea	= 11843 // kcal/kmol
rocp	= 500 // kcal/(m3ºC)
Tin	= 25+273 // K
CAin	= 10 // kmol/m3
UA_V	= 150 // kcal/(m3ºC hr)
Tc	= 25+273 // K
R = 1.987 // kcal/(kmol.K)

function f = rcai(y)
  CA = y(1); T = y(2)
  f(1) = (CAin - CA)/tau - ko*exp(-Ea/R./T).*CA
  f(2) = (Tin - T)/tau +(mDH)/rocp*ko*exp(-Ea/R./T).*CA-UA_V/rocp*(T-Tc)
  f=f'
endfunction
//y0 = [9 300]; por ejemplo
CA0 = input('punto inicial de la iteración, CA (kmol/m3) =  ')
T0 = input('punto inicial de la iteración, T (K) =  ')
y0 = [CA0 T0]

y = fsolve(y0,rcai);

CAs = y(1); Ts = y(2)
disp(CAs,'CAs =  '); disp(Ts,'Ts =  ')

ks = ko*exp(-Ea/R/Ts)
ksprima = ks*Ea/R/Ts^2

A(1,1) = -1/tau - ks
A(1,2) = -ksprima*CAs
A(2,1) = ks*mDH/rocp
A(2,2) = -1/tau - UA_V/rocp + ksprima*mDH/rocp*CAs

spec(A); disp(spec(A),'valores propios: ') // calcula valores propios
if real(spec(A))<zeros(2,1)
  printf('punto estable')
else 
  printf('punto inestable')
end

B(1) = 0;
B(2) = UA_V/rocp;

C = eye(2,2);

sl = syslin('c',A,B,C);
g = ss2tf(sl);disp(g,'g(s) = ')
