// REACTOR RCAI NO ISOTÉRMICO - DIAGRAMA DE FASE

// 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(t,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)
endfunction
clf
plot2d([2.359 8.5636],[368.1 311.2],rect=[0,250,15,400],style=[-2]) //puntos estables
plot2d([5.518],[339.1],rect=[0,250,15,400],style=[-9])  //punto inestable
xlabel('CA (kmol/m3)')
ylabel('T (K)')

t = 0:0.1:20
for i = 1:20
[b,x10,x20] = xclick()
       x = ode([x10;x20],0,t,rcai)
    plot2d(x(1,:),x(2,:),rect = [0,250,15,400])
end
