%
% Project 2
%
%                       Problem 2
%
%
  clear all
%
% Input polarization values and constants.  Construct intervals with per=20% perturbation values.
%

  Pf = 0.8;
  P2 = -Pf:.01:Pf;
  P = 0:.05:Pf;
  alpha_1 = -389.4;
  alpha_11 = 761.3;
  alpha_111 = 61.46;

  per = 0.2;   
  Alpha_l = [alpha_1 - per*alpha_1; alpha_11 - per*alpha_11; alpha_111 - per*alpha_111]; 
  Alpha_r = [alpha_1 + per*alpha_1; alpha_11 + per*alpha_11; alpha_111 + per*alpha_111];

  psi = alpha_1*P2.^2 + alpha_11*P2.^4 + alpha_111*P2.^6;

%
% Plot the Helmholtz energy psi.
%

  figure(1)
  plot(P2,psi,P2,0*P2,'k','linewidth',3)
  axis([-Pf Pf -60 80])
  set(gca,'Fontsize',[20]);
  xlabel('Polarization P')
  ylabel('Helmholtz Energy \psi')

%
% Compute the derivatives with respect to alpha_1, alpha_11, alpha_111 and
% and compute the sensitivity matrix and Fisher information matrix V.  Compute
% the rank of V.
%

  psi_alpha_1 = P.^2;
  psi_alpha_11 = P.^4;
  psi_alpha_111 = P.^6;

  sensmat = [psi_alpha_1; psi_alpha_11; psi_alpha_111];
  V = sensmat*sensmat';
  rank(V)

%
%  Compute Morris indices
%

  r = 50;                  % Number of random samples
  delta = 1/20;            % Stepsize
  for i = 1:r
    xs = rand(3,1);
    x1 = xs + delta*[1;0;0];
    x2 = xs + delta*[0;1;0];
    x3 = xs + delta*[0;0;1];
    qs = Alpha_l + (Alpha_r - Alpha_l).*xs;
    q1 = Alpha_l + (Alpha_r - Alpha_l).*x1;
    q2 = Alpha_l + (Alpha_r - Alpha_l).*x2;
    q3 = Alpha_l + (Alpha_r - Alpha_l).*x3;
    diff1 = q1 - qs;
    diff2 = q2 - qs;
    diff3 = q3 - qs;
    psi_s = qs(1)*(1/3)*Pf^3 + qs(2)*(1/5)*Pf^5 + qs(3)*(1/7)*Pf^7;
    psi_1 = q1(1)*(1/3)*Pf^3 + q1(2)*(1/5)*Pf^5 + q1(3)*(1/7)*Pf^7;
    psi_2 = q2(1)*(1/3)*Pf^3 + q2(2)*(1/5)*Pf^5 + q2(3)*(1/7)*Pf^7;
    psi_3 = q3(1)*(1/3)*Pf^3 + q3(2)*(1/5)*Pf^5 + q3(3)*(1/7)*Pf^7;
    d1 = (psi_1 - psi_s)/diff1(1);
    d2 = (psi_2 - psi_s)/diff2(2);
    d3 = (psi_3 - psi_s)/diff3(3);
    d(i,:) = [d1 d2 d3]; 
    absd(i,:) = [abs(d1) abs(d2) abs(d3)];
  end
  
  mu1 = (1/r)*sum(d(:,1));
  mu2 = (1/r)*sum(d(:,2));
  mu3 = (1/r)*sum(d(:,3));

  mus1 = (1/r)*sum(absd(:,1));
  mus2 = (1/r)*sum(absd(:,2));
  mus3 = (1/r)*sum(absd(:,3));

  sigma1 = (1/(r-1))*sum((d(:,1)-mu1).^2);
  sigma2 = (1/(r-1))*sum((d(:,2)-mu2).^2);
  sigma3 = (1/(r-1))*sum((d(:,3)-mu3).^2);

%
% Output the Morris indices mu and sigma.
%

  mu = [mus1 mus2 mus3]
  sigma = [sqrt(sigma1) sqrt(sigma2) sqrt(sigma3)]

%
% Use the Saltelli algorithm to compute Sobol indices.
% 

  M = 10000;            % Number of random samples
  X = rand(M,3);
  Y= rand(M,3);
  Z1 = Y;
  Z2 = Y;
  Z3 = Y;
  Z1(:,1) = X(:,1);
  Z2(:,2) = X(:,2);
  Z3(:,3) = X(:,3);

  for j = 1:M
    A(j,:) = Alpha_l' + X(j,:).*(Alpha_r' - Alpha_l');
    B(j,:) = Alpha_l' + Y(j,:).*(Alpha_r' - Alpha_l');
    C1(j,:) = Alpha_l' + Z1(j,:).*(Alpha_r' - Alpha_l');
    C2(j,:) = Alpha_l' + Z2(j,:).*(Alpha_r' - Alpha_l');
    C3(j,:) = Alpha_l' + Z3(j,:).*(Alpha_r' - Alpha_l');
  end

  y_A = A(:,1)*(1/3)*Pf^3 + A(:,2)*(1/5)*Pf^5 + A(:,3)*(1/7)*Pf^7; 
  y_Aabb = A(:,1)*(1/3)*Pf^3 + A(:,2)*(1/5)*Pf^5;
  y_B = B(:,1)*(1/3)*Pf^3 + B(:,2)*(1/5)*Pf^5 + B(:,3)*(1/7)*Pf^7;
  y_C1 = C1(:,1)*(1/3)*Pf^3 + C1(:,2)*(1/5)*Pf^5 + C1(:,3)*(1/7)*Pf^7;
  y_C2 = C2(:,1)*(1/3)*Pf^3 + C2(:,2)*(1/5)*Pf^5 + C2(:,3)*(1/7)*Pf^7;
  y_C3 = C3(:,1)*(1/3)*Pf^3 + C3(:,2)*(1/5)*Pf^5 + C3(:,3)*(1/7)*Pf^7;
  
  f02 = ((1/M)^2)*sum(y_A)*sum(y_B);

  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);
  S3 = ((1/M)*y_A'*y_C3 - 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);
  ST3 = 1 - ((1/M)*y_B'*y_C3 - f02)/((1/M)*y_A'*y_A - f02);

%
% Output the Sobol indices.
%

  S = [S1 S2 S3]
  ST = [ST1 ST2 ST3]
   
%