assume(%rho>0,P>0,a>=0);

U : [%rho, %rho*v[x], %rho*v[x]^2/2+P/(%gamma-1)];

W : [%rho, v[x], P];

dU_dW : jacobian(U,W);

dW_dU : invert(dU_dW);

F : [%rho*v[x], %rho*v[x]^2+P, (%rho*v[x]^2/2+P*%gamma/(%gamma-1))*v[x]];

dF_dW : jacobian(F,W);

dW_dU . dF_dW$
ratsimp(%)$
dF_dW_ql : %;

/* calculate the eigenvectors of the quasilinear flux jacobian wrt primitives */""$

load(eigen);

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

perm : permutationMatrix([1,3,2]);

eigenvalueListFromResults(results) := results[1][1];

eigenvalueMatrixFromResults(results) := block([
    eigenvalueList : eigenvalueListFromResults(results)],
    genmatrix(lambda([i,j], if i=j then eigenvalueList[i] else 0), length(eigenvalueList), length(eigenvalueList))
);

eigenvectorMatrixFromResults(results) := transpose(apply(matrix, map(lambda([l], l[1]), results[2])));

block([
    results : eigenvectors(dF_dW_ql)
],
    %Lambda : perm . eigenvalueMatrixFromResults(results) . perm,
    evr : eigenvectorMatrixFromResults(results) . perm)$

%Lambda;

evr;

evl : invert(evr);

evr . %Lambda . evl$
ratsimp(%);
% - dF_dW_ql;

P_from_a : P = a^2*%rho/%gamma;

a_from_P : solve(P_from_a,a)[2];

dF_dW_ql$
subst([P_from_a],%)$
%;

%Lambda$
subst([P_from_a],%)$
ratsimp(%)$
%;

evr$
subst([P_from_a],%)$
ratsimp(%)$
%;

evl$
subst([P_from_a],%)$
ratsimp(%)$
%;

/* how about the eigenvectors of the flux derivative wrt conservative, then convert them to the quasilinear form ... will that be the same? */""$

dF_dW . dW_dU$
ratsimp(%)$
dF_dU : %;

block([
    results : eigenvectors(dF_dU),
    x
],
    x : perm . eigenvalueMatrixFromResults(results) . perm,
    x : subst([P_from_a],x),
    x : ratsimp(x),
    %Lambda : x,
    x : eigenvectorMatrixFromResults(results) . perm,
    x : subst([P_from_a],x),
    x : ratsimp(x),
    evr : x
)$

%Lambda;

evr;

invert(evr)$
ratsimp(%)$
evl : %;

/* verify reconstruction matches original matrix */
evr . %Lambda . evl$
ratsimp(%)$
subst([a_from_P],%)$
ratsimp(%)$
% - dF_dU;

/* calc eigenvectors, transformed by dU/dW and dW/dU: */
dW_dU . evr$ ratsimp(%);
evl . dU_dW$ ratsimp(%);