function f = cstr_ss(x)
%
% 2-state diabatic CSTR model
% Use fsolve to find steady-state concentration and temperature
%
% The call statement from the command window is:
%
%   x = fsolve('cstr_ss',x0);
% 
% where x is the solution vector (x = [conc;temp])
% and x0 is an initial guess vector
%
% define the following global parameter so that values
% can be entered in the command window
%     global CSTR_PAR
% parameters (case 2)
% CSTR_PAR(1):  frequency factor (9703*3600)
% CSTR_PAR(2):  heat of reaction (5960)
% CSTR_PAR(3):  activation energy (11843)
% CSTR_PAR(4):  density*heat cap. (500)
% CSTR_PAR(5):  heat trans coeff * area (150)
% CSTR_PAR(6):  ideal gas constant (1.987)
% CSTR_PAR(7):  reactor volume (1)
% CSTR_PAR(8):  feed flowrate (1)
% CSTR_PAR(9):  feed concentration (10)
% CSTR_PAR(10): feed temperature (298)
% CSTR_PAR(11): jacket temperature (298)
%
% states 1 and 2
% ca     = concentration of A, kgmol/m^3
% Temp   = reactor temperature, deg K
%
% 16 July 96 - (c) B.W. Bequette
% revised 2 January 1997 - added CSTR_PAR
%
  f = zeros(2,1);
%
  global CSTR_PAR
%
% below we use common notation for the states to
% improve code readibility
%
  ca     = x(1);
  Temp   = x(2);
%
% parameter values, case 2 shown in parentheses:
%
  k0    = CSTR_PAR(1);  % frequency factor (9703*3600)
  H_rxn = CSTR_PAR(2);  % heat of reaction (5960)
  E_act = CSTR_PAR(3);  % activation energy (11843)
  rhocp = CSTR_PAR(4);  % density*heat cap. (500)
  UA    = CSTR_PAR(5);  % ht trans coeff * area (150)
  R     = CSTR_PAR(6);  % gas constant (1.987)
  V     = CSTR_PAR(7);  % reactor volume (1)
  F     = CSTR_PAR(8);  % feed flowrate (1)
  caf   = CSTR_PAR(9);  % feed concentration (10)
  Tf    = CSTR_PAR(10); % feed temperature (298)
  Tj    = CSTR_PAR(11); % jacket temperature (298)
%
% ratios used in the equations
%
  fov   = F/V;
  UAoV  = UA/V;
%
% modeling equations:
%
  rate = k0*exp(-E_act/(R*Temp))*ca;
%
  dcadt = fov*(caf - ca) - rate;
  dTdt  = fov*(Tf - Temp) + (H_rxn/rhocp)*rate - (UAoV/rhocp)*(Temp - Tj);
%
% below we convert back to function notation.  fsolve
% seeks values of x(1) and x(2) to drive f(1) and f(2) to zero.
%
  f(1) = dcadt;
  f(2) = dTdt;