/*
 * 4006.maxima
 *
 * (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
 */

A: matrix([ a+b*%i,      1,       0,      0 ],
          [      0, a+b*%i,       0,      0 ],
          [      0,      0,  a-b*%i,      1 ],
          [      0,      0,       0, a-b*%i ]);

expand(charpoly(A, x));

S: (1/sqrt(2)) * matrix([ 1, -%i,  0,   0 ],
                        [ 0,   0,  1, -%i ],
                        [ 1,  %i,  0,   0 ],
                        [ 0,   0,  1,  %i ]);

B: expand(invert(S).A.S);


C: subst(2, a, B);
C: subst(3, b, C);
A: subst(2, a, A);
A: subst(3, b, A);

U: matrix([ 1, 0, 1, 0 ],
          [ 0, 1, 1, 2 ],
          [ 0, 0, 1, 0 ],
          [ 0, 0, 0, 1 ]);
V: matrix([ 1, 0, 0, 0 ],
          [ 0, 1, 0, 0 ],
          [ 0, 1, 1, 0 ],
          [ 1, 0, 0, 1 ]);
T: U.V;
invert(T);

D: T.C.invert(T);

p: expand(charpoly(D, x));

factor(p);

lambda: 2+3*%i;

Dlambda: ratsimp(expand(D - lambda * identfor(D)));
rank(Dlambda);
/* D2: expand(Dlambda.Dlambda); */
/* rank(D2); */

load(functs);

/*
E: Dlambda;
E[1]: (rational(1/E[1,1]))*E[1]$
E[2]: E[2] - E[2,1] * E[1]$
E[3]: E[3] - E[3,1] * E[1]$
E[4]: E[4] - E[4,1] * E[1]$
E: ratsimp(E)$

E[2]: (rational(1/E[2,2])) * E[2]$
E[3]: E[3] - E[3,2] * E[2]$
E[4]: E[4] - E[4,2] * E[2]$
E: ratsimp(E)$

E[3]: (rational(1/E[3,3])) * E[3]$
E[4]: E[4] - E[4,3] * E[3]$
E: ratsimp(E)$

E[2]: E[2] - E[2,3] * E[3]$
E[1]: E[1] - E[1,3] * E[3]$
E: ratsimp(E)$

E[1]: E[1] - E[1,2] * E[2]$
E: ratsimp(E)$

E;
*/

b1: matrix([1+%i],[2+2*%i],[%i],[1]);
ratsimp(D.b1 - lambda*b1);

G: Dlambda;
G: addcol(G, b1);
G[1]: (rational(1/G[1,1]))*G[1]$
G[2]: G[2] - G[2,1] * G[1]$
G[3]: G[3] - G[3,1] * G[1]$
G[4]: G[4] - G[4,1] * G[1]$
G: ratsimp(G)$

G[2]: (rational(1/G[2,2])) * G[2]$
G[3]: G[3] - G[3,2] * G[2]$
G[4]: G[4] - G[4,2] * G[2]$
G: ratsimp(G)$

G[3]: (rational(1/G[3,3])) * G[3]$
G[4]: G[4] - G[4,3] * G[3]$
G: ratsimp(G)$

G[2]: G[2] - G[2,3] * G[3]$
G[1]: G[1] - G[1,3] * G[3]$
G: ratsimp(G)$

G[1]: G[1] - G[1,2] * G[2]$
G: ratsimp(G)$

G;

b2: matrix([ G[1,5] ], [ G[2,5] ], [ G[3,5] ], [ G[4,5] ]);

expand(D.b2 - lambda * b2 - b1);

c1: 2 * realpart(b1);
d1: 2 * imagpart(b1);
c2: 2 * realpart(b2);
d2: 2 * imagpart(b2);

D.c1 - 2 * c1 + 3 * d1;
D.d1 - 3 * c1 - 2 * d1;
D.c2 - 2 * c2 + 3 * d2 - c1;
D.d2 - 3 * c2 - 2 * d2 - d1;