function [Tjs,Tfs,concs] = cstr_io(Tempvec)
%
% Generate the steady-state i/o curve for the diabatic CSTR.
% Tempvec is a vector of reactor temperatures generated
% before this routine is called.
%
% This routine solves for vector of jacket temperatures
% that would yield the required reactor temperature vector
% (assuming a constant feed temperature).
%
% In addition, this routine solves for the vector of
% feed temperatures that would yield the required reactor
% temperature vector (assuming constant jacket temperature).
% 
% For consistency with cstr_ss.m and cstr_dyn.m, we use
% a global parameter vector. The elements of this vector 
% should be set in the command window before this routine
% is run.
%   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)
%
% 16 July 96 - (c) B.W. Bequette
% revised 2 January 1997 - added CSTR_PAR
%
  global CSTR_PAR
%
% 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;
%
for i = 1:length(Tempvec);
  Temp   = Tempvec(i);
  k      = k0*exp(-E_act/(R*Temp));
  ca     = fov*caf/(fov + k);
  rate   = k0*exp(-E_act/(R*Temp))*ca;
  stuff1 = fov*rhocp*(Temp - Tf);
  stuff2 = H_rxn*rate;
  Tjs(i) = Tempvec(i) + (stuff1 - stuff2)/UAoV;
  concs(i) = ca;
  stuff3 = UAoV*(Temp - Tj);
  stuff4 = H_rxn*rate;
  Tfs(i) = Tempvec(i) + (stuff3 - stuff4)/(fov*rhocp);
end;