%% **** Program noxselox_b.m
%% ****
%% **** 
%% ---- Read input files ----

%% ... Compounds symbols
i = 1;
fid = fopen ('compounds.dat');
while ( ~ feof(fid))
  %%A{i,1} = fgetl (fid);
  A{i,1} = fgetl(fid);
  i  = i + 1; 
end % while
fclose(fid);
%% ... FLOW RATE MATRIX
%%     [kmol\s]
F = dlmread('flow.dat');
%% ... WEIGTHS  
W = dlmread('weights.dat');
%% ... ELEMENT MATRIX
E = dlmread('element.dat');

%%
%% ---- BEGIN executable statements ----
%%
%% ... check sizes
[na,trs]  = size(A);
[nc,trs]  = size(F);
[tmp,ne]  = size(E);
[sw,trs]  = size(W);
%% ...
if (nc ~= tmp)
    fprintf('N. of row of flow matrix F=%i and element matrix E=%i differs\n',nc,tmp);
    return
elseif(sw ~= nc)
    fprintf('N. of elements of vector W=%i is not equal to nc=%i \n',sw0,nc);
    return
elseif(na ~= nc)
    fprintf('N. of elements of vector A=%i is not equal to nc=%i\n',na,nc);
    return
else
    %% ... Allocate arrays to make the code faster
    D  = zeros(nc,nc);
    b  = zeros(ne,1);
    z  = zeros(ne,1);
    cF = zeros(nc,2);
end %if

%% ... Create D matrix and evaluate P
for i = 1 : nc
  D(i,i) = (1./W(i,1) + 1./W(i,2));
end %for 
P = E' * D * E;

%% ... b = 2 sum[e(ij)*(F0(j) - F1(j))]
for i = 1 : ne
  fsum = 0;
  for j = 1 : nc
    fsum = fsum + E(j,i)*(F(j,1) - F(j,2));
  end %for
  b(i) = 2 * fsum;
end %for

%% ... Solving P.z = b
z = P \ b;

%% ... Evaluate cF0 and cF1
for i = 1 : nc
  cF(i,1) = F(i,1) - 0.5*(E(i,:)*z(:)) / W(i,1);
  cF(i,2) = F(i,2) + 0.5*(E(i,:)*z(:)) / W(i,2);
end % for

%%
%% ---- Print results
%%
fprintf('\n           **Measured Values**  | **Corrected Values**\n');
fprintf('          Inlet      Outlet     |    Inlet      Outlet   \n');
fprintf('          -----------------------------------------------\n');
for i = 1 : nc
    fprintf('%8s %10.3e %10.3e  |   %10.3e %10.3e \n', A{i,1}(:), F(i,:), cF(i,:));
end %for
return;
