% programme provisoire de démo UQSU
% JLBK, 12 février 2008 pour le Groupe Autosurveillance GRAIE
clear p2 yi n m F Q R M b e Sr C;
format short eng;
% lecture données (hi, Si)
disp(' Nom du fichier .csv contenant les données (hi, Si)');
disp(' (ne pas écrire l''extension): ');
NomFichier = input(' ','s');
Origine = input(' Passage du polynôme par l''origine ? o/n : ', 's');
p2=dlmread([NomFichier,'.csv'],';',1,0);
yi = p2(:,2);
n = size(p2,1);
m = 3;
% construction matrice F
if Origine == 'o'
F(:,1)=p2(:,1);
F(:,2)=p2(:,1).^2;
F(:,3)=p2(:,1).^3;
else
F(1:size(p2,1),1)=ones;
F(:,2)=p2(:,1);
F(:,3)=p2(:,1).^2;
F(:,4)=p2(:,1).^3;
end
% décomposition QR de F
[Q,R]=qr(F,0);
M=inv(R)*Q';
b=M*yi;
disp(' ');
disp('b = ');
disp(num2str(b));  
disp(' ');
% calcul des résidus
e=F*b-yi;
% calcul de la somme des carrés des résidus
Sr = e'*e;
disp('Sr = ');
disp(num2str(Sr));  
disp(' ');
% matrice variance/covariances
C=Sr/(n-m-1)*M*M';
disp('C = ');
disp(num2str(C));  
disp(' ');
% calculs de S, Q et u(Q)
d = [1 1 1 1 ];
while d(1) > 0
disp(['valeurs de [h u(h) U u(U)] entre crochets pour le calcul ?' ;
      '(pour arrêter: taper 0 (zéro) pour h)                     ']);
disp(' ');
d = input('');
if d(1) > 0
h = d(1); U = d(3);
if Origine == 'o'
    vh = [h; h^2; h^3];
else
    vh = [1; h; h^2; h^3];
end
S = vh'*b;
Q = S*U;
u1Q2 = d(4)^2*S^2;
if Origine == 'o'
    u2Q2 = d(2)^2*U^2*(b(1)+2*b(2)*h+3*b(3)*h^2)^2;
else
    u2Q2 = d(2)^2*U^2*(b(2)+2*b(3)*h+3*b(4)*h^2)^2;
end
u3Q2 = Sr/(n-m-1)*vh'*M*M'*vh;
uQ = sqrt(u1Q2+u2Q2+u3Q2);
DQsurQ = 2*uQ/Q*100;
disp(' ');
disp(['section mouillée S (m2) :                ',num2str(S)]);
disp(['débit (m3/s) :                           ',num2str(Q)]);
disp(['incertitude type u(Q) (m3/s) :           ',num2str(uQ)]);
disp(['incertitude relative élargie (%) :       ',num2str(DQsurQ)]);
disp(['contributions u(U), u(h) et u(S) en % :  ',...
         num2str(100*u1Q2/uQ^2),' , ',num2str(100*u2Q2/uQ^2),' , ',num2str(100*u3Q2/uQ^2)]);
disp(' ');
end
end
disp fin
