%
%   Code for the Saltelli implementation of Sobol sensitivity indices.
%

  clear all

%
% Construct the matrices for the integral response.
%

  M = 1e+6;
  A = rand(M,2);
  B = rand(M,2);
  C1 = B;
  C2 = B;
  C1(:,1) = A(:,1);
  C2(:,2) = A(:,2);

  for j=1:M
    params = A(j,:);
    y_A(j) = sqrt(params(2)/params(1));
%    y_A(j) = cos(sqrt(params(1)/params(2))*pi/2);
    
    params = B(j,:);
    y_B(j) = sqrt(params(2)/params(1));
%    y_B(j) = cos(sqrt(params(1)/params(2))*pi/2);

    params = C1(j,:);
    y_C1(j) = sqrt(params(2)/params(1));
%    y_C1(j) = cos(sqrt(params(1)/params(2))*pi/2);

    params = C2(j,:);
    y_C2(j) = sqrt(params(2)/params(1));
%    y_C2(j) = cos(sqrt(params(1)/params(2))*pi/2);
  end

  f02 = ((1/M)^2)*sum(y_A)*sum(y_B);

%
% Construct the Sobol indices for m and k.
%

  S1 = ((1/M)*y_A*y_C1' - f02)/((1/M)*y_A*y_A' - f02);
  S2 = ((1/M)*y_A*y_C2' - f02)/((1/M)*y_A*y_A' - f02);

  ST1 = 1 - ((1/M)*y_B*y_C1' - f02)/((1/M)*y_A*y_A' - f02);
  ST2 = 1 - ((1/M)*y_B*y_C2' - f02)/((1/M)*y_A*y_A' - f02);

  S = [S1 S2]
  ST = [ST1 ST2]