%
%  Linear Data Reconciliation
%
%     NV = number of variables,  NE = number of equations
%
%  Data file must contain (NE+3) rows of (NV) columns with the
%  following format:
%
%    Row 1 to NE :     equation coefficients
%    Row NE + 1  :     right hand side plus (NV-NE) zeros
%    Row NE + 2  :     standard deviations of variable errors
%    Row NE + 3  :     measurements of variables
%  
%  Bob Newell, November 1994
%
%--------------------------------------------------------------
disp( ' ' )
disp( 'Linear Data Reconciliation' )
%
%  load the problem data from nominated file
%
file = input( 'enter the data file name (in single quotes) ' );
disp( ' ' )
disp( [ 'loading ' file ' ...' ] )
eval( [ 'load ' file ] )
data = [];
for i = [1 : length(file)];
  if strcmp( file(i), '.' ) == 1; break; end;
  data = [ data file(i) ];
end
eval( [ 'rec = ' data ';' ] );
eval( [ 'clear ' data ] );
%--------------------------------------------------------------
%  find the dimensions of the problem
%
[nr nc] = size( rec );
nv = nc;
disp( ' ' )
disp( [ 'number of variables = ' int2str(nv) ] )
ne = nr - 3;
disp( [ 'number of balances = ' int2str(ne) ] )
input( 'press enter for the balances ...' );
%--------------------------------------------------------------
%  now get the model equations
%
a = rec( [1:ne], : );
b = rec( ne+1, [1:ne] )';
disp( ' ' )
disp( 'the A matrix and b vector are:' )
[ a b ]
input( 'press enter for the measurement data ...' );
%--------------------------------------------------------------
%  now get the measurement stdev and values
%
xs = rec( ne+2, : )';
xm = rec( ne+3, : )';
disp( ' ' )
disp( 'the measurement stdev and values are:' )
[ xs' ; xm']
input( 'press enter for adjustments and error tests ...' );
%--------------------------------------------------------------
%  now for the calculations
%
q = inv( diag( xs .* xs ) );
aa = [ 2*q a' ; a zeros(ne) ];
bb = [ 2*q*xm ; b ];
xx = inv( aa ) * bb;
x = xx( [1:nv] );
disp( ' ' )
disp( 'new_value, adjustment, statistic and gross_error_flag are:' )
e = xm - x;
% z = abs(e) / std(e);
z = abs(e) ./ xs;
[ x e z (z > 1.96) ]
%--------------------------------------------------------------
% yep that's all
