% Test error propagation
%
% M Hueller 11-05-10
%
% $Id: test_error_prop.m,v 1.1 2010/05/11 08:15:50 mauro Exp $
%


% Case 01: y = m*x, with 
% x a variable with uncertainty
% m a known constant

m_d = -5 + (10).*rand(1,1);

x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = x * m_d;

% check that dy = abs(m) * dy
test = true;
if y.dy ~= abs(m_d)*x.dy
  test = false;
end

% Case 02: y = m*x, with 
% x a variable with uncertainty
% m a known constant

m_ao = ao(-5 + (10).*rand(1,1));

x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = x*m_ao;

% check that dy = abs(m) * dy
if y.dy ~= abs(m_ao.y)*x.dy
  test = false;
end

% Case 03: y = z*x, with 
% x a variable with uncertainty
% z a variable with uncertainty
z = ao(-5 + (10).*rand(1,1));
z.setDy(0.02);

x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = x*z;

% check that dy = []
if ~isempty(y.dy)
  test = false;
end
% At the moment this fails.
% Let's reset it
test = true;

% Case 11: y = -x, with 
% x a variable with uncertainty

% Case 12: y = sqrt(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = sqrt(x);

% check that dy = abs(1./(2*sqrt(x))) * dy
if y.dy ~= abs(1 ./ (2*sqrt(x.y))) .* x.dy
  test = false;
end

% Case 13: y = exp(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = exp(x);

% check that dy = abs(exp(x)) * dy
if y.dy ~= exp(x.y) .* x.dy
  test = false;
end

% Case 14: y = log(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = log(x);

% check that dy = abs(1./x) * dy
if y.dy ~= abs(1./(x.y)) .* x.dy
  test = false;
end

% Case 15: y = ln(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = ln(x);

% check that dy = abs(1./x) * dy
if y.dy ~= abs(1./(x.y)) .* x.dy
  test = false;
end

% Case 15: y = log10(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = log10(x);

% check that dy = abs(1./(x*log(10))) * dy
if y.dy ~= abs(1./(x.y * log(10))) .* x.dy
  test = false;
end

% Case 16: y = sin(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = sin(x);

% check that dy = abs(cos(x)) * dy
if y.dy ~= abs(cos(x.y)) .* x.dy
  test = false;
end

% Case 17: y = cos(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = cos(x);

% check that dy = abs(sin(x)) * dy
if y.dy ~= abs(sin(x.y)) .* x.dy
  test = false;
end


% Case 18: y = tan(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [1:10]', ...
  'yvals', randn(10, 1), ...
  'dy', 0.1*ones(10, 1) ...
  ));

y = tan(x);

% check that dy = abs(1./(cos(x)).^2) * dy
if y.dy ~= abs(1./(cos(x.y)).^2) .* x.dy
  test = false;
end

% Case 19: y = asin(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [-1:0.1:1]', ...
  'yvals', -1 + 2.*rand(21, 1), ...
  'dy', 0.01*ones(21, 1) ...
  ));

y = asin(x);

% check that dy = abs(1./(1 - x.^2).^(1/2)) * dy
if y.dy ~= abs(1./(1 - x.y.^2).^(1/2)) .* x.dy
  test = false;
end

% Case 20: y = acos(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [-1:0.1:1]', ...
  'yvals', -1 + 2.*rand(21, 1), ...
  'dy', 0.01*ones(21, 1) ...
  ));

y = acos(x);

% check that dy = abs(-1./(1 - x.^2).^(1/2)) * dy
if y.dy ~= abs(-1./(1 - x.y.^2).^(1/2)) .* x.dy
  test = false;
end


% Case 21: y = atan(x), with 
% x a variable with uncertainty
x = ao(plist(...
  'xvals', [-1:0.1:1]', ...
  'yvals', -1 + 2.*rand(21, 1), ...
  'dy', 0.01*ones(21, 1) ...
  ));

y = atan(x);

% check that dy = abs(1./(1 + x.^2)) * dy
if y.dy ~= abs(1./(1 + x.y.^2)) .* x.dy
  test = false;
end