%% Basic linear algebra operations
%  
%  File: ccs_TP_linalg.m
%  Directory: 4_gyujtemegy/11_CCS/_1_ccs/ccs_2018/TP_2018_09_18
%  Author: Peter Polcz (ppolcz@gmail.com) 
%  
%  Created on 2018. September 18.
%
%% Subspace, image and kernel space
% Given an orthonormal basis of the image space of matrix
% 
% $$A = \pmqty{ -2 & -1 & 7 \\ -3 & -5 & 14 \\ 1 & 2 & -5 }.$$

A = [
   -2    -1     7
   -3    -5    14
    1     2    -5
    ];

Im = orth(A);

figure;
pcz_plotvec(Im, 'LineWidth', 2, 'DisplayName', 'Basis of Im(A)'), hold on;
pcz_plotplane(Im, 'DisplayName', 'Im(A)', 'FaceAlpha', 0.5)
shading interp
grid on
L = legend;
L.Interpreter = 'latex';
L.FontSize = 12;


%%%
% Give on orthonormal basis of the kernel (or null) space of matrix $A$.
Ker = null(A);
pcz_plotvec(Ker, 'LineWidth', 2, 'DisplayName', 'Basis of Ker(A)');

%%%
% Give on orthonormal basis of $\rm{Im}(A)^\perp$, the orthogonal
% complement of the image space of matrix $A$.

ImOrthCo = null(Im');
pcz_plotvec(ImOrthCo, 'LineWidth', 2, 'DisplayName', 'Basis of Im(A)$^\perp$')

%%%
% Give on orthonormal basis of the kernel space of matrix $A^T$.
KerT = null(A');
pcz_plotvec(KerT, 'LineWidth', 2, 'DisplayName', 'Basis of Ker(A$^T$)');

%%%
% Given on orthonormal basis of the image space of matrix $A^T$.
ImT = orth(A');
pcz_plotvec(ImT, 'LineWidth', 2, 'DisplayName', 'Basis of Im(A$^T$)');
pcz_plotplane(ImT, 'DisplayName', 'Im(A$^T$)', 'FaceAlpha', 0.5)
shading interp

%%
% Compute the singular value decomposition of matrix $A$.
[U,S,V] = svd(A)

%% 
% Compute the rank of diagonal matrix $S$ and matrix $A$.
rank(S), rank(A)

%% 
% Compute the intersection of $V = \rm{Im}(A)$ and $W = \rm{Im}(A^T)$.
% Note that $V \cap W = \Big(V^\perp \cup W^\perp\Big)^\perp$

V = orth(A);
W = orth(A');

OrthV = null(V');
OrthW = null(W');

Union_OrthV_OrthW = orth( [ OrthV , OrthW ] );

Intersection_V_W = null(Union_OrthV_OrthW');

figure,
pcz_plotvec(Intersection_V_W, 'LineWidth', 2), hold on;
pcz_plotplane(V, 'FaceAlpha', 0.8);
pcz_plotplane(W, 'FaceAlpha', 0.8);
view([ 160.1 , 29.2 ])
shading interp

%% Partial fractional decomposition
% *Using |residue|.*
% 
% Example 1.

syms s

b_sym = s;
a_sym = s^2 + 3*s + 2;

H1_sym = b_sym/a_sym;

a = sym2poly(a_sym);
b = sym2poly(b_sym);

[r,p,k] = residue(b,a)

H1_pfrac_sym = sum(r ./ (s - p));

%%
pcz_align('H(s) = %s = %s', latex(H1_sym), latex(H1_pfrac_sym))


%%
% Example 2.

syms s

b_sym = s^3 + 3*s + 1;
a_sym = (s - 2)*(s + 3)*(s^2 + 2*s + 2)*(s^2 + 4);

H2_sym = b_sym/a_sym;

a = sym2poly(a_sym);
b = sym2poly(b_sym);

[r,p,k] = residue(b,a)

H2_pfrac_sym = vpa(sum(r ./ (s - p)),2)

%%
pcz_align('H(s) = %s = %s', latex(H2_sym), latex(H2_pfrac_sym))


%%
% *Using |partfrac|.*
% 
% Example 1.

H1_pfrac_sym = partfrac(H1_sym, s);

%%
pcz_align('H(s) = %s = %s', latex(H1_sym), latex(H1_pfrac_sym))

%%
% Example 2.

H2_pfrac_sym = partfrac(H2_sym, s, 'factormode', 'rational');

%%
pcz_align('H(s) = %s = %s', latex(H2_sym), latex(H2_pfrac_sym))

%%
H2_pfrac_sym = partfrac(H2_sym, s, 'factormode', 'real');

%%
pcz_align('H(s) = %s = %s', latex(H2_sym), latex(H2_pfrac_sym))

%%
H2_pfrac_sym = partfrac(H2_sym, s, 'factormode', 'complex');

%%
pcz_align('H(s) = %s = %s', latex(H2_sym), latex(H2_pfrac_sym))

%%
H2_pfrac_sym = partfrac(H2_sym, s, 'factormode', 'full');

%%
pcz_align('H(s) = %s = %s', latex(H2_sym), latex(H2_pfrac_sym))
