% S-Funktion zur Modellierung des Ruehrkessels aus dem
% GMA-Beispiel (K.-U. Klatt: 5/92 )
% Hierbei gilt:   x(1) = ca
%                 x(2) = cb
%                 x(3) = Temp
%                 x(4) = Tk
%                 u(1) = ca0
%                 u(2) = f = F/V
%                 u(3) = Qk
%                 u(4) = Tzu
%
% Änderung (JOT : 7/94 )
%
%	Kr=[k1  EA1  dH1
%           k2  EA2  dH2
%           k3  EA3  dH3  ]
%z.B.:   
%     Kr=[
%        1.287e12  -9758.3   4.2;
%        1.287e12  -9758.3  -11.0;
%        9.043e09  -8560    -41.85;]
%
%    Konstante=[
%        rho*Cp = 0.9342*3.01;
%        kw*Ar  = 4032.0*0.215;
%        1/(mk*Cpk) =0.1;
%        Vr  = 10;
%    ]

% LACIP/UFRGS: JOT 


function [sys, x0] = cstr2(t,x,u,flag,Kr,Konstante,inicio)


  if abs(flag) == 1
  % Rueckgabe der Zustandsaenderungen

     r=Kr(:,1).*exp(Kr(:,2)./(x(3)+273.15)).*[x(1); x(2);  x(1)^2];
     sys(1) = u(2)*(u(1) - x(1)) - r(1) - r(3);
     sys(2) = -u(2)*x(2) + r(1) - r(2);
     sys(3) = u(2)*(u(4)-x(3)) + (Konstante(2)*(x(4)-x(3))/Konstante(4) - r'*Kr(:,3))/Konstante(1);
     sys(4) = Konstante(3)*(u(3) + Konstante(2)*(x(3)-x(4)));
  
    %   dca = u(2).*(u(1) - x(1)) - k1*x(1) - k3*(x(1).^2);
    %   dcb = -u(2).*x(2) + k1*x(1) - k2*x(2);
    %   dT  = u(2).*(Tzu-x(3)) - (k1*Hab*x(1)+k2*Hbc*x(2)+k3*Had*(x(1).^2))/(rho*Cp) + kw*Ar*(x(4)-x(3))/(rho*Cp*Vr);
    %   dTk = 0.1*(u(3) + kw*Ar*(x(3)-x(4)));

    %

  elseif flag == 3
  % Rueckgabe der Ausgangsgroessen
      sys=x;
  
  elseif flag == 0
  % Rueckgabe der Anfangswerte und der Systemparameter
    sys = [4;0;4;4;0;0];
    x0  = [1.235;0.9;134.14;128.95];
    %x0  = fsolve('cstr_x',x0,[],[],Kr,Konstante,inicio);
	 options = optimset('LargeScale','off','MaxFunEvals',1000000); %vrs 5.3
    x0 = fsolve('cstr_x',x0,options,Kr,Konstante,inicio);% vrs 5.3   
 else
    sys = [];

  end
