/*
My other Maxwell eigenvector is based on Trangenstein's book, section 4.3.2.
Its U is [epsilon_0 * E, B]
Subsequently the eigenvectors look a bit more symmetric.
*/;

load("diag")$
load("vect")$

/*
This one is based on 2014 Abgrall, Kumar.
Its U is [B, E].
Also included are potential terms psi and phi for B and E respectively.
*/
/* fluxMethod : "Abgrall Kumar 2014" */

/*
This one is based on Munz 2004
*/
/*fluxMethod : "Munz 2004"*/

/*
This is based on the Munz 2004, except using D instead of E, so the flux looks similar to Trangenstein's, and works in arbitrary medium. 
*/
fluxMethod : "D and B"

;

assume(%gamma > 0, chi > 0, c > 0);

/* returns the matrix skew(A) . B == cross(A, B) */
skew(A) := matrix(
    [0, -A[3], A[2]],
    [A[3], 0, -A[1]],
    [-A[2], A[1], 0]
)$

D : [D_x, D_y, D_z];

H : [H_x, H_y, H_z];

E : [E_x, E_y, E_z];

B : [B_x, B_y, B_z];

I : ident(3);

if fluxMethod = "D and B" then
    U : transpose(matrix(append(D, B, [phi, psi])))
else
    U : transpose(matrix(append(E, B, [phi, psi])));

F : if fluxMethod = "Abgrall Kumar 2014" then
        /* here is the 2014 Abgrall, Kumar flux for perfectly hyperbolic Maxwell equations */ 
        makelist(
            diag([c^2, c^2, c^2, 1, 1, 1, chi, c^2 * %gamma]) . (
                transpose(matrix(append( express( -I[i] ~ B ), express(I[i] ~ E ), [0,0] ))) 
                + diag([chi * phi,  chi * phi,  chi * phi, psi * %gamma, psi * %gamma, psi * %gamma, 0, 0]) . U
                + transpose(matrix([0,0,0,0,0,0, E[i], B[i]]))
            )
            , i, 1, 3)
    else if fluxMethod = "Munz 2004" then
        /* here is the 2000 Munz flux for perfectly hyperbolic Maxwell equations */
        makelist(
            diag([c^2, c^2, c^2, 1, 1, 1, chi, c^2 * %gamma]) . (
                transpose(matrix(append( -express(I[i] ~ B ), express( I[i] ~ E ), [0,0] ))) 
                + transpose(matrix( append( I[i] * phi * chi, I[i] * psi * %gamma, [E[i], B[i]]) ))
            )
            , i, 1, 3)
    else if fluxMethod = "D and B" then
        makelist(
            diag([1, 1, 1, 1, 1, 1, chi, %gamma]) . (
                transpose(matrix(append( express( -I[i] ~ B ) / mu, express( I[i] ~ D ) / epsilon, [0,0] ))) 
                + transpose(matrix( append( I[i] * phi * chi, I[i] * psi * %gamma, [D[i], B[i]]) ))
            )
            , i, 1, 3)

;

A : makelist(
        jacobian(transpose(F[i])[1], transpose(U)[1])
    , i, 1, 3);

results : makelist(eigenvectors(A[i]), i, 1, 3);

perm(L) := apply(matrix, 
    makelist(
        makelist(
            if j = L[i] then 1 else 0
        , j, 1, length(L))
    , i, 1, length(L))
);

P : if fluxMethod = "Abgrall Kumar 2014" then
    /* 2014 Abgrall, Kumar */
    [
        perm([3,6,1,2,7,8,4,5]),
        perm([3,6,1,2,7,8,4,5]),
        perm([3,6,1,2,7,8,4,5])
    ]
    else if fluxMethod = "Munz 2004" then
    /* 2000 Munz */
    [
        transpose(perm([1, 3, 5, 6, 7, 8, 4, 2])),
        transpose(perm([1, 3, 5, 6, 7, 8, 4, 2])),
        transpose(perm([1, 3, 5, 6, 7, 8, 4, 2]))
    ]
    else if fluxMethod = "D and B" then
    [
        transpose(perm([1,3,5,6,7,8,2,4])),
        transpose(perm([1,3,5,6,7,8,2,4])),
        transpose(perm([1,3,5,6,7,8,2,4]))
    ]
;

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

/* maxima leaves off the nullspace, and there's no guarantee on how the remaining nullspace eigenvectors match up with the eigenvalues ... */
evr : if fluxMethod = "Abgrall Kumar 2004" then
    /* 2014 Abgrall, Kumar */
    makelist(
        transpose(apply(matrix, 
            append(
                apply(append, 
                    subst([rhs(lambda_A_def)=lhs(lambda_A_def)],results[k][2])
                ),
                [
                    [
                        [1,0,0,0,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,0,0,1,0,0],
                        [0,0,0,0,0,0,1,0]
                    ]
                ][k]
            )
        )) . P, k, 1, 3)
    else if fluxMethod = "Munz 2004" then
        /* 2000 Munz */
        makelist(
            transpose(apply(matrix, 
                apply(append, results[k][2])
            )) . P[k], k, 1, 3)
    else if fluxMethod = "D and B" then
        [
            transpose(matrix(results[1][2][1][1], [0,1,0,0,0,0,0,0], results[1][2][2][1], [0,0,0,0,0,1,0,0], results[1][2][3][1], results[1][2][4][1], results[1][2][5][1], results[1][2][6][1] )) . P[1],
            transpose(matrix(results[2][2][1][1], [1,0,0,0,0,0,0,0], results[2][2][2][1], [0,0,0,0,0,1,0,0], results[2][2][3][1], results[2][2][4][1], results[2][2][5][1], results[2][2][6][1] )) . P[2],
            transpose(matrix(results[3][2][3][1], results[3][2][4][1], results[3][2][5][1], results[3][2][6][1], [1,0,0,0,0,0,0,0], [0,1,0,0,0,0,0,0], [0,0,0,1,0,0,0,0], [0,0,0,0,1,0,0,0] )) . P[3]
        ]
;

makelist(determinant(evr[k]), k, 1, 3);

evl : makelist(ratsimp(invert(evr[k])), k, 1, 3);

evl : if fluxMethod = "Abgrall Kumar 2014" then
     subst([
        rhs(lambda_B_def)=lhs(lambda_B_def)
    ], evl)
    else 
        evl
;

Areconstruct : makelist(ratsimp(evr[k] . evs[k] . evl[k]), k, 1, 3);

is(Areconstruct = A);

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

makelist(
    evl[k] . vs
    , k, 1, 3);

makelist(
    evr[k] . vs
    , k, 1, 3);