clear all
clc
%definição das variáveis simbólicas
syms Ca Cb T F Tk k10 k20 k30 DH1 DH2 DH3 V rho Cp Kw Ar T0 Ca0 E1_R E2_R E3_R mk cpk Qk

%equações diferenciais
df(1) = (F/V)*(Ca0 - Ca) - k10.*exp(-E1_R/T)*Ca - k30.*exp(-E3_R/T)*(Ca^2);
df(2) = -(F/V)*Cb +  k10.*exp(-E1_R/T)*Ca - k20.*exp(-E2_R/T)*Cb;
df(3) = (1/(rho*Cp))*(k10.*exp(-E1_R/T)*Ca*DH1 + k20.*exp(-E2_R/T)*Cb*DH2 + k30.*exp(-E3_R/T)*Ca*Ca*DH3) +  (F/V)*(T0 - T) + ((Kw*Ar)/(rho*Cp*V))*(Tk - T);
df(4) = (1/(mk*cpk))*(Qk+Kw*Ar*(T-Tk)); 

%matriz A do espaço de estados

%primeira linha: derivadas em relação a primeira equação
A(1,1) = diff(df(1),Ca);
A(1,2) = diff(df(1),Cb);
A(1,3) = diff(df(1),T);
A(1,4) = diff(df(1),Tk);

%segunda linha: derivadas em relação a segunda equação
A(2,1) = diff(df(2),Ca);
A(2,2) = diff(df(2),Cb);
A(2,3) = diff(df(2),T);
A(2,4) = diff(df(2),Tk);

%terceira linha: derivadas em relação a terceira equação
A(3,1) = diff(df(3),Ca);
A(3,2) = diff(df(3),Cb);
A(3,3) = diff(df(3),T);
A(3,4) = diff(df(3),Tk);

%terceira linha: derivadas em relação a terceira equação
A(4,1) = diff(df(4),Ca);
A(4,2) = diff(df(4),Cb);
A(4,3) = diff(df(4),T);
A(4,4) = diff(df(4),Tk);

%Matriz B
B(1,1) = diff(df(1),F);
B(1,2) = diff(df(1),Qk);

B(2,1) = diff(df(2),F);
B(2,2) = diff(df(2),Qk);

B(3,1) = diff(df(3),F);
B(3,2) = diff(df(3),Qk);

B(4,1) = diff(df(4),F);
B(4,2) = diff(df(4),Qk);


% % %Calculando os valores de estado estacionário para u7 = 65 h-1
% % %F = 650;
% % F = 1500;
% % %F = 100;
% % 
% % V = 10;
% % Tk = 128.95 + 273.15;
% % y0 = [1 4 Tk];
% % rpar = [F Tk];
% % estados_est = fsolve('vanvusse',y0,[],rpar);
% % Ca = estados_est(1);
% % Cb = estados_est(2);
% % T =  estados_est(3);
% % 
% % %parâmetros
% % 
% % rho = 0.9342;
% % Cp = 3.01;
% % Kw = 4032;
% % Ar = 0.215;
% % T0 = 130 + 273.15;
% % Ca0 = 5.1;
% % DH1 = -4.2;
% % DH2 = 11;
% % DH3 = 41.85;
% % k10 = 1.287e12;
% % k20 = 1.287e12;
% % k30 = 9.043e9;
% % E1_R = 9758.3;
% % E2_R = 9758.3;
% % E3_R = 8560.0;
% % 
% % %valores numéricos para as matrizes A B C e D
% % An = eval(A);
% % Bn = eval(B);
% % Cn = eye(3);
% % Dn = zeros(size(B));
% % 
% % %modelo em espaço de estados
% % modelo_EE = ss(An,Bn,Cn,0);
% % 
% % %funções de transferência para a entrada 1: F
% % [num1,den1]=ss2tf(An,Bn,Cn,Dn,1);
% % 
% % %funções de transferÊncia para a entrada 2: Tk
% % 
% % [num2,den2]=ss2tf(An,Bn,Cn,Dn,2);






