%% Physics 233, MATLAB Exercises for Hughes and Hase: Chapter 4
% MATLAB activities to complement _Measurements and their Uncertainties_,
% Chapter 4, "Error propagation." Section numbers in the book
% correspond to those listed below.
%
% Author: J. S. Dodge, 2016

%% 4.1 Propagating the error in a single-variable function
% Reproduce Fig. 4.1(a). The caption says that this shows Bragg's law for a
% silicon single crystal. This is incorrect: silicon has a room temperature
% lattice constant of $a = 0.5431$ nm, which gives $\lambda = 2a = 1.1$ nm
% for $\theta = \pi/2$, as shown in the plot. But silicon's diamond
% structure _suppresses_ Bragg scattering from planes with this spacing;
% Bragg scattering _is_ allowed, however, for lattice planes $d =
% a/\sqrt{3} = 0.3$ nm, which would give $\lambda = 2a/\surd{2} = 0.63$ nm
% for $\theta = \pi/2$. We'll ignore this oversight so that we can compare
% results with the book.

d = 0.5431;  % (nm)
thetaRad = linspace(0,pi/2);    % MATLAB uses radians for angles
thetaDeg = rad2deg(thetaRad);   % Convert to degrees for plotting

plot(thetaDeg, 2*d*sin(thetaRad), 'k-')
xlabel('Angle (degrees)')
xlim([0 90])
ylabel('\lambda (nm)')
box off
text(0.95,0.8,'(a)',...
    'Units','Normalized')

%% 4.1.1 The functional approach for single-variable functions
% Consider a fixed angular uncertainty of $\alpha_\theta = 0.050^\circ$.

alphaThetaDeg = 0.05;
alphaThetaRad = deg2rad(alphaThetaDeg);

%%
% Confirm the wavelength uncertainties given in the text.

thetaEvalDeg = [15, 75];
thetaEvalRad = deg2rad(thetaEvalDeg);

alphaLambda = 2*d*(sin(thetaEvalRad+alphaThetaRad) - sin(thetaEvalRad));

disp(1e3*alphaLambda)   % Display in units of pm

%%
% Reproduce Fig. 4.1(b) for $\bar{\theta} = 75^\circ$, where the |sin|
% function is more nonlinear. Use the |hold on| and |hold off| commands to
% make MATLAB overlay the results of different plot commands. This makes it
% easier to control the formatting of the solid line, the dashed lines, the
% open circles, and the closed circle. Graphical elements in MATLAB are
% represented as objects, with properties that control their appearance and
% behavior. The command |gca| stands for "get current axes," and returns a
% pointer to the current "Axes" object. Use the |get| command to get the   
% y-axis limits of the plot, given by the |YLim| property of the current
% axes.

thetaBarDeg = thetaEvalDeg(2);
thetaBarRad = deg2rad(thetaBarDeg);
thetaIntervalDeg = thetaBarDeg+alphaThetaDeg*[-1,1];
thetaIntervalRad = deg2rad(thetaIntervalDeg);
thetaLim = thetaBarDeg + 2.5*alphaThetaDeg*[-1, 1];
thetaDeg = linspace(thetaLim(1),thetaLim(2));
thetaRad = deg2rad(thetaDeg);

% Plot theory curve
figure
plot(thetaDeg, 2*d*sin(thetaRad), 'k-')
xlim(thetaLim)

% Get y-axis limits
lambdaLim = get(gca,'YLim');

% Hold plot
hold on

% Place closed marker at measurement point
plot(thetaBarDeg,2*d*sin(thetaBarRad),'ko',...
    'MarkerFaceColor','k');

% Place open markers at uncertainty limits
plot(thetaIntervalDeg,2*d*sin(thetaIntervalRad),'ko')

% Connect markers with axes with dashed lines
markerPosX = [thetaBarDeg, thetaIntervalDeg];
markerPosY = 2*d*sin([thetaBarRad, thetaIntervalRad]);

for iMarker = 1:3
vLineX = markerPosX(iMarker)*[1,1];
vLineY = [lambdaLim(1),markerPosY(iMarker)];
hLineX = [thetaLim(1),markerPosX(iMarker)];
hLineY = markerPosY(iMarker)*[1,1];
plot(hLineX,hLineY,'k--',...
    vLineX,vLineY,'k--')
end

hold off

box off
text(0.95,0.8,'(b)',...
    'Units','Normalized')

%% 4.1.2 A calculus-based approximation for single-variable functions
% Compare the calculus-based method with the functional approach for the
% Bragg's law example.

alphaLambdaCalc = abs(2*d*cos(thetaEvalRad))*alphaThetaRad;

disp('Uncertainty, functional approach (pm):')
disp(1e3*alphaLambda)
disp('Uncertainty, calculus-based approximation (pm):')
disp(1e3*alphaLambdaCalc)

%%
% Increase angular uncertainty by a factor of 100, from 0.05 degrees to 5
% degrees, and compare again. Note that the calculus-based estimate
% overestimates the uncertainty because the |sin| function is
% sub-linear---that is, as the argument increases, the slope decreases.

alphaThetaDeg = 100*0.05;
alphaThetaRad = deg2rad(alphaThetaDeg);

alphaLambda = 2*d*(sin(thetaEvalRad+alphaThetaRad) - sin(thetaEvalRad));
alphaLambdaCalc = abs(2*d*cos(thetaEvalRad))*alphaThetaRad;

disp('Uncertainty, functional approach (pm):')
disp(1e3*alphaLambda)
disp('Uncertainty, calculus-based approximation (pm):')
disp(1e3*alphaLambdaCalc)

%% 4.1.3 Look-up table for common single-variable functions
% Let's continue this comparison for the functions shown in Table 4.1. Set
% $A = 1$ and $\alpha_A = 0.01$, and calculate $\alpha_Z$ using both the
% functional approach and the calculus-based approach. It's convenient to
% put this in a loop, where the function changes in each iteration but
% otherwise the calculation is the same. We can do this using two new
% MATLAB programming concepts: structures and anonymous functions.

%%
% Let's start with <http://www.mathworks.com/help/matlab/structures.html
% structures> , which are a way to store data of different types in a
% single array. We'll begin by making a structure array |S| that contains
% the function names in the field |name|. 

S(1).name = '1/A';
S(2).name = 'exp(A)';
S(3).name = 'ln A';
S(4).name = 'log A';
S(5).name = 'A^n';
S(6).name = '10^A';
S(7).name = 'sin A';
S(8).name = 'cos A';
S(9).name = 'tan A';

disp(S)

disp(S(1).name)

%%
% Note that MATLAB displays structures differently than arrays. You can
% also open |S| in the variable browser to see the contents of each
% structure array element and field.
%
% MATLAB provides several ways to define functions, as described in the
% <http://www.mathworks.com/help/matlab/function-basics.html
% documentation>. To define simple functions like those in Table 4.1, we'll
% use
% <http://www.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html
% anonymous functions>. Let's start by defining the function $Z(A) = 1/A$.

zFun = @(a) 1./a;
a = (1:5)';
disp([a zFun(a)])

%%
% Now we can define new fields in |S|, |function| and |derivative|, that
% contain anonymous functions for the function and its derivative. We
% choose $n=3$ for the function $Z(A) = A^n$.

S(1).function = @(a) 1./a;
S(2).function = @(a) exp(a);
S(3).function = @(a) log(a);
S(4).function = @(a) log10(a);
S(5).function = @(a) a.^3;
S(6).function = @(a) 10.^a;
S(7).function = @(a) sin(a);
S(8).function = @(a) cos(a);
S(9).function = @(a) tan(a);

S(1).derivative = @(a) -1./a.^2;
S(2).derivative = @(a) exp(a);
S(3).derivative = @(a) 1./a;
S(4).derivative = @(a) 1./(log(10)*a);
S(5).derivative = @(a) 3*a.^2;
S(6).derivative = @(a) 10.^a*log(10);
S(7).derivative = @(a) cos(a);
S(8).derivative = @(a) -sin(a);
S(9).derivative = @(a) 1 + tan(a).^2;

%%
% The structure now has nine elements with three fields each, and two of
% the fields contain anonymous functions.

disp(S)
disp(S(1).function)
disp(S(1).derivative)

%%
% Now we can loop over the elements of the structure |S| to compare the
% functional and calculus-based approaches to error propagation methods. To
% put it in a nice tabular form, we use another new function,
% <http://www.mathworks.com/help/matlab/ref/fprintf.html |fprintf|>,
% which is designed to write formatted data to files but can also direct
% its output to the MATLAB command window. The output is a string that may
% include special strings called
% <http://www.mathworks.com/help/matlab/ref/fprintf.html#inputarg_formatSpec
% _formatting operators_>, which start with a percent sign, |%|.


A = 1;
alphaA = 0.01;
fprintf('A = %.4f, dA = %.4f\n',A,alphaA);

fprintf('Z(A)           Z(A+dA)-Z(A)   |dZ/dA|alphaA\n')
for iZ = 1:length(S)
    alphaZfun = abs(S(iZ).function(A+alphaA) - S(iZ).function(A));
    alphaZcal = abs(S(iZ).derivative(A))*alphaA;
    fprintf('%-8s       %9.6f      %9.6f\n',...
        S(iZ).name,alphaZfun,alphaZcal);
end

%%
% Let's take a moment to examine this data and how it is generated.
%
% First, we set |A = 1| and |dA = 0.01|, and display this with the
% |fprintf| function. By calling |fprintf| without an initial |fileID|
% parameter, we direct the output to the command window. The first argument
% to |fprintf| is the format string |'A = %.4f, dA = %.4f\n'|, which
% includes two formatting operators that convert the the next two
% arguments, |A| and |dA|, into formatted strings. The format operator
% |%.4f| will convert a number into a floating-point number with four
% digits of precision. The format string ends with |\n|, which is a newline
% string that indicates the end of the current line and the start of a new
% one.
%
% The next group of statements creates the table. The initial call to
% |fprintf| creates the header line, which is just a string and a newline
% (|\n|) without any additional formatting. After the header, we have a
% |for| loop that creates the data table.  For each element of the
% structure array |S|, the statements inside the loop compute $\alpha_Z$
% using the functional approach, then compute it using the calculus-based
% approach, then prints the function name and both values for $\alpha_Z$ to
% the command window using the format string |'%-8s       %9.6f
% %9.6f\n'|: the first format operator, |%-8s|, formats a string input
% (|S(iZ).name| here) so that it is left-justified within a field eight
% characters wide; the second and third format operators, |%9.6f|, format
% numerical inputs (|alphaZfun| and |alphaZcal|) as floating-point numbers
% with six digits of precision and nine characters wide. The white space
% between the format operators ensures that the tabular information lines
% up with the header columns, and the whole format string ends with a
% newline (|\n|) so that each iteration of the loop creates a new line in
% the table.
%
% With minor modifications, we can use the same commands to direct the
% formatted text to a file. To do this, we need to use the MATLAB functions
% <http://www.mathworks.com/help/matlab/ref/fopen.html |fopen|> and
% <http://www.mathworks.com/help/matlab/ref/fclose.html |fclose|>. To
% create a new file for reading and writing, we use the syntax |fileID =
% fopen(filename,permission)|, where the
% <http://www.mathworks.com/help/matlab/ref/fopen.html#inputarg_permission
% |permission|> string is set to |'w+'|. Note that this will overwrite any
% data that already exists in a file with the same name, so be careful when
% you use this statement! The |fopen| statement returns a file identifier,
% |fID|; by including this as the first argument in an |fprintf| statement,
% you direct its output to the file. Once you are finished with the file,
% you can close it with |fclose(fID)|.
%
% The lines below increase |alphaA| by a factor of ten, to |alphaA = 0.1|,
% and write the resulting error propagation values to |TestTable.txt|. Once
% you run these commands, you can view the results in a text editor;
% alternatively, you can view the file in the command window using the
% <http://www.mathworks.com/help/matlab/ref/type.html |type|> statement.

fID = fopen('TestTable.txt','w+');

A = 1;
alphaA = 0.1;
fprintf(fID, 'A = %.4f, dA = %.4f\n',A,alphaA);

fprintf(fID, 'Z(A)           Z(A+dA)-Z(A)   |dZ/dA|alphaA\n');
for iZ = 1:length(S)
    alphaZfun = abs(S(iZ).function(A+alphaA) - S(iZ).function(A));
    alphaZcal = abs(S(iZ).derivative(A))*alphaA;
    fprintf(fID, '%-8s       %9.6f      %9.6f\n',...
        S(iZ).name,alphaZfun,alphaZcal);
end

fclose(fID);

type('TestTable.txt')

%% 4.1.4 Worked example---single variable function
% Check the worked examples using |S(6)|, defined above. The book discusses
% the functional and calculus-based approach for |alphaA = 0.1| and
% |alphaA = 0.4|. The |%5.1g| format operator displays the input in
% floating-point format to one ''significant figure'', as appropriate for
% uncertainties (|%5.1f| will produce one digit after the decimal place,
% regardless of its significance).

A = 2.3;
Zbar = S(6).function(A);
fprintf('Zbar = %5.1f\n', Zbar);

alphaA = 0.1;

alphaZfunP = abs(S(6).function(A+alphaA) - S(6).function(A));
alphaZfunM = abs(S(6).function(A) - S(6).function(A-alphaA));
alphaZcal = abs(S(6).derivative(A))*alphaA;

fprintf('Error, functional approach (alphaA = %.1f)\n',alphaA);
fprintf('alphaZ = (%5.1g, %5.1g)\n\n',...
    -alphaZfunM,alphaZfunP);
fprintf('Error, calculus-based approach (alphaA = %.1f)\n',alphaA);
fprintf('alphaZ = (%5.1g, %5.1g)\n\n',...
    -alphaZcal,alphaZcal);

alphaA = 0.4;

alphaZfunP = abs(S(6).function(A+alphaA) - S(6).function(A));
alphaZfunM = abs(S(6).function(A) - S(6).function(A-alphaA));
alphaZcal = abs(S(6).derivative(A))*alphaA;

fprintf('Error, functional approach (alphaA = %.1f)\n',alphaA);
fprintf('alphaZ = (%5.1g, %5.1g)\n\n',...
    -alphaZfunM,alphaZfunP);
fprintf('Error, calculus-based approach (alphaA = %.1f)\n',alphaA);
fprintf('alphaZ = (%5.1g, %5.1g)\n\n',...
    -alphaZcal,alphaZcal);

%% 4.2 Propagating the error through a multi-variable function
%% 4.2.1 The functional approach for multi-variable functions
%% 4.2.2 Worked example---functional approach for multi-variable functions
% Confirm the results from this section, using an anonymous function of
% ''two'' variables. In this example, note the scope of the variables: a,
% b, and R are not defined as function variables, so the function inherits
% the values assigned to them in the
% <http://www.mathworks.com/help/matlab/matlab_prog/base-and-function-workspaces.html
% base workspace>. Note also that the function is
% <http://www.mathworks.com/help/matlab/matlab_prog/vectorization.html
% vectorized> using element-wise
% <http://www.mathworks.com/help/matlab/ref/times.html multiplication> and 
% <http://www.mathworks.com/help/matlab/ref/rdivide.html division>; these
% are both examples of a larger class of
% <http://www.mathworks.com/help/matlab/matlab_prog/array-vs-matrix-operations.html
% array operators>. Finally, note the unit definitions included beside all
% of the constants, and that the anonymous function includes multiplication
% by $10^{-6}$ to convert the units from Pa to MPa.

VmBar = 2.000e-4;       % m^3/mol
alphaV = 0.003e-4;

TBar = 298.0;           % K
alphaT = 0.2;

a = 1.408e-1;           % (m^6 Pa)/mol^2
b = 3.913e-5;           % m^3/mol
R = 8.3145;             % J/K = (m^3 Pa)/K

Pfun = @(Vm,T) 1e-6*(R*T./(Vm - b) - a./Vm.^2);    % MPa

fprintf('P(VmBar,Tbar) = %.3f MPa\n', Pfun(VmBar,TBar));

%%
% Vectorizing the function can simplify the code. Here, we form the vectors
% $[\bar{V}_m, \bar{V}_m + \alpha_V]$ and $[\bar{T}, \bar{T} + \alpha_T]$
% and use the |diff| function to find differences in the function
% evaluations.

VmRange = VmBar + [0, 1]*alphaV;
TRange = TBar + [0, 1]*alphaT;

alphaPT = abs(diff(Pfun(VmBar,TRange)));
alphaPV = abs(diff(Pfun(VmRange,TBar)));
alphaP = sqrt(alphaPT^2 + alphaPV^2);

fprintf(...
    'alphaPT = %.3f MPa, \t alphaPV = %.3f MPa, \t alphaP = %.3f MPa\n',...
    alphaPT,alphaPV,alphaP);

%% 4.2.3 A calculus approximation for multi-variable functions
%% 4.2.4 A look-up table for multi-variable functions
% Practice your multivariable error propagation skills by using Eq. 4.16 to
% reproduce each line in Table 4.2 by direct calculation.

%% 4.2.5 Comparison of methods
% Confirm the results of this section. Note that the error is given as a
% percentage of the measured value.

Zfun = @(A,B) (A-B)./(A+B);

Abar = 1000;    % Dimensionless
Bbar = 80;

alphaA = 0.01*Abar;
alphaB = 0.01*Bbar;

Arange = Abar + [0,1]*alphaA;
Brange = Bbar + [0,1]*alphaB;

Zbar = Zfun(Abar,Bbar);

fprintf('Zbar = %.3f\n',Zbar);

dZdA = @(A,B) 2*B./(A+B).^2;
dZdB = @(A,B) -2*A./(A+B).^2;

alphaZcalc = sqrt((dZdA(Abar,Bbar).*alphaA).^2 + ...
    (dZdB(Abar,Bbar).*alphaB).^2);
alphaZfun = sqrt((diff(Zfun(Arange,Bbar))).^2 + ...
    (diff(Zfun(Abar,Brange))).^2);

fprintf('alphaZcalc = %.5f,\t alphaZfun = %.5f\n', alphaZcalc, alphaZfun);

%% 4.2.6 Percentage errors---dominant error
% This is an important idea for any experimentalist: find the ''dominant''
% error and focus on it. Uncertainties add in quadrature: $\alpha_Z =
% \sqrt{\alpha_A^2 + \alpha_B^2}$ for $Z = A \pm B$, so if $\alpha_A =
% 2\alpha_B$ the square operation makes it ''four times'' more important to
% $\alpha_Z$. Confirm this using the example in this section. Let
% |alphaRel| be a vector of the relative uncertainties for four quantities,
% and compute the uncertainty in their product by taking the square root of
% the sum of their squares. Confirm that reducing the first source of
% uncertainty is more important than reducing the second.

alphaRel = [5 1 1 1]*0.01;
alphaTot = sqrt(sum(alphaRel.^2));
fprintf('Initial: %.5f\n',alphaTot)

alphaRel = [5/2 1 1 1]*0.01;
alphaTot = sqrt(sum(alphaRel.^2));
fprintf('Reduce first term: %.5f\n',alphaTot)

alphaRel = [5 1/2 1 1]*0.01;
alphaTot = sqrt(sum(alphaRel.^2));
fprintf('Reduce second term: %.5f\n',alphaTot)

%% 4.2.7 Using the look-up tables
%% 4.2.8 Using the look-up tables---health warning
% We can illustrate the conclusions of this section with a simulation. Use
% the random number generator to simulate the results of |Nmc = 1000|
% measurements of $A$, $B$, and $Z$, with $Z = (A-B)/(A+B)$, and confirm
% that the standard deviation of $Z$ is consistent with $\alpha_Z \approx
% 0.002$ calculated above.

rand('seed',0)
%rng('default')
Nmc = 1000;

Amc = Abar + alphaA*randn(Nmc,1);
Bmc = Bbar + alphaB*randn(Nmc,1);
Zmc = Zfun(Amc,Bmc);
stdZ = std(Zmc);
fprintf('stdZ = %.5f,\t alphaZcalc = %.5f\n', stdZ, alphaZcalc);

%%
% Now plot $A-B$ versus $A+B$ to demonstrate that they are
% correlated. To understand why, note that $\alpha_A = 10$ and $\alpha_B =
% 0.8$, so the scatterplot of $B$ versus $A$ is much narrower along $B$
% than along $A$---set the axis limits to have the same width to show this
% clearly. The plot of $A-B$ versus $A+B$ just rotates the distribution by
% 45 degrees.

figure;
plot(Amc+Bmc, Amc-Bmc, '.')
xlabel('A+B')
ylabel('A-B')

figure;
plot(Amc, Bmc, '.')
xlabel('A')
ylabel('B')
ylim([40 120])

%%
% The <http://www.mathworks.com/help/stats/scatterhist.html |scatterhist|>
% plot function in the Statistics and Machine Learning Toolbox makes the
% point of this section even more clearly.

figure;
scatterhist(Amc+Bmc, Amc-Bmc, 'Marker', '.')
xlabel('A+B')
ylabel('A-B')

%%
% The scatterplots show that large values of $A-B$ are correlated with
% large values of $A+B$; this makes the fluctuations in the ratio $Z
% = (A-B)/(A+B)$ ''smaller'' than what they would be if the numerator and
% denominator were uncorrelated. We can confirm this with another
% simulation: let $Z' = X/Y$, where $\bar{X} = \bar{A}-\bar{B}$, $alpha_X =
% alpha_{A-B}$, $\bar{Y} = \bar{A}+\bar{B}$, $\alpha_Y =
% \alpha_{A+B}$, but $X$ and $Y$ are ''uncorrelated''.

Xbar = Abar - Bbar;
Ybar = Abar + Bbar;
alphaX = sqrt(alphaA^2 + alphaB^2);
alphaY = sqrt(alphaA^2 + alphaB^2);
alphaZprime = Zbar*sqrt((alphaX/Xbar)^2 + (alphaY/Ybar)^2);

Xmc = Xbar + alphaX*randn(Nmc,1);
Ymc = Ybar + alphaY*randn(Nmc,1);
Zprimemc = Xmc./Ymc;

stdZprime = std(Zprimemc);
fprintf('stdZprime = %.5f,\t alphaZprime = %.5f\n', ...
    stdZprime, alphaZprime);

figure;
scatterhist(Xmc,Ymc, 'Marker', '.')
xlabel('X')
ylabel('Y')

figure;
histogram(Zmc,(0.8:0.001:0.9))
hold on
histogram(Zprimemc,(0.8:0.001:0.9))
hold off
xlabel('Z')
ylabel('Frequency')
legend('Z','Z\prime')

%% 4.3 Propagating errors in functions---a summary
%% 4.4 Experimental strategy based on error analysis
%% 4.4.1 Experimental strategy for reducing the dominant error
%% 4.5 Combined experiments---the weighted mean
% Simulate the result of combining two consistent measurements to improve a
% result. In the first experiment, Alice performs five measurements of
% the acceleration of gravity with a pendulum, and reports the mean and
% standard error of her results.

g = 9.8;        % m/s^2
alphag = 0.01;

n1 = 5;
gm1 = g + alphag*randn(n1,1);

gBar1 = mean(gm1);
gErr1 = std(gm1)/sqrt(n1);

fprintf('Alice measures: (%.3f +/- %.3f) m/s^2\n', gBar1, gErr1);

%%
% Now Bob does another five measurements. He also reports the mean and
% standard error of his results.

n2 = 5;
gm2 = g + alphag*randn(n2,1);

gBar2 = mean(gm2);
gErr2 = std(gm2)/sqrt(n2);

fprintf('Bob measures: (%.3f +/- %.3f) m/s^2\n', gBar2, gErr2);

%%
% The standard errors are different, so to combine them we must take the
% weighted mean using (4.30), and use (4.32) to compute the combined
% standard error.

w1 = 1/gErr1^2;
w2 = 1/gErr2^2;

gCE = (w1*gBar1 + w2*gBar2)/(w1+w2);
gErrCE = sqrt(1/(w1+w2));

fprintf('The combined measurement is: (%.3f +/- %.3f) m/s^2\n', ...
    gCE, gErrCE);

%%
% Of course we could also just take the mean and standard error of ''all''
% of their measurements, since we have the raw data---the weighted mean is
% really only necessary if we must rely on the reported means and standard
% errors. The results of the two methods may be slightly different, but
% they should be close if the original measurements are limited only by
% statistical uncertainty.

gAll = mean([gm1;gm2]);
gErrAll = std([gm1;gm2])/sqrt(n1+n2);
fprintf('Combining the raw data gives: (%.3f +/- %.3f) m/s^2\n', ...
    gAll, gErrAll);

%% 4.5.1 The error in the mean---a special case of the weighted mean

%%
% $Id: P233HH3.m 151 2016-02-02 18:38:44Z physcmpt $