%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Matlab/Simulink program for Results comparison %%%%%%%
%%%%%%% between full ANSYS and reduced MOR simulations %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Pascal Maglie %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% maglie@iwf.mavt.ethz.ch %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% IWF - Institut für Werkzeugmaschinen und Fertigung %%%%%
%%%%% ETHZ - Eidgenössische Technische Hochschule Zürich %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Transformation from Matrix Market Format into Matlab format
[K, rowsK, colsK, entriesK] = MOR2MATLAB('mor.K');
[M, rowsM, colsM, entriesM] = MOR2MATLAB('mor.M');
[B, rowsB, colsB, entriesB] = MOR2MATLAB('mor.B');
[C, rowsC, colsC, entriesC] = MOR2MATLAB('mor.C');


n = 20;     % Order of the reduced model
b = 1;      % Dimension of the Input vector u
c = 2;      % Dimension of the Output vector y

% nxn identity matrix
IDn = eye(n);
% nxn null matrix
ONn = zeros(n,n);
% nxb null matrix
OBn = zeros(n,b);
% cxn null matrix
OCn = zeros(c,n);

% Damping coefficients (alpha & beta)
alpha = 0;
beta = 0.00318;

% Rayleigh damping matrix D
D = alpha*M + beta*K;

% Augmented M matrix
M_aug = [M ONn ; ONn IDn];
% Inverse of the augmentes M matrix
M_inv = inv(M_aug);

% State-Space A Matrix
A_aug = M_inv*[-D -K ; IDn ONn];
% State-Space B Matrix
B_aug = M_inv*[B OBn ; OBn OBn];
% State-Space C Matrix
C_aug = [OCn C];
% State-Space D Matrix
D_aug = [0 0 ; 0 0];