%% Bilance simultannich reakci
%% --------------------------------------------------------------
%% nc     .... (in)  number of components 
%% nr     .... (in)  number of reactions
%% mnu    .... (in)  stoi. matrix, dimension  [nr,nc] 
%% n0     .... (in)  vector of initial composition [mol/time] 
%% outmix .... (out) vector of equilibrium composition [mol/time] 
%% kx     .... (in)  vector of equilibrium constants, dimension [nr] 
%% alfa   .... (in-out) init. guess and resulting value of conversion 
%%                      of key component j 
%% key    .... (in)  ID of key component <1;nc>
%% -------------------------------------------------------------- ....
%% Clean arguments
%%fprintf ('%s', program_name ());
%%arg_list = argv ();
%%for i = 1:nargin
%%  fprintf ( '%s', arg_list{i});
%%end %for
%%printf ('\n');
%%
%% ......... Global variables ...............
global key nr nc mnu n0 p kx
%%
%%............ read input data ..............
%%
fid=fopen('input.dat','r');
%% Initial comment "%"
trash=fgetl(fid);
%% Read "nc", "nr", "j" a "p"
nx=fscanf(fid,'%i\n',Inf);
%% Put vector valuest to the scalar variables
%% avoids future probles ...
nc = nx(1);
nr = nx(2);
key= nx(3);
%% Comment line
trash=fgetl(fid);
%% Read pressure
p=fscanf(fid,'%f\n',Inf);
%% Comment line
trash=fgetl(fid);
%% Read Id. of compounds -  (!UWAGA!)
idl = fgetl(fid);
ii=1;
while (any(idl))
  [temp,idl] = strtok(idl,' ');
  id{ii} = temp;
  ii=ii+1;
end   %while
%%idl = fgetl(fid);
%%id  = cellstr(split(idl," ",nc));
%% Comment
trash=fgetl(fid);
%% Read vector of eq. constants
kx=fscanf(fid,'%f\n',Inf);
%% Commnet
trash=fgetl(fid);
%% Read vector of init. concentrations
n0=fscanf(fid,'%f\n',Inf);
%% Comment
trash=fgetl(fid);
%% Read vector of init. reaction extent
ksi0=fscanf(fid,'%f\n',Inf);
%% Comment
trash=fgetl(fid);
%% Read stoi. matrix
tmp=fscanf(fid,'%f',[nc,nr]);
fclose(fid);
for i=1:nr
  for j=1:nc
    mnu(i,j)=tmp(j,i);
  end    %for
end   %for
tol = zeros(nr,1);
%%.... Rank of stoi. matrix
rankmnu=rank(mnu);
if (nr > rankmnu)
  fprintf('Rank of stoi. matrix is less then NR\n')
else
  fprintf('Rank of stoi. matrix is %i\n',rankmnu)
end %if
%%
%%................ Solution ..............
%%
wrgstoi = zeros(nr,1);
prtrcell(nr, nc, mnu, id, wrgstoi)

for i = 1 : nr
    tol(i) = 1e-6;
end
[x, eps, info] = newrap(@fkxi, ksi0, tol);
fprintf('\n')
fprintf('Total number of newrap iteration = %i\n',info)
for i = 1 : nr
    fprintf('xi(%i) = %10.4e with final eps = %10.6e \n', i,x(i),eps(i));
end
%% .... Solve system of NEQ
%opt=optimset('MaxIter',1000000,'TolX',1e-7,'TolFun',1e-7);
%[x, info, msg]=fsolve(@fkxi,ksi0,opt);
%if ( msg == 1)
%    fprintf('The fsolve routine converge normally and found:\n');
%    for i = 1 : nr
%        fprintf('Reaction extent (reaction %i) = %10.4e\n',i,x(i));
%    end
%else
%    fprintf('The fsolve routine could not converge\n');
%    fprintf('Error message code %i\n',msg);
%end

%% .... Init. molar frac
x0 = (n0/sum(n0))';
%% .... Outlet mol frac
xq = fcompm(nr,nc,mnu,n0,x);
fprintf('\n');
for i=1:nc
  fprintf('%s\t', id{i}(:))
end %for
fprintf('\n')
for i=1:nc
  fprintf('%6.4f\t', x0(i))
end
fprintf('\n')
for i=1:nc
  fprintf('%6.4f\t', xq(i))
end
fprintf('\n')
fprintf('Conversion of key component %s\n',id{key}(:))
