function [c,C,bo,cor,X,r,cond] = solsys(A,b,tol,X,C)
%  function [c,C,e95,cor,X,r,cond] = solsys(A,b,tol,X,C)
%
% Least squares solution of an over determined system Ac = b,
% where A (m,n), c(n,1) and b(m,1). 
% The method of solution is through the Pseudo-inverse of A, obtained by
% SVD with tolerance tol. The LS coefficients (c), the covariance matrix
% of the coefficients (C) and the 95% confidence limits on the fit
% (e95), are obtained as output. If tol=0, tol=max(size(A))*s(1)*eps 
% is used.
%
%                         Julio Candela VII/93.
 igr = 0;
 [m,n] = size(A); b=b(:);
%
% calculate Pseudo-inverse of A.
%
 if nargin == 3,
    [X,V,S,r,cond] = pinvvs(A,tol,igr);
    C = V(:,1:r) * (S.*S) * (V(:,1:r))';
 if igr ;  clg ; 
 subplot(211),contour(C), title('covariance matrix of coefficients'),
 subplot(212), mesh(C), subplot, end
 end,
%
% solved for coefficients and calculate statistics of fit.
%
 c = X * b; bo = A*c ; d = b - bo ; v = d'*d;
%
% calculate integral (~cumulative sum) of autocorrelation function
% to obtained unbiased variance estimate. Also calculates correlation
% coefficient of the fit.
%
%% cb = xcov(b,'coeff'); cb = cb(m:m*2-1); indp = max(cumsum(cb(1:m/4)));
%% dof = (m/indp) - n; var = v/dof;
 corre = correl(real(b),real(bo)); 
% cori = correl(imag(b),imag(bo)); if isnan(cori), cori=0; end;
 cori=0;
 cor = corre + i*cori; 
%
% if 'b' is assume to be a random variable normally distributed about
% a true value 'bm' with varince 'v', the coefficients c will also be 
% random variables with variance given by:
%
% cv = v.*diag(C);
%
% Because the coefficients are random variables the estimates 'bo' are 
% also randomly distributed about a true value 'bom' with variance:
%
% bov = v .* diag(A*C*A');
% 
% the 95% confidence limits on true mean value at the observation
% locations are given considering a value of t = 2.086 for dof=20.
%
% t = 2.086; e95 = zeros(m,1); s = sqrt(var); st = s*t;
% e95 = s * sqrt(diag(A*C*A'));
% ss = ['independence every ',num2str(indp),' points'],
% ss = ['dof = ',num2str(dof),'  unbiased std =',num2str(s)],
 
