assume(%gamma>0,%rho>0,P>0,c[s]>=0,c[ax]>=c[s],c[f]>=c[ax],c[a]>=0,c[star]>0,c[snd]>0);

assume(c[f]>=c[snd],c[snd]>=c[s]);

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

/* the symmetrized flux derivative wrt primitives (dW/dU . dF/dU . dU/dW = dW/dU . dF/dW) according to Athena (for W ordered rho,vx,vy,vz,P,by,bz)*/""$

dF_dW_sym_Athena_P5th_7x7 : matrix(
[v[x], %rho, 0, 0, 0, 0, 0],
[0, v[x], 0, 0, 1/%rho, b[y]/%rho, b[z]/%rho],
[0, 0, v[x], 0, 0, -b[x]/%rho, 0],
[0, 0, 0, v[x], 0, 0, -b[x]/%rho],
[0, %rho*c[snd]^2, 0, 0, v[x], 0, 0],
[0, b[y], -b[x], 0, 0, v[x], 0],
[0, b[z], 0, -b[x], 0, 0, v[x]]);

/* Athena dF/dW, Athena's order of primitives with bx added on the end */""$

dF_dW_sym_Athena_P5th_bxLast_8x8 : genmatrix(lambda([i,j], if i<8 and j<8 then dF_dW_sym_Athena_P5th_7x7[i,j] /*else if i = j then v[x]*/ else 0), 8, 8);

/* Athena dF/dW, Ravi's order (rho,vx,vy,vz,bx,by,bz,P) */""$

permute_P5th_to_P8th : perm([1,2,3,4,8,6,7,5]);

dF_dW_sym_Athena_P8th : transpose(permute_P5th_to_P8th) . dF_dW_sym_Athena_P5th_bxLast_8x8 . permute_P5th_to_P8th;

/* Ravi's flux eigenvalues left and right for the flux derivative wrt primitives */""$

/* the eigenvectors of dF/dW in Ravi's code: */
sbx_leftEigenvectors_dF_dW_Ravi_row34 : sbx;
""$

/* trying to fix the rest ... */
sbx_rightEigenvectors_dF_dW_Ravi_row34 : sbx;
sbx_leftEigenvectors_dF_dW_Ravi_col34 : sbx;
/* in Ravi's code: */
sbx_rightEigenvectors_dF_dW_Ravi_67x34 : sbx;
""$

rightEigenvectors_dF_dW_sym_Ravi:matrix(
[1,0,0,0,%rho*%alpha[f],%rho*%alpha[f],%rho*%alpha[s],%rho*%alpha[s]],
[0,0,0,0,%alpha[f]*c[f],-%alpha[f]*c[f],%alpha[s]*c[s],-%alpha[s]*c[s]],
[0,0,-%beta[z]/sqrt(2),-%beta[z]/sqrt(2),-%beta[y]*%alpha[s]*c[s]*sbx_rightEigenvectors_dF_dW_Ravi_row34,%beta[y]*%alpha[s]*c[s]*sbx_rightEigenvectors_dF_dW_Ravi_row34,%beta[y]*%alpha[f]*c[f]*sbx_rightEigenvectors_dF_dW_Ravi_row34,-%beta[y]*%alpha[f]*c[f]*sbx_rightEigenvectors_dF_dW_Ravi_row34],
[0,0,%beta[y]/sqrt(2),%beta[y]/sqrt(2),-%beta[z]*%alpha[s]*c[s]*sbx_rightEigenvectors_dF_dW_Ravi_row34,%beta[z]*%alpha[s]*c[s]*sbx_rightEigenvectors_dF_dW_Ravi_row34,%beta[z]*%alpha[f]*c[f]*sbx_rightEigenvectors_dF_dW_Ravi_row34,-%beta[z]*%alpha[f]*c[f]*sbx_rightEigenvectors_dF_dW_Ravi_row34],
[0,1,0,0,0,0,0,0],
[0,0,sbx_rightEigenvectors_dF_dW_Ravi_67x34*%beta[z]*sqrt(%rho/2),sbx_rightEigenvectors_dF_dW_Ravi_67x34*-%beta[z]*sqrt(%rho/2),%beta[y]*%alpha[s]*sqrt(%rho)*c[snd],%beta[y]*%alpha[s]*sqrt(%rho)*c[snd],-%beta[y]*%alpha[f]*c[snd]*sqrt(%rho),-%beta[y]*%alpha[f]*c[snd]*sqrt(%rho)],
[0,0,sbx_rightEigenvectors_dF_dW_Ravi_67x34*-%beta[y]*sqrt(%rho/2),sbx_rightEigenvectors_dF_dW_Ravi_67x34*%beta[y]*sqrt(%rho/2),%beta[z]*%alpha[s]*sqrt(%rho)*c[snd],%beta[z]*%alpha[s]*sqrt(%rho)*c[snd],-%beta[z]*%alpha[f]*c[snd]*sqrt(%rho),-%beta[z]*%alpha[f]*c[snd]*sqrt(%rho)],
[0,0,0,0,%rho*%alpha[f]*c[snd]^2,%rho*%alpha[f]*c[snd]^2,%rho*%alpha[s]*c[snd]^2,%rho*%alpha[s]*c[snd]^2]);

leftEigenvectors_dF_dW_sym_Ravi:matrix(
[1,0,0,0,0,0,0,-1/c[snd]^2],
[0,0,0,0,1,0,0,0],
[0,0,-%beta[z]/sqrt(2),%beta[y]/sqrt(2),0,%beta[z]*sbx_leftEigenvectors_dF_dW_Ravi_row34/sqrt(2*%rho),-%beta[y]*sbx_leftEigenvectors_dF_dW_Ravi_row34/sqrt(2*%rho),0],
[0,0,-%beta[z]/sqrt(2),%beta[y]/sqrt(2),0,-%beta[z]*sbx_leftEigenvectors_dF_dW_Ravi_row34/sqrt(2*%rho),%beta[y]*sbx_leftEigenvectors_dF_dW_Ravi_row34/sqrt(2*%rho),0],
[0,c[f]*%alpha[f]/(2*c[snd]^2),-%beta[y]*%alpha[s]*c[s]*sbx_leftEigenvectors_dF_dW_Ravi_col34/(2*c[snd]^2),-%beta[z]*%alpha[s]*c[s]*sbx_leftEigenvectors_dF_dW_Ravi_col34/(2*c[snd]^2),0,%beta[y]*%alpha[s]/(2*sqrt(%rho)*c[snd]),%beta[z]*%alpha[s]/(2*sqrt(%rho)*c[snd]),%alpha[f]/(2*%rho*c[snd]^2)],
[0,-c[f]*%alpha[f]/(2*c[snd]^2),%beta[y]*%alpha[s]*c[s]*sbx_leftEigenvectors_dF_dW_Ravi_col34/(2*c[snd]^2),%beta[z]*%alpha[s]*c[s]*sbx_leftEigenvectors_dF_dW_Ravi_col34/(2*c[snd]^2),0,%beta[y]*%alpha[s]/(2*sqrt(%rho)*c[snd]),%beta[z]*%alpha[s]/(2*sqrt(%rho)*c[snd]),%alpha[f]/(2*%rho*c[snd]^2)],
[0,c[s]*%alpha[s]/(2*c[snd]^2),%beta[y]*%alpha[f]*c[f]*sbx_leftEigenvectors_dF_dW_Ravi_col34/(2*c[snd]^2),%beta[z]*%alpha[f]*c[f]*sbx_leftEigenvectors_dF_dW_Ravi_col34/(2*c[snd]^2),0,-%beta[y]*%alpha[f]/(2*sqrt(%rho)*c[snd]),-%beta[z]*%alpha[f]/(2*sqrt(%rho)*c[snd]),%alpha[s]/(2*%rho*c[snd]^2)],
[0,-c[s]*%alpha[s]/(2*c[snd]^2),-%beta[y]*%alpha[f]*c[f]*sbx_leftEigenvectors_dF_dW_Ravi_col34/(2*c[snd]^2),-%beta[z]*%alpha[f]*c[f]*sbx_leftEigenvectors_dF_dW_Ravi_col34/(2*c[snd]^2),0,-%beta[y]*%alpha[f]/(2*sqrt(%rho)*c[snd]),-%beta[z]*%alpha[f]/(2*sqrt(%rho)*c[snd]),%alpha[s]/(2*%rho*c[snd]^2)]);

eigenvalues_Ravi:genmatrix(lambda([i,j],if i=j then [v[x],v[x],v[x]+c[ax],v[x]-c[ax],v[x]+c[f],v[x]-c[f],v[x]+c[s],v[x]-c[s]][i] else 0),8,8);

csndSq_def : c[snd]^2 = %gamma * P / %rho;

csnd_def : sqrt(csndSq_def);

caSq_def : c[a]^2 = (b[x]^2+b[y]^2+b[z]^2)/%rho;

ca_def : sqrt(caSq_def);

caxSq_def : c[ax]^2 = b[x]^2/%rho;

cax_def : sqrt(caxSq_def);

cstarSq_def : c[star]^2 = (c[a]^2+c[snd]^2)/2;

cstar_def : sqrt(cstarSq_def);

cfs_discr : c[star]^4-c[snd]^2*c[ax]^2;

cfSq_def : c[f]^2 = c[star]^2 + sqrt(cfs_discr);

cf_def : sqrt(cfSq_def);

subst([cstarSq_def^2,cstarSq_def,caSq_def^2,csndSq_def^2,caSq_def,csndSq_def],cf_def)$
ratsimp(%)$
cf_def_wrt_prims : %;

csSq_def : c[s]^2 = c[star]^2 - sqrt(cfs_discr);

cs_def : sqrt(csSq_def);

c[f]^2-c[s]^2;
subst([cfSq_def,csSq_def],%);

subst([cstarSq_def^2,cstarSq_def,caSq_def^2,csndSq_def^2,caSq_def,csndSq_def],cs_def)$
ratsimp(%)$
cs_def_wrt_prims : %;

cfSq_minus_csndSq_times_csndSq_minus_csSq_def : (c[f]^2-c[snd]^2)*(c[snd]^2-c[s]^2);

cfSq_minus_csndSq_times_csndSq_minus_csSq_def$
ratsimp(%)$
subst([cs_def],%)$
ratsimp(%)$
subst([cf_def],%)$
ratsimp(%)$
subst([cstar_def],%)$
ratsimp(%)$
subst([caSq_def,caxSq_def],%)$
ratsimp(%)$
factor(%)$
cfSq_minus_csndSq_times_csndSq_minus_csSq_def=%$
%/(c[snd]^2-c[s]^2)$
sqrt(%)$
sqrt_cfSq_minus_csndSq_def:%;

cf_times_cs_def : c[f]*c[s];

cf_times_cs_def$
%^2$
subst([cf_def,cs_def],%)$
ratsimp(%)$
sqrt(%)$
cf_times_cs_def=%$
%/c[s]$
cf_from_cs_csnd_cax_def : %;

/* reconstruct flux.  it matches the Athena papers flux wrt primitives up to a permutation.
Except one detail: Ravi puts ETotal last, while Athena puts this 5th.
*/""$

rightEigenvectors_dF_dW_sym_Ravi . eigenvalues_Ravi . leftEigenvectors_dF_dW_sym_Ravi$
ratsimp(%)$
subst([sbx^2=1],%)$
ratsimp(%)$
subst([%alpha[f] = sqrt( (c[snd]^2-c[s]^2)/(c[f]^2-c[s]^2) )],%)$
ratsimp(%)$
subst([%alpha[s] = sqrt( (c[f]^2-c[snd]^2)/(c[f]^2-c[s]^2) )],%)$
ratsimp(%)$
subst([%beta[y] = b[y] / sqrt(b[y]^2 + b[z]^2)],%)$
ratsimp(%)$
subst([%beta[z] = b[z] / sqrt(b[y]^2 + b[z]^2)],%)$
ratsimp(%)$
subst([sqrt_cfSq_minus_csndSq_def],%)$
ratsimp(%)$
subst([cf_from_cs_csnd_cax_def],%)$
ratsimp(%)$
subst([cax_def],%)$
ratsimp(%)$
subst([sbx*abs(b[x])=b[x]],%)$
ratsimp(%)$
dF_dW_sym_Ravi:%;

/* Ravi's eigenvectors and values (which use primitives rho,vx,vy,vz,bx,by,bz,P) reconstruct close to the Athena flux with its primitives correctly rearranged)
Now to look at the differences ... */""$

dF_dW_sym_Athena_P8th;

/* differences ... */
dF_dW_sym_Ravi - dF_dW_sym_Athena_P8th;

ratioOfMatrices(A,B) := genmatrix(lambda([i,j],if B[i,j] = 0 then (if A[i,j] = 0 then 1 else undefined) else A[i,j]/B[i,j] ),length(A),length(A[1]));

/* ratios ... */
ratioOfMatrices(dF_dW_sym_Athena_P8th, dF_dW_sym_Ravi);

/* derivative of conservative wrt primitives, in Ravi's order (8 variables, P last, ETotal last) */""$

U_Ravi:transpose(matrix([%rho,%rho*v[x],%rho*v[y],%rho*v[z],b[x],b[y],b[z],%rho*(v[x]^2+v[y]^2+v[z]^2)/2+P/(%gamma-1)+(b[x]^2+b[y]^2+b[z]^2)/2]));

W_Ravi:transpose(matrix([%rho,v[x],v[y],v[z],b[x],b[y],b[z],P]));

dU_dW_Ravi:genmatrix(lambda([i,j], diff(U_Ravi[i][1],W_Ravi[j][1])),8,8);

invert(dU_dW_Ravi)$
ratsimp(%)$
dW_dU_Ravi:%;

/* Athena's 8-variable flux vector, with ETotal 5th (rho,mx,my,mz,ETotal,bx,by,bz), with the E[total] flux terms of H[total] */""$

F_Athena_8Var_ETotal5th_wrt_HTotal : transpose(matrix([
%rho*v[x],
%rho*v[x]^2+P+b^2/2-b[x]^2,
%rho*v[x]*v[y]-b[x]*b[y],
%rho*v[x]*v[z]-b[x]*b[z],
H[total]*v[x]-b_dot_v*b[x],
0,
b[y]*v[x]-b[x]*v[y],
b[z]*v[x]-b[x]*v[z]]));

/* Athena's 8-variable flux vector with ETotal 8th, with E[total] flux in terms of H[total] */""$

F_Athena_8Var_ETotal8th_wrt_HTotal : perm([1,2,3,4,8,5,6,7]).F_Athena_8Var_ETotal5th_wrt_HTotal;

/* Athena's 8-variable flux vector, ETotal 8th, with E[total] flux in terms of primitives.  Matches what I have derived as well. */""$

F_Athena_8Var_ETotal8th_wrt_HTotal$
subst([H[total]=E[total]+P[total]],%)$
ratsimp(%)$
subst([E[total]=%rho*v^2/2+P/(%gamma-1)+b^2/2],%)$
ratsimp(%)$
subst([P[total]=P+b^2/2],%)$
ratsimp(%)$
subst([b_dot_v=b[x]*v[x]+b[y]*v[y]+b[z]*v[z]],%)$
ratsimp(%)$
subst([b^2=b[x]^2+b[y]^2+b[z]^2],%)$
ratsimp(%)$
subst([v^2=v[x]^2+v[y]^2+v[z]^2],%)$
ratsimp(%)$
F_Athena_8Var_ETotal8th : %;

/* dF/dW of Athena's / mine 8-variable flux vector (with ETotal 8th and P 8th).  Athena's matches my derivation of dF/dW.  That's reassuring. */""$

genmatrix(lambda([i,j], diff(F_Athena_8Var_ETotal8th[i][1], W_Ravi[j][1])),8,8)$
ratsimp(%)$
dF_dW_Athena_8Var_P8th_fromFlux : %;

/* dW/dU.dF/dW of Athena's 8-variable flux vector */""$

dW_dU_Ravi . dF_dW_Athena_8Var_P8th_fromFlux$
ratsimp(%)$
subst([v[y]^2=v^2-v[x]^2-v[z]^2],%)$
ratsimp(%)$
subst([b[y]^2=b^2-b[x]^2-b[z]^2],%)$
ratsimp(%)$
subst([P=c[snd]^2*%rho/%gamma],%)$
ratsimp(%)$
dF_dW_sym_Athena_P8th_fromFlux : %;

/* compare to the flux they provided ... the 5th col has differences */
dF_dW_sym_Athena_P8th;

/* derivative of flux wrt conservatives, converted from the flux wrt primitives that Athena provides, converted to Ravi's order (P last).
This is what I want in my simulation for verification that the eigenbasis is reconstructing correctly.
Athena provided dW/dU.dF/dW -> dF/dU */""$

dU_dW_Ravi . dF_dW_sym_Athena_P8th . dW_dU_Ravi$
ratsimp(%)$
subst([v[y]^2=v^2-v[x]^2-v[z]^2],%)$
ratsimp(%)$
subst([b[y]^2=b^2-b[x]^2-b[z]^2],%)$
ratsimp(%)$
subst([c[snd]^2=(%gamma-1)*(h[total]-v^2/2-b^2/%rho)],%)$
ratsimp(%)$
subst([2*b[x]*b[z]*v[z]=2*b[x]*(b_dot_v-b[x]*v[x]-b[y]*v[y])],%)$
ratsimp(%)$
dF_dU_Athena_P8th : %;

/* Athena provided flux ... F -> dF/dW -> dW/dU.dF/dW -> dF/dU
To make sure their flux matches their flux, and that they aren't pulling a fast one on me. */""$

dU_dW_Ravi . dF_dW_sym_Athena_P8th_fromFlux . dW_dU_Ravi$
ratsimp(%)$
subst([v[y]^2=v^2-v[x]^2-v[z]^2],%)$
ratsimp(%)$
subst([b[y]^2=b^2-b[x]^2-b[z]^2],%)$
ratsimp(%)$
subst([c[snd]^2=(%gamma-1)*(h[total]-v^2/2-b^2/%rho)],%)$
ratsimp(%)$
subst([2*b[x]*b[z]*v[z]=2*b[x]*(b_dot_v-b[x]*v[x]-b[y]*v[y])],%)$
ratsimp(%)$
dF_dU_Athena_P8th_fromFlux : %;

/* difference of flux matrix from Athena's flux vector, and Athena's flux matrix */
dF_dU_Athena_P8th_fromFlux - dF_dU_Athena_P8th$
ratsimp(%);

/* derivative of flux wrt conservative, converted from the flux reconstructed from Ravi's eigenvectors & values, in Ravi's order (8 vars, P and ETotal last) */""$

dU_dW_Ravi . dF_dW_sym_Ravi . dW_dU_Ravi$
ratsimp(%)$
subst([v[y]^2=v^2-v[x]^2-v[z]^2],%)$
ratsimp(%)$
dF_dU_Ravi : %;

/* rearranging Ravi's order of eigenvalues to Athena's order of eigenvalues */""$

permute_eigenvalues_Ravi_to_Athena : perm([4,5,7,2,8,1,6,3]);

eigenvalues_Athena : permute_eigenvalues_Ravi_to_Athena . eigenvalues_Ravi . transpose(permute_eigenvalues_Ravi_to_Athena);

/* to convert Ravi's eigenvectors of dF/dW, based on Ravi's order of eigenvalues, to eigenvectors based on Athena's order of eigenvalues
you would use the following eigenvectors:
(these should match the Athena paper ... and they almost do, up to some sbx and 2's vs sqrt(2)'s) */""$

permute_P5th_to_P8th . rightEigenvectors_dF_dW_sym_Ravi . transpose(permute_eigenvalues_Ravi_to_Athena) . transpose(perm([1,2,3,4,8,5,6,7]))$
rightEigenvectors_dF_dW_sym_Ravi_inAthenaOrder : %;

perm([1,2,3,4,8,5,6,7]) . permute_eigenvalues_Ravi_to_Athena . leftEigenvectors_dF_dW_sym_Ravi . transpose(permute_P5th_to_P8th)$
leftEigenvectors_dF_dW_sym_Ravi_inAthenaOrder : %;

/* eigenvectors of dF/dU, converted from the dW/dU.dF/dW, with eigenvalues in Athena's order */""$

dU_dW_Ravi . rightEigenvectors_dF_dW_sym_Ravi . transpose(permute_eigenvalues_Ravi_to_Athena)$
ratsimp(%)$
rightEigenvectors_dF_dU_Ravi_AthenaEigenvalueOrder : %;

permute_eigenvalues_Ravi_to_Athena . leftEigenvectors_dF_dW_sym_Ravi . dW_dU_Ravi$
ratsimp(%)$
leftEigenvectors_dF_dU_Ravi_AthenaEigenvectorOrder : %;

/* Athena's eigenvalues of adiabatic MHD flux wrt primitives: */""$

rightEigenvectors_dF_dW_sym_Athena_inAthenaOrder : matrix(
[%rho*%alpha[f],0,%rho*%alpha[s],1,%rho*%alpha[s],0,%rho*%alpha[f],0],
[-%alpha[f]*c[f],0,-%alpha[s]*c[s],0,%alpha[s]*c[s],0,%alpha[f]*c[f],0],
[%alpha[s]*c[s]*sbx*%beta[y],-%beta[z],-%alpha[f]*c[f]*sbx*%beta[y],0,%alpha[f]*c[f]*sbx*%beta[y],%beta[z],-%alpha[s]*c[s]*sbx*%beta[y],0],
[%alpha[s]*c[s]*sbx*%beta[z],%beta[y],-%alpha[f]*c[f]*sbx*%beta[z],0,%alpha[f]*c[f]*sbx*%beta[z],-%beta[y],-%alpha[s]*c[s]*sbx*%beta[z],0],
[%rho*%alpha[f]*c[snd]^2,0,%rho*%alpha[s]*c[snd]^2,0,%rho*%alpha[s]*c[snd]^2,0,%rho*%alpha[f]*c[snd]^2,0],
[sqrt(%rho)*%alpha[s]*c[snd]*%beta[y],-%beta[z]*sbx*sqrt(%rho),-sqrt(%rho)*%alpha[f]*c[snd]*%beta[y],0,-sqrt(%rho)*%alpha[f]*c[snd]*%beta[y],-%beta[z]*sbx*sqrt(%rho),sqrt(%rho)*%alpha[s]*c[snd]*%beta[y],0],
[sqrt(%rho)*%alpha[s]*c[snd]*%beta[z],%beta[y]*sbx*sqrt(%rho),-sqrt(%rho)*%alpha[f]*c[snd]*%beta[z],0,-sqrt(%rho)*%alpha[f]*c[snd]*%beta[z],%beta[y]*sbx*sqrt(%rho),sqrt(%rho)*%alpha[s]*c[snd]*%beta[z],0],
[0,0,0,0,0,0,0,1]);

leftEigenvectors_dF_dW_sym_Athena_inAthenaOrder : matrix(
[0,-(%alpha[f]*c[f])/(2*c[snd]^2),(%alpha[s]*c[s]*sbx*%beta[y])/(2*c[snd]^2),(%alpha[s]*c[s]*sbx*%beta[z])/(2*c[snd]^2),%alpha[f]/(2*%rho*c[snd]^2),(%alpha[s]*%beta[y])/(2*sqrt(%rho)*c[snd]),(%alpha[s]*%beta[z])/(2*sqrt(%rho)*c[snd]),0],
[0,0,-%beta[z]/2,%beta[y]/2,0,-(sbx*%beta[z])/(2*sqrt(%rho)),(sbx*%beta[y])/(2*sqrt(%rho)),0],
[0,-(%alpha[s]*c[s])/(2*c[snd]^2),-(%alpha[f]*c[f]*sbx*%beta[y])/(2*c[snd]^2),-(%alpha[f]*c[f]*sbx*%beta[z])/(2*c[snd]^2),%alpha[s]/(2*%rho*c[snd]^2),-(%alpha[f]*%beta[y])/(2*sqrt(%rho)*c[snd]),-(%alpha[f]*%beta[z])/(2*sqrt(%rho)*c[snd]),0],
[1,0,0,0,-1/c[snd]^2,0,0,0],
[0,(%alpha[s]*c[s])/(2*c[snd]^2),(%alpha[f]*c[f]*sbx*%beta[y])/(2*c[snd]^2),(%alpha[f]*c[f]*sbx*%beta[z])/(2*c[snd]^2),%alpha[s]/(2*%rho*c[snd]^2),-(%alpha[f]*%beta[y])/(2*sqrt(%rho)*c[snd]),-(%alpha[f]*%beta[z])/(2*sqrt(%rho)*c[snd]),0],
[0,0,%beta[z]/2,-%beta[y]/2,0,-(sbx*%beta[z])/(2*sqrt(%rho)),(sbx*%beta[y])/(2*sqrt(%rho)),0],
[0,(%alpha[f]*c[f])/(2*c[snd]^2),-(%alpha[s]*c[s]*sbx*%beta[y])/(2*c[snd]^2),-(%alpha[s]*c[s]*sbx*%beta[z])/(2*c[snd]^2),%alpha[f]/(2*%rho*c[snd]^2),(%alpha[s]*%beta[y])/(2*sqrt(%rho)*c[snd]),(%alpha[s]*%beta[z])/(2*sqrt(%rho)*c[snd]),0],
[0,0,0,0,0,0,0,1]);

/* difference between Athena's eigenvectors of dW/dU.dF/dW and Ravi's (rearranged to Athena's order) */""$

rightEigenvectors_dF_dW_sym_Athena_inAthenaOrder - rightEigenvectors_dF_dW_sym_Ravi_inAthenaOrder;

leftEigenvectors_dF_dW_sym_Athena_inAthenaOrder - leftEigenvectors_dF_dW_sym_Ravi_inAthenaOrder;

/* ratios */""$

ratioOfMatrices(rightEigenvectors_dF_dW_sym_Athena_inAthenaOrder, rightEigenvectors_dF_dW_sym_Ravi_inAthenaOrder);

ratioOfMatrices(leftEigenvectors_dF_dW_sym_Athena_inAthenaOrder, leftEigenvectors_dF_dW_sym_Ravi_inAthenaOrder);

/* Let's see what kind of flux derivative matrix the Athena eigenvectors recreate ... 
... looks similar to Ravi's except for the Alfven vector:
 sqrt(2)'s twice are replaced with 2's once, some sign changes, and some added sqrt(rho)'s to fix the 6-7x6-7 cells.
The 3-4x6-7 and 6-7x3-4 cells don't match Ravi's (by a factor of sqrt(rho)), and don't simplify (easily) down to the correct answer though ... */""$

/* Athena eigenvectors, rearranged to match P being last */""$

transpose(permute_P5th_to_P8th) . rightEigenvectors_dF_dW_sym_Athena_inAthenaOrder . perm([1,2,3,4,8,5,6,7]);

transpose(perm([1,2,3,4,8,5,6,7])) . leftEigenvectors_dF_dW_sym_Athena_inAthenaOrder . permute_P5th_to_P8th;

permute_P5th_to_P8th . rightEigenvectors_dF_dW_sym_Athena_inAthenaOrder . eigenvalues_Athena . leftEigenvectors_dF_dW_sym_Athena_inAthenaOrder . transpose(permute_P5th_to_P8th)$
ratsimp(%)$
subst([sbx^2=1],%)$
ratsimp(%)$
subst([%alpha[f] = sqrt( (c[snd]^2-c[s]^2)/(c[f]^2-c[s]^2) )],%)$
ratsimp(%)$
subst([%alpha[s] = sqrt( (c[f]^2-c[snd]^2)/(c[f]^2-c[s]^2) )],%)$
ratsimp(%)$
subst([%beta[y] = b[y] / sqrt(b[y]^2 + b[z]^2)],%)$
ratsimp(%)$
subst([%beta[z] = b[z] / sqrt(b[y]^2 + b[z]^2)],%)$
ratsimp(%)$
subst([sqrt_cfSq_minus_csndSq_def],%)$
ratsimp(%)$
subst([cf_from_cs_csnd_cax_def],%)$
ratsimp(%)$
/*subst([cax_def],%)$
ratsimp(%)$
subst([sbx*abs(b[x])=b[x]],%)$
ratsimp(%)$
*/
AwA:%;