function [Phi,decayrates,B,Tmean] = solvePODModes(self,snapshots,coefstruct)
% solvePODModes - Internal function to compute POD modes.
%

% Copyright 2021-2022 The MathWorks, Inc.

if size(snapshots,1) ~= size(self.Mesh.Nodes,2)
    error(message('pde:ThermalModel:incorrectSnapshotsSize'));
end

if isvector(snapshots)
    error(message('pde:ThermalModel:insufficientSnapshots'));
end

T = snapshots';
Tmean = mean(T).';
T = T - Tmean.';
numSnapshots = size(T,1);
[~,S,Phi] = svd(T/sqrt(numSnapshots-1));
diagS = diag(S);
percentEnergy = diagS/sum(diagS);
energeticModes = find(cumsum(percentEnergy) < self.SolverOptions.CumulativeEnergyPercentage/100);
Phi = Phi(:,1:nnz(energeticModes));

% Now compute the decay rates using the smaller system obtained by
% projecting the full FE model onto POD modes.
[p,e,t] = self.Mesh.meshToPet();
u0 = pdeuxpd(p,0,self.PDESystemSize);
tdummy = 0:1;
femat = pde.DynamicDiscretizedPDEModel(self,p,e,t,coefstruct,u0,tdummy,isa(self,'pde.StructuralModel'));
K = femat.K;
A = femat.A;
Q = femat.Q;
H = femat.H;
M = femat.Mass;

% Impose Dirichlet BC.
[B,~]=pdenullorth(H);
K = B'*(K+Q+A)*B;
M = B'*M*B;
PhiB = B'*Phi;

Kr = PhiB'*K*PhiB;
Mr = PhiB'*M*PhiB;

% If the asymmetry in the modal system is on the order of 10 times the
% numerical round-off error, then make the matrices symmetric.
Kr = self.checkForSymmetry(Kr);
Mr = self.checkForSymmetry(Mr);

[V,decayrates] = eig(Kr,Mr);
decayrates = diag(decayrates); 

% Make POD modes mass-orthogonal.
R = chol(V.'*Mr*V);
V = V/R;
Phi = Phi*V;

end
