function hdif=sepen(s,tmp,in,out,t)
% sepen - General purpose separator module with energy balances
%         and multiple inlets and outlets.
%
% function hdif=sepen(s,tmp,in,out,t)
% Arguments:
%      s       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).
%      t       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.
%  
% sepen returns the sum of the enthalpies of the exit streams
%     minus the sum of the enthalpies of the inlet streams.
%
% Example:
%   >> sepen(['vvv',[0,55,100],1:2,3:5,[.25 .5 .75; .5 .25 0]);
%
%     this simulates a module where streams 1 & 2 are separated into
%     streams 3-5, with respective temperatures 0, 55, and 100 Tdeg,
%     all being vapor.  The 3 compounds in the product are divided 
%     with compound A going to outlets in a ratio of 25:50:50; B with
%     a ratio of 50:25:25; and C with a ratio of 75:0:25.
%
% JWD, SHD & OKB; help revised by MAH


global ne
nstr=length(out);
if nstr==1
% If there is only one exit stream, use mixe.
   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=mixe(s,tmp,in,out)
   return
end
if nargin<=4
   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 fifth 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
nout=sum([ne(in,4:(nc+3));zeros(1,nc)]);
k=find(abs(nout)>eps);
% This is the beginning of a patch to fix sepen, so it does the 
% separation correctly
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));

