disp(' '); disp(' Start here'); disp(' ')
clear
echo on
% This is a MATLAB script to study error analysis of two types,
% absolute errors and root-mean square errors.
% The script assumes that lift coefficient is to be computed
% from measured values of weight (W), air density (rho),
% speed (v), and wing area (S). 
%
% CL=W/(.5*pho*v^2S)
%
% x(1)=W=weight in pounds
% x(2)=rho=air density in slugs per ft^3
% x(3)=v=speed in ft/sec
% x(4)=S=wing area in ft^2
%
% Unfortunately, the measured values of 
% these four quantities are all slightly in error. These errors mean that the 
% computed values for lift coefficient will also be in error.
% This script studies the relationship between the errors
% in the four inputs (W, rho, v S) and the error in the output (CL).
echo off

%W=  x(1);
%rho=x(2);
%v=  x(3);
%S=  x(4);
% Define the 'true' values of weight, density, speed and wing area
Wtrue=5500         % Weight pounds
rhotrue=0.0023769  % Sea level density, slugs per cubic foot
vKnots=70      % speed knots
vtrue=1.687808*vKnots % convert knots to ft/sec
Strue=199.2        % ft squared

% Set things up in vector form where W, rho, v, and S are placed on one vector
x0=[Wtrue rhotrue vtrue Strue];
CLtrue=CLfun(x0)    % Compute the true value of lift coefficient using the function CLfun

error=[7 .00001 .1 .1]   % Assume each x0 has errors of this size
% Perform an error analysis using function err_analysis and function CLfun
[abs_error,rss_error,err_vec]=err_analysis('CLfun',x0,error)

% Let's examine the contribution to the error in CL from each of the individual errors
ifig=0;
ifig=ifig+1; figure(ifig)
bar(err_vec)
title('Individual Errors')
text2(.5,.95,'Errors for W(lbf)   rho(s/f3)   v(f/s)   S(ft)')
xlabel(' Error number 1=W(lbf)   2=rho(s/f3)  3= v(f/s)   4=S(ft)')
ylabel('Error in CL')

echo on
% Analysis using absolute errors
% Assume that  weight, density, speed and wing area are all
% measured with errors.
% Assume that the errors are uniformly distributed around
% the true values plus and minus the absolute error.
% The vector error is interpreted as the +- (absolute) errors 
% of a uniformly distributed random variable.
echo off

np=100; % number of test points
Wuniform=uniform(np,error(1));
rhouniform=uniform(np,error(2));
vuniform=uniform(np,error(3));
Suniform=uniform(np,error(4));
ifig=ifig+1; figure(ifig); plot(Wuniform)
str1=['Uniformly distributed errors in weight, abs error= ',num2str(error(1)),' pounds'];
title(str1); ylabel('Weight error (lbf)')
str2=['Notice that no error is larger then than +- ',num2str(error(1)),' pounds'];
text2(.1,.95,str2)
xlabel('Data point number')

Wvector=Wtrue+Wuniform;
rhovector=rhotrue+rhouniform;
vvector=vtrue+vuniform;
Svector=Strue+Suniform;

abs_error=zeros(size(Wvector));
rss_error=zeros(size(Wvector));
err_vec=zeros(length(Wvector),4);
CLvector=zeros(size(Wvector));

for i=1:length(Wvector)
	x0=[Wvector(i) rhovector(i) vvector(i) Svector(i)];
	CLvector(i)=CLfun(x0);
	[abs_error(i),rss_error(i),err_vec(i,:)]=err_analysis('CLfun',x0,error);
end

ifig=ifig+1; figure(ifig)
titlestring='CL computed using multiple data points with absolute errors. Error bars are +- 1 abs error';
ylabelstring='CL';
ErrorPlot(CLvector,abs_error,1,titlestring,ylabelstring)
text2(.1,.95,'Notice that no errors are outside of the absolute error bars')
hold on; plot([1 length(CLvector)],[CLtrue CLtrue],'-'); text(np/2,CLtrue,'true value'); hold off

echo on
% Analysis using root-mean square (Gaussian) errors
% Assume that  weight, density, speed and wing area are all
% measured with errors.
% Assume that the errors are Gausian distributed around
% the true values with a sigma as specified with the array error.
% The vector error is interpreted as the sigma or roor-mean square of a 
% Gausian distributed random variable.
echo off

WGaussian=error(1)*randn(1,np);
rhoGaussian=error(2)*randn(1,np);
vGaussian=error(3)*randn(1,np);
SGaussian=error(4)*randn(1,np);
ifig=ifig+1; figure(ifig); plot(WGaussian)
str1=['Gausian distributed errors in weight, root-mean square or sigma error= ',num2str(error(1)),' pounds'];
title(str1); ylabel('Weight error (lbf)')
str2=['Notice that many errors (32%) are larger then than +- ',num2str(error(1)),' pounds (the rss error bounds)'];
text2(.03,.95,str2)
str3=['Notice that some errors (5%) are larger then than +- ',num2str(2*error(1)),' pounds (twice the rss error bounds)'];
text2(.03,.85,str3)
xlabel('Data point number')

Wvector=Wtrue+WGaussian;
rhovector=rhotrue+rhoGaussian;
vvector=vtrue+vGaussian;
Svector=Strue+SGaussian;

abs_error=zeros(size(Wvector));
rss_error=zeros(size(Wvector));
err_vec=zeros(length(Wvector),4);
CLvector=zeros(size(Wvector));

for i=1:length(Wvector)
	x0=[Wvector(i) rhovector(i) vvector(i) Svector(i)];
	CLvector(i)=CLfun(x0);
	[abs_error(i),rss_error(i),err_vec(i,:)]=err_analysis('CLfun',x0,error);
end

% Plot the results with +- 1 sigma bounds
ifig=ifig+1; figure(ifig)
titlestring='CL computed using multiple data points with root-mean square errors. Error bars are +- 1sigma';
ylabelstring='CL';
ErrorPlot(CLvector,rss_error,1,titlestring,ylabelstring)
str2=['Notice that many errors (~32%) are larger then than +- ',num2str(rss_error(1)),' (the rss error bounds)'];
text2(.1,.95,str2)
hold on; plot([1 length(CLvector)],[CLtrue CLtrue],'-'); text(np/2,CLtrue,'true value'); hold off

% Plot the results with +- 2 sigma bounds
ifig=ifig+1; figure(ifig)
titlestring='CL computed using multiple data points with root-mean square errors. Error bars are +- 2sigma';
ylabelstring='CL';
ErrorPlot(CLvector,rss_error,2,titlestring,ylabelstring)
str3=['Notice that some errors (~5%) are larger then than +- ',num2str(2*rss_error(1)),'  (twice the rss error bounds)'];
text2(.03,.95,str3)
hold on; plot([1 length(CLvector)],[CLtrue CLtrue],'-'); text(np/2,CLtrue,'true value'); hold off





