assume(%gamma>0,P>0,%rho>0);

vSq_def : v^2=v[x]^2+v[y]^2+v[z]^2;

[%rho, v[x], v[y], v[z], P]$
W : %;

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

[%rho * v[x], %rho * v[x]^2 + P, %rho * v[x] * v[y], %rho * v[x] * v[z], (%rho * v^2/2 + %gamma/(%gamma-1) * P) * v[x]]$
subst([vSq_def],%)$
Fx : %;

jacobian(U,W)$
ratsimp(%)$
dU_dW:%;

invert(dU_dW)$
ratsimp(%)$
dW_dU : %;

jacobian(Fx,W)$
ratsimp(%)$
dFx_dW:%;

dW_dU . dFx_dW$
ratsimp(%)$
dFx_dW_quasi : %;

/*acoustic tensor:*/
A : dFx_dW_quasi - ident(5)*v[x];

results : eigenvectors(A);

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

diag(l) := genmatrix(lambda([i,j],if i=j then l[i] else 0),length(l), length(l));

eigenvalues : eigenvalueList(results);

eigenvalueMatrix : diag(eigenvalues);

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

rightEigenvectors : rightEigenvectorMatrix(results);

invert(rightEigenvectors)$
ratsimp(%)$
leftEigenvectors : %;

/* verify reconstruction of A = dFx/dW - I*v[x] */
rightEigenvectors . eigenvalueMatrix . leftEigenvectors$
ratsimp(%);
is(%=A);

vzSq_from_vSq : v[z]^2=v^2-v[x]^2-v[y]^2;

/* verify reconstruction of dFx/dU */
dFx_dW_quasi . dW_dU$
ratsimp(%)$
subst([vzSq_from_vSq],%)$
ratsimp(%)$
%;
/* TODO verify this equals dFx/dU */""$

/* entropy */
S_def : S = S[0] + c[v] * log(P/P[0]*(%rho[0]/%rho)^%gamma);

S0_from_S : solve(S_def,S[0])[1];

/* gradient of the entropy wrt primtives */
'jacobian([S],W)$
subst([S_def],%)$
%,nouns$
dS_dW : %;

/* entropy flux function */
Fs_def : F[S] = S * %rho;

/* gradient of entropy flux function wrt primitives */
'jacobian([F[S]],W)$
subst([Fs_def,S_def],%)$
%,nouns$
subst([S0_from_S],%)$
dFs_dW : %;

/* gradient of the entropy flux function wrt primitives, times the quasilinear flux */
dFs_dW . dFx_dW_quasi$
ratsimp(%)$
dFs_dW_times_dFx_dW_quasi : %;

/* show that the gradient of the entropy flux function wrt the primitives, times the quasilinear flux, equals the gradient of the flux function times the velocity */
'jacobian([F[S]*v[x]],W)$
subst([Fs_def,S_def],%)$
%,nouns$
subst([S0_from_S],%)$
ratsimp(%)$
%;
is(%=dFs_dW_times_dFx_dW_quasi);

/* second derivatives of entropy flux function wrt primitives */
'jacobian(dFs_dW[1], W)$
subst([Fs_def,S_def],%)$
%,nouns$
d2Fs_dW2 : %;

eigenvectors(%);