%% **** Program noxselox_b.m
%% ****
%% **** 
%% ... Read input data - E matrix, F matrix  
%%clear;clc;
E = dlmread('element.dat');
F = dlmread('flow.dat');
%% ... Compounds
i = 0;
fid = fopen ('compounds.dat');
while ( ~ feof(fid))
    i  = i + 1;
    C{i,1} = fgetl (fid);
end % while
fclose(fid);
%% ... check input data (matrix dimensions), tr - trasch values (not used)
[tr,nc ] = size(F);
[tmp,ne] = size(E);
[ncc,tr] = size(C);
if (nc ~= tmp)
    fprintf('Dimensions of matrix F=%i and element matrix E=%i differs\n',nc,tmp);
    return
elseif(nc ~= ncc)
    fprintf('Dimensions of matrix F=%i and compound matrix C=%i differs\n',nc,ncc);
    return
end
%%.. nk = number of knowns is given by
nk = nc - ne;
%%.. nu = number of unknowns is given by
nu = ne;
%% ... "allocate" F1 and F2 arrayfor faster program run
FU = zeros(nu,1); % Delta unknowns fluxes
FK = zeros(nk,1); % Delta knowns fluxes
EU = zeros(nu,ne); % element matrix of unknowns
EK = zeros(nk,ne); % element matrix of knowns
unk = zeros(nu,1); % index of unknowns
knw = zeros(nk,1); % index of knowns
iou = zeros(nu,1); % input or output unknown variable
iok = zeros(nu,1); % input or output known variable

%% ... locate unknown fluxes ( fi < 0)
j = 0;
k = 0;
for i = 1 : nc
    if (F(1,i) < 0 )      % unknown inlet
        j       = j + 1;
        F(1,i)  = 0;
        EU(j,:) = E(i,:);
        unk(j)  = i;
        iou(j)  = 0;
    elseif (F(2,i) < 0 )  % unknown outlet
        j       = j + 1;
        F(2,i)  = 0;
        EU(j,:) = E(i,:);
        unk(j)  = i;
        iou(j)  = 1;
    else
        k       = k + 1;
        FK(k)   = F(2,i) - F(1,i); %% F(out) - F(in)
        EK(k,:) = E(i,:);
        knw(k)  = i;
        
    end
end
%% ... test if ther is no more unknowns then knowns
if (k ~= nk)
    fprintf('\n Number of unknowns (%i) differs from the number of equations (%i)\n',j,nk);
    return
end
%% ... solve invese matrix (E2^T)^(-1) - not efficient!
%FU = - (inv(EU') * (EK') * FK);
%% ... Solve linear system by LU decomposition
% A  = EU';
% b  = -EK'*FK;
% FU = A \ b
%% ... or one line
FU = EU'\(-EK'*FK);
%% ... Print output
fprintf('\n');
fprintf('%10s\t%9s %9s\n','','Inlet','Outlet');

%fprintf('\n Knowns variables\n');
for i = 1: nk
    fprintf('FK(%6s)\t%9.3f %9.3f\n',C{knw(i),1}(:), F(2,knw(i))-FK(i),F(1,knw(i))+FK(i));
end
%fprintf('\n Unknowns variables\n')
for i = 1: nu
    if(iou(i) == 0)  % unknown inlet 
        fprintf('FU(%6s)\t%9.3f %9.3f\n',C{unk(i),1}(:), F(2,unk(i))-FU(i),F(2,unk(i)));
    else             % unknown outlet  
        fprintf('FU(%6s)\t%9.3f %9.3f\n',C{unk(i),1}(:), F(1,unk(i)),F(1,unk(i))+FU(i));
    end
end