evR:matrix(
    [1,1,0,0,1],
    [vx-nx*Cs, vx, ax, bx, vx+nx*Cs],
    [vy-ny*Cs, vy, ay, by, vy+ny*Cs],
    [vz-nz*Cs, vz, az, bz, vz+nz*Cs],
    [hTotal-v_n * Cs, vSq/2, v_a, v_b, hTotal+v_n*Cs]
);

results : ratsimp(invert(evR));

let(ax*by,nz+ay*bx);
let(ay*bz,nx+az*by);
let(az*bx,ny+ax*bz);
let(ay*nx,bz+ax*ny);
let(az*ny,bx+ay*nz);
let(ax*nz,by+az*nx);
let(bx*ny,az+by*nx);
let(by*nz,ax+bz*ny);
let(bz*nx,ay+bx*nz);
let(nz*vz,v_n-nx*vx-ny*vy);
let(az*vz,v_a-ay*vy-ax*vx);
let(bz*vz,v_b-by*vy-bx*vx);
let(az^2,1-ax^2-ay^2);
let(bz^2,1-bx^2-by^2);
let(v_n^2,vSq-v_a^2-v_b^2);
let(nx*v_n,vx-bx*v_b-ax*v_a);
let(ny*v_n,vy-by*v_b-ay*v_a);
let(nz*v_n,vz-bz*v_b-az*v_a);

results$
scanmap(letsimp,%)$
ratsimp(%)$
subst([hTotal=h+vSq/2,h=%gamma*eInt,eInt=P/((%gamma-1)*%rho),P=Cs^2*%rho/%gamma,vSq=2*eKin,%gamma-1=%gamma_1],%)$
ratsimp(%)$
evL:%;