load("diag");

assume(detgvar>0, epsilon_0>0, mu_0>0, g_xx>0, g_yy>0, g_zz>0, gUxx>0, gUyy>0, gUzz>0);

g : matrix(
    [g_xx, g_xy, g_xz],
    [g_xy, g_yy, g_yz],
    [g_xz, g_yz, g_zz]
);

detg : determinant(g);

ggU : ratsimp(detg * invert(g));

gUvars : matrix(
    [gUxx, gUxy, gUxz],
    [gUxy, gUyy, gUyz],
    [gUxz, gUyz, gUzz]
);

matrix(
    [0,0,0,zeta/sqrt(detgvar),0,0],
    [0,0,0,0,0,1/(sqrt(detgvar)*mu_0)],
    [0,0,0,0,-1/(sqrt(detgvar)*mu_0),0],
    [zeta/sqrt(detgvar),0,0,0,0,0],
    [0,0,-1/(sqrt(detgvar)*epsilon_0),0,0,0],
    [0,1/(sqrt(detgvar)*epsilon_0),0,0,0,0]
)$
subst([zeta=0],%)$
flatM : %;

m : flatM
 .
matrix(
    [g_xx, g_xy, g_xz,0,0,0],
    [g_xy, g_yy, g_yz,0,0,0],
    [g_xz, g_yz, g_zz,0,0,0],
    [0,0,0,g_xx, g_xy, g_xz],
    [0,0,0,g_xy, g_yy, g_yz],
    [0,0,0,g_xz, g_yz, g_zz]
);
/*
m : matrix(
    [0, 0, 0, 0, 0, 0],
    [0, 0, 0, g_zx, g_zy, g_zz]/mu_0,
    [0, 0, 0, -g_yx, -g_yy, -g_zy]/mu_0,
    [0, 0, 0, 0, 0, 0],
    [-g_zx, -g_zy, -g_zz, 0, 0, 0]/epsilon_0,
    [g_yx, g_yy, g_zy, 0, 0, 0]/epsilon_0
) / sqrt(detgvar);
*/;

results : eigenvectors(m);

results$
subst([g_zx=g_xz, g_zy=g_yz, g_yx=g_xy], %)$
subst(
    apply(append, apply(append, makelist(makelist([ggU[i,j]=gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(2*ggU[i,j]) = 2*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(4*ggU[i,j]) = 4*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(-4*ggU[i,j]) = -4*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(4*g_zz*ggU[i,j]) = 4*g_zz*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
ratsimp(%)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(2*ggU[i,j]) = 2*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ggU[i,j]=gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(-ggU[i,j]) = -gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(-detgvar*ggU[i,j]) = -detgvar*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(epsilon_0*ggU[i,j]) = epsilon_0*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$

ratsimp(%)$
results : %;

P : matrix(
    [0,1,0,0,0,0],
    [1,0,0,0,0,0],
    [0,0,0,0,1,0],
    [0,0,0,0,0,1],
    [0,0,1,0,0,0],
    [0,0,0,1,0,0]
);

diag(apply(append, makelist(makelist(results[1][1][i], k, 1, results[1][2][i]), i, 1, length(results[1][1]))))$
P . % . transpose(P)$
Lambda : %;

apply(append, makelist(results[2][i], i, 1, length(results[2])))$
/* maxima is missing some eigenvectors ... */
append(%, [
        /*
        [sqrt(gUxx/mu_0),0,0,sqrt(detgvar/epsilon_0)*g_xx, g_xz*sqrt(detgvar/epsilon_0), -g_xy*sqrt(detgvar/epsilon_0)],
        [sqrt(gUxx/mu_0),0,0,sqrt(detgvar/epsilon_0)*-g_xx,-g_xz*sqrt(detgvar/epsilon_0), g_xy*sqrt(detgvar/epsilon_0)]
        */
        [-1,0,0,sqrt(detgvar/gUxx)*g_xx,0,0],
        [1,0,0,sqrt(detgvar/gUxx)*g_xx,0,0]
    ])$
transpose(apply(matrix, %))$
% . transpose(P)$
% . diag([
        sqrt(gUxx/(2*mu_0)),
        -sqrt(gUxx/(2*mu_0)),
        sqrt(gUxx/(2*mu_0)),
        sqrt(gUxx/(2*mu_0)),
        sqrt(gUxx/(2*mu_0)),
        -sqrt(gUxx/(2*mu_0))
    ])$
ratsimp(%)$

subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(sqrt(2)*epsilon_0*ggU[i,j]) = sqrt(2)*epsilon_0*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
ratsimp(%)$
% * sqrt(mu_0 * epsilon_0)$
evR : %$

/*
evR[4,1] : sqrt(detgvar/(2*epsilon_0))*g_xz$
evR[4,6] : -sqrt(detgvar/(2*epsilon_0))*g_xz$
evR[4,2] : sqrt(detgvar/(2*epsilon_0))*g_xy$
evR[4,5] : -sqrt(detgvar/(2*epsilon_0))*g_xy$
*/
evR;

evR$
invert(%)$
ratsimp(%)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(sqrt(2)*detgvar*ggU[i,j]) = sqrt(2)*detgvar*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(2*gUxx*ggU[i,j]) = 2*gUxx*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(2*g_xx*ggU[i,j]) = 2*g_xx*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(sqrt(2)*ggU[i,j]) = sqrt(2)*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(-sqrt(2)*ggU[i,j]) = -sqrt(2)*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(epsilon_0*ggU[i,j]) = epsilon_0*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(-epsilon_0*ggU[i,j]) = -epsilon_0*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
evL : %;

/* hmm, this doesn't reconstruct everything ... */
evR . Lambda . evL$
ratsimp(%)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(detgvar*ggU[i,j]) = detgvar*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
ratsimp(%)$
subst(
    apply(append, apply(append, makelist(makelist([ ratsimp(-detgvar*ggU[i,j]) = -detgvar*gUvars[i,j]/detgvar],i,1,3),j,1,3)))
, %)$
ratsimp(%)$

%;

test_A : matrix([0,0],[1,1]);

test_results : eigenvectors(test_A);

test_Lambda : diag([0,1]);

test_evR : transpose(matrix([1,-1],[0,1]));

test_evL : ratsimp(invert(test_evR));

test_evR . test_Lambda . test_evL;