assume(%epsilon[0]>0, %mu[0]>0);

perm(L) := genmatrix(lambda([i,j], if i=L[j] then 1 else 0 ), length(L), length(L));

eigenvalueList(results) := apply(append,makelist(makelist(results[1][1][j],i,1,results[1][2][j]),j,1,length(results[1][1])));

rightEigenvectorMatrix(results) := transpose(apply(matrix, apply(append, results[2])));

load("diag");

/* constant susceptibility */
epsilon_mat : ident(3) * epsilon_0;
mu_mat : ident(3) * mu_0;
/* matrix susceptibility
epsilon_mat : matrix(
    [epsilon_1, -%i * epsilon_2, 0],
    [%i * epsilon_2, epsilon_1, 0],
    [0, 0, %epsilon_z]    
);
mu_mat : ident(3) * mu_0;
*/;

determinant(epsilon_mat);

determinant(mu_mat);

epsilon_inv : ratsimp(invert(epsilon_mat));

mu_inv : ratsimp(invert(mu_mat));

eps_mu_inv_xform : genmatrix(lambda([i,j], 
        if i<=3 and j<=3 then
            epsilon_inv[i,j]
        else if i >= 4 and j >= 4 then
            mu_inv[i-3,j-3]
        else 0
    ), 6, 6);

/* using eps0 E as a conserved quantity... */
[matrix(
[0,0,0,0,0,0],
[0,0,0,0,0,1],
[0,0,0,0,-1,0],
[0,0,0,0,0,0],
[0,0,-1,0,0,0],
[0,1,0,0,0,0]
),
matrix(
[0,0,0,0,0,-1],
[0,0,0,0,0,0],
[0,0,0,1,0,0],
[0,0,1,0,0,0],
[0,0,0,0,0,0],
[-1,0,0,0,0,0]
),
matrix(
[0,0,0,0,1,0],
[0,0,0,-1,0,0],
[0,0,0,0,0,0],
[0,-1,0,0,0,0],
[1,0,0,0,0,0],
[0,0,0,0,0,0]
)];
map(lambda([x], x . eps_mu_inv_xform), %)$
F : %;

results : map(eigenvectors,F);

eigenvalues : map(eigenvalueList, results);

P : perm([1,2,5,6,3,4]) . diag([sqrt(1/2),sqrt(1/2),1,1,sqrt(1/2),sqrt(1/2)]);

eigenvalueMatrix : map(lambda([x], invert(P) . diag(x) . P), eigenvalues);

rightEigenvectors : map(lambda([x], rightEigenvectorMatrix(x) . P), results);

leftEigenvectors : map(lambda([x], ratsimp(invert(x))), rightEigenvectors);

/* verify reconstruction */
makelist([ is(rightEigenvectors[i] . eigenvalueMatrix[i] . leftEigenvectors[i] = F[i]) ],i,1,3);

makelist([ is(rightEigenvectors[i] . leftEigenvectors[i] = ident(6))], i, 1, 3);

makelist(makelist(
    is(F[i] . transpose(transpose(rightEigenvectors[i])[j]) = eigenvalueMatrix[i][j,j] * transpose(transpose(rightEigenvectors[i])[j])),
    j,1,6), i,1,3);

vs : transpose(makelist(v[i-1], i, 1, 6));

makelist( leftEigenvectors[i] . vs, i, 1, 3);

makelist( rightEigenvectors[i] . vs, i, 1, 3);

makelist (F[i] . vs, i, 1, 3);