%
%                   Covariance_construct.m 
%
%  Code to construct a covariance function for the random field associated with the Helmholtz
%  energy
%
%  \alpha(x,\omega) = \alpha_1(\omega) x^2 + \alpha_{11} x^4 + \alpha_{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 = 5000;                    % 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
  
  alpha_1 = -389.4;              % Mean value of alpha_1
  alpha_11 = 761.3;              % Mean value of alpha_11    
  alpha_111 = 61.46;             % Mean value of alpha_111

  per = 0.2;
  alpha1_sample = (alpha_1 - per*alpha_1) + 2*per*alpha_1*rand(Nmc,1);
  alpha11_sample = (alpha_11 - per*alpha_11) + 2*per*alpha_111*rand(Nmc,1);
  alpha111_sample = (alpha_111 - per*alpha_111) + 2*per*alpha_111*rand(Nmc,1);

%
% Plot the Helmholtz energy psi.
%
  P = 0:.05:Pf;
  psi = alpha_1*P.^2 + alpha_11*P.^4 + alpha_111*P.^6;

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

%
%  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 = alpha1_sample*(x.^2) + alpha11_sample*(x.^4) + alpha111_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

  [V,D] = eig(A);
  Phi = W2inv*V;
  
%
% Plot the normalized eigenvalues
%

  figure(2)
  eigenvalue = sort(abs(diag(D)),'descend');
  semilogy(eigenvalue/eigenvalue(1),'x','linewidth',3);
  set(gca,'Fontsize',[20]);
  xlabel('n')
  ylabel('\lambda_n/\lambda_1')
  
%
%  Construct the first three KL modes Q_n(\omega)
% 

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



%
%