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

load("eigen");

/* full metric */;

F : matrix(
[0,1,0,0,0],
[-vx*vx + gxx*(%gamma-1)*vSq/2, 2*vx - (%gamma-1)*gxx*vLx, -(%gamma-1)*gxy*vLy, -(%gamma-1)*gxz*vLz, gxx*(%gamma-1)],
[-vx*vy + gxy*(%gamma-1)*vSq/2, vx - (%gamma-1)*gxy*vLx, vx - (%gamma-1)*gyy*vLy, -(%gamma-1)*gyz*vLz, gxy*(%gamma-1)],
[-vx*vz + gxz*(%gamma-1)*vSq/2, vx - (%gamma-1)*gxz*vLx, -(%gamma-1)*gyz*vLy, vx - (%gamma-1)*gzz*vLz, gxz*(%gamma-1)],
[vx*((%gamma-1)*vSq/2 - hTotal), -(%gamma-1)*vx*vLx + hTotal, -(%gamma-1)*vx*vLy, -(%gamma-1)*vx*vLz, %gamma*vx]);

results : eigenvectors(F);

/* diagonal-only metric */;

matrix(
[0,1,0,0,0],
[-vx*vx + 1/2*(%gamma-1)*vSq/g[xx], vx + vx - (%gamma-1)*vx, -(%gamma-1)*vy*g[yy]/g[xx], -(%gamma-1)*vz*g[zz]/g[xx], (%gamma-1)/g[xx]],
[-vx*vy, vy, vx, 0, 0],
[-vx*vz, vz, 0, vx, 0],
[vx*((%gamma-1)*vSq/2 - hTotal), -(%gamma-1)*vx*vx/g[xx] + hTotal, -(%gamma-1)*vx*vy/g[yy], -(%gamma-1)*vx*vz/g[zz], %gamma*vx]
)$
subst([hTotal=eTotal + P/%rho],%)$
subst([P=(%gamma-1)*%rho*eInt],%)$
subst([eInt=eTotal-eKin],%)$
subst([eKin=vSq/2],%)$
subst([vSq=vx^2/g[xx]+vy^2/g[yy]+vz^2/g[zz]],%)$
ratsimp(%)$
F : %;

results : eigenvectors(F);

/* no metric */;

matrix(
[0,1,0,0,0],
[-vx*vx + (%gamma-1)*vSq/2, 2*vx - (%gamma-1)*vx, - (%gamma-1)*vy, - (%gamma-1)*vz, (%gamma-1)],
[-vx*vy, vy, vx, 0, 0],
[-vx*vz, vz, 0, vx, 0],
[vx*((%gamma-1)*vSq/2 - hTotal), -(%gamma-1)*vx*vx + hTotal, -(%gamma-1)*vx*vy, -(%gamma-1)*vx*vz, %gamma*vx]
)$
subst([hTotal=eTotal + P/%rho],%)$
subst([P=(%gamma-1)*%rho*eInt],%)$
subst([eInt=eTotal-eKin],%)$
subst([eKin=vSq/2],%)$
subst([vSq=vx^2+vy^2+vz^2],%)$
ratsimp(%)$
F : %;

results : eigenvectors(F);