%
%                   Covariance_construct.m 
%
%  Code to construct a covariance function for the random field associated with the Helmholtz
%  energy
%
%  \alpha(x,\omega) = a_1(\omega) x^2 + a_{11} x^4 + a_{111} x^6
%
%
  clear all
%

%
% Input polarization values and constants.  Construct intervals with per=20% perturbation values.
%

  Pf = 1.0;
  Nquad = 50;                     % Number of composite trapezoid or midpoint quadrature points
  Nmc = 50000;                    % Number of Monte Carlo sampling points
  method = 2;                     % method = 1: Composite trapezoid
                                  % method = 2; Composite midpoint 
%
%
  if method==1
    h = Pf/(Nquad - 1);
    x = 0:h:Pf;
    w = h*ones(1,Nquad);
    w(1) = 0.5*h;
    w(Nquad) = 0.5*h;
  else
    h = Pf/Nquad;
    w = h*ones(1,Nquad);
    for j=1:Nquad
      values(j) = 2*j-1;
      x(j) = (h/2)*values(j);
    end
  end
  
  a_1 = -389.4;              % Mean value of a_1
  a_11 = 761.3;              % Mean value of a_11    
  a_111 = 61.46;             % Mean value of a_111

  per = 0.2;
  a1_sample = (a_1 - per*abs(a_1)) + 2*per*abs(a_1)*rand(Nmc,1);
  a11_sample = (a_11 - per*abs(a_11)) + 2*per*abs(a_11)*rand(Nmc,1);
  a111_sample = (a_111 - per*abs(a_111)) + 2*per*abs(a_111)*rand(Nmc,1);

%
% Plot the Helmholtz energy psi.
%
  P = 0:.05:Pf;
  psi = a1_sample(1,1)*P.^2 + a11_sample(1,1)*P.^4 + a111_sample(1,1)*P.^6;

  figure(1)
  plot(P,psi,P,0*P,'k','linewidth',3)
  set(gca,'Fontsize',[20]);
  xlabel('Polarization x')
  ylabel('Helmholtz Energy \alpha')

%
%  Compute the matrices of random field evaluations and covariance matrix C as well as diagonal matrix W of quadrature weights
%
%    Alpha: Matrix of random field evaluations -- alpha(x_j,omega^k)
%    Alpha_c: Matrix of centered process -- alpha_c(x_i,omega^k) = alpha(x_i,omega^k) - (1/Nmc)sum_{m=1}^Nmc alpha(x_i,omega^m)
%    C: Covariance matrix

  Alpha = a1_sample*(x.^2) + a11_sample*(x.^4) + a111_sample*(x.^6);
  Alpha_c = Alpha - (1/Nmc)*ones(Nmc,1)*sum(Alpha);
  C1 = (1/(Nmc-1))*(Alpha_c')*Alpha_c;    % C1 can be slightly nonsymmetric due to numerical issues
  C = max(C1,C1.');                       % This guarantees that C is symmetric
  W = diag(w);
  W2 = sqrtm(W);
  W2inv = inv(W2);
  A = W2*C*W2;

%
% Solve the symmetric eigenvalue problem and compute eigenvectors

  [V1,D1] = eig(A);               

  D=diag(sort(diag(D1),'descend'));       % The following commands sort D and V from largest to smallest
  [c, ind]=sort(diag(D),'descend'); 
  V = V1(:,ind);
  Phi = W2inv*V;
  eigenvalues = sort(abs(diag(D)),'descend');

%
%  Compute the bound b_N given by (6.21)
%

  b_N = sum(eigenvalues(4:30))/sum(eigenvalues(1:3))

  
%
% Plot the normalized eigenvalues
%

  figure(2)
%  semilogy(eigenvalues/eigenvalues(1),'x','linewidth',3);
  semilogy(eigenvalues(1:30),'+','linewidth',3);
  axis([1 30 1e-16 2000])
  set(gca,'Fontsize',[20]);
  xlabel('n')
%  ylabel('\lambda_n/\lambda_1')
  ylabel('\lambda_n')
  
%
%  Construct the first three KL modes Q_n(\omega)
% 

  clear Alpha Alpha_c

  omega_vals = 0:.01:1;
  N_omega = length(omega_vals);
  a1_s = (a_1 - per*abs(a_1)) + 2*per*abs(a_1)*omega_vals';
  a11_s = (a_11 - per*abs(a_11)) + 2*per*abs(a_11)*omega_vals';
  a111_s = (a_111 - per*abs(a_111)) + 2*per*abs(a_111)*omega_vals';

  Alpha = a1_s*(x.^2) + a11_s*(x.^4) + a111_s*(x.^6);
  Alpha_c = Alpha - (1/N_omega)*ones(N_omega,1)*sum(Alpha);


  for k=1:N_omega
    sum1 = 0;
    sum2 = 0;
    sum3 = 0;
    for j=1:Nquad
      sum1 = sum1 + w(j)*Alpha_c(k,j)*Phi(j,1);
      sum2 = sum2 + w(j)*Alpha_c(k,j)*Phi(j,2);
      sum3 = sum3 + w(j)*Alpha_c(k,j)*Phi(j,3);
    end
    Qn(k,1) = (1/sqrt(eigenvalues(1)))*sum1;
    Qn(k,2) = (1/sqrt(eigenvalues(2)))*sum2;
    Qn(k,3) = (1/sqrt(eigenvalues(3)))*sum3;
  end

  figure(3)
  plot(Qn(:,1))

%
% Compute and plot the covariance function c(x,y)
%

  x_small = x(1:2:end);
  C_small = C(1:2:end,1:2:end);
  figure(3)
  mesh(x_small,x_small,C_small,'linewidth',2) 
  brighten(-.2)
  set(gca,'Fontsize',[20]);
  xlabel('x-axis')
  ylabel('y-axis')
  zlabel('c(x,y)')

  
%
%