function hdif=reacten(s,tmp,in,out,rs,t)
% reacten - General purpose reactor module with energy balances
%           and multiple inlets and outlets
% 
% reacten is a general purpose reactor module with energy balances.
% The number of inlet and exit streams can be any integer >0.
%
% function hdif=reacten(s,tmp,in,out,rs,t)
% Arguments:
%    s       vector of states of the exit stream(s).
%            ('v', 'l', or 's')
%    tmp     temperatures of exit streams in Tdeg.
%    in      the indices of the inlet stream(s).
%    out     the indices of the exit stream(s).
%    rs      vector of reaction rates.
%    t       vector of the fraction of each compound in the product
%            that goes into each stream, except that the last stream
%            need not be given. Each row should correspond to one
%            exit stream.  If tmp and out are scalars, t does not need
%            to be given.
%
%  reacten returns the sum of the enthalpies of the 
%     exit streams minus the sum of the enthalpies 
%     of the inlet streams.
%
% Example: 
%     >> reacten('vs',[312 25],1:2,3:4,[.94 0.21],[1 1 0]);
%
% JWD, SHD & OKB, TYLC; help revised by MAH


global ne stoic
nstr=length(out);
if nstr==1
% If there is only one exit stream, use reacte.
   if length(tmp)~=1
% Check first to see that there is exactly one exit temperature
      disp('The number of streams out is one, but the length of tmp is:')
      disp(length(tmp))
      return
   end
   hdif=reacte(s,tmp,in,out,rs)
   return
end
if nargin<=5
   disp('Your number of exit streams is:')
   disp(nstr)
   disp('You must give the fraction of each compound out for all except')
   disp('the last of those streams as a sixth argument.')
end
[nrt,nct]=size(t);
if nrt==nstr-1 & nrt>0
% the zeros in the line below are included so that sum sums columns
% when t is a vector and there are two exits.
   t=[t;1-sum([t;zeros(1,nct)])];
end
sne=size(ne);
nc=sne(2)-3;
% these zeros are needed for the same reason as above
NS=sum([ne(in,4:(nc+3));zeros(1,nc)]);
nout=NS+(rs*stoic);
k=find(abs(nout)>eps);
str=ne(in(1),:);
setne(s(1),in(1),tmp(1),nout(k),k)
if nstr>1
   for j=2:nstr
      sepe(s(j),s(1),tmp(j),tmp(1),in(1),out([j 1]),t(j,:));
   end
   rem=nout.*t(1,:);
   l=find(abs(rem)>eps);
   setne(s(1),out(1),tmp(1),rem(l),l);
   ne(in(1),:)=str;
end
hdif=sum(ne(out,3))-sum(ne(in,3));
