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

load("diag");

W : transpose(matrix([%rho, vx, vy, vz, P]));

U : transpose(matrix([%rho, %rho * vx, %rho * vy, %rho * vz, 1/2 * %rho * (vx^2 + vy^2 + vz^2) + %rho * (P / (%rho * (%gamma - 1)))]));

dU_dW : jacobian(transpose(U)[1], transpose(W)[1]);

dW_dU : ratsimp(invert(dU_dW));

/* flat space */;

A : matrix([0, %rho, 0, 0, 0], [0, 0, 0, 0, 1/%rho], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, %gamma * P, 0, 0, 0]);

results : eigenvectors(A);

A_Lambda : diag([results[1][1][1], results[1][1][3], results[1][1][3], results[1][1][3], results[1][1][2]]);

A_evR : transpose(matrix(
        results[2][1][1],
        results[2][3][1],
        results[2][3][2],
        results[2][3][3],
        results[2][2][1]
    ));

invert(A_evR)$ ratsimp(%)$ A_evL : %;

A_evR . A_Lambda . A_evL$ ratsimp(%);

F_evR : dU_dW . A_evR;

A_evL . dW_dU$ ratsimp(%)$ F_evL : %;

F_Lambda : diag([vx, vx, vx, vx, vx]) + A_Lambda;

F_evR . F_Lambda . F_evL$ ratsimp(%)$ 
subst(vz^2=v^2-vx^2-vy^2,%)$ ratsimp(%)$
%;

/* 



curved space 



*/;

side : 1;

U : transpose(matrix([%rho, %rho * vx, %rho * vy, %rho * vz, 1/2 * %rho * (vx * vxL + vy * vyL + vz * vzL) + %rho * (P / (%rho * (%gamma - 1)))]));

jacobian(transpose(U)[1], transpose(W)[1])$
dU_dW : %$
/* I'm struggling with subst, so I'll just do that manually ... */
dU_dW[5,1] : vSq/2$
/* the result doesn't incorporate dvL/dv , so I'll manually replace that here: */
dU_dW[5,2] : %rho * vxL$ dU_dW[5,3] : %rho * vyL$ dU_dW[5,4] : %rho * vzL$
dU_dW;

invert(dU_dW)$
ratsimp(%)$
subst((2 * %gamma - 2) * vy * vyL = (2 * %gamma - 2) * (vSq - vx * vxL - vz * vzL), %)$
ratsimp(%)$
dW_dU : %;
/* once again, manual simplification ... */
/* dW_dU : %$ dW_dU[5,1] : ratsimp((%gamma-1)*vSq + (1-%gamma)/2*vSq)$ dW_dU; */;

A : [
    /* x direction */
    matrix([0, %rho, 0, 0, 0], [0, 0, 0, 0, gUxx/%rho], [0, 0, 0, 0, gUxy/%rho], [0, 0, 0, 0, gUxz/%rho], [0, %gamma * P, 0, 0, 0]),
    /* y direction */
    matrix([0, 0, %rho, 0, 0], [0, 0, 0, 0, gUxy/%rho], [0, 0, 0, 0, gUyy/%rho], [0, 0, 0, 0, gUyz/%rho], [0, 0, %gamma * P, 0, 0]),
    /* z direction */
    matrix([0, 0, 0, %rho, 0], [0, 0, 0, 0, gUxz/%rho], [0, 0, 0, 0, gUyz/%rho], [0, 0, 0, 0, gUzz/%rho], [0, 0, 0, %gamma * P, 0])
][side];

eigenvectors(A)$
subst(P = %rho * c_s^2 / %gamma, %)$
results : %;

A_Lambda : diag([results[1][1][1], results[1][1][3], results[1][1][3], results[1][1][3], results[1][1][2]]);

A_evr_order : [ [1,2,3], [2,1,3], [2,3,1] ][side]$
A_evR : transpose(matrix(
        results[2][1][1],
        results[2][3][A_evr_order[1]] / %rho * (if side = 1 then %rho else 1),
        results[2][3][A_evr_order[2]] / %rho * (if side = 2 then %rho else 1),
        results[2][3][A_evr_order[3]] / %rho * (if side = 3 then %rho else 1),
        results[2][2][1]
    ));

invert(A_evR)$ ratsimp(%)$ A_evL : %;

A_evR . A_Lambda . A_evL$ subst(c_s = sqrt(%gamma * P / %rho), %)$ ratsimp(%)$ is(% = A);;

dU_dW . A_evR$ 

/* v^x = g^xj v_j */
subst( (c_s * gUxz * vzL) / sqrt(gUxx) = vxpart1, %)$
subst( (c_s * gUxy * vyL) / sqrt(gUxx) = vxpart2, %)$
subst( c_s * sqrt(gUxx) * vxL = vxpart3, %)$
subst( vxpart3 = c_s / sqrt(gUxx) * vx - vxpart1 - vxpart2, %)$

/* v^y = g^yj v_j */
subst( (c_s * gUyz * vzL) / sqrt(gUyy) = vypart1, %)$
subst( (c_s * gUxy * vxL) / sqrt(gUyy) = vypart2, %)$
subst( c_s * sqrt(gUyy) * vyL = vypart3, %)$
subst( vypart3 = c_s / sqrt(gUyy) * vy - vypart1 - vypart2, %)$

/* v^z = g^zj v_j */
subst( (c_s * gUxz * vxL) / sqrt(gUzz) = vzpart1, %)$
subst( (c_s * gUyz * vyL) / sqrt(gUzz) = vzpart2, %)$
subst( c_s * sqrt(gUzz) * vzL = vzpart3, %)$
subst( vzpart3 = c_s / sqrt(gUzz) * vz - vzpart1 - vzpart2, %)$

subst( c_s^2 / (%gamma - 1) = h_total - vSq / 2, % )$
F_evR : %;

A_evL . dW_dU$
F_evL : %;

F_Lambda : diag([vx, vx, vx, vx, vx]) + A_Lambda;

F_evR . F_Lambda . F_evL$
subst(c_s = sqrt(%gamma * P / %rho), %)$
ratsimp(%)$
%;