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

depends(E, W);

/* primitives */
W_from_prim : W=[rho, u, P];

E_from_P : E = P/(%gamma-1)+u^2/2;

P_from_E : solve(E_from_P, P)[1];

H_from_E_P : H = E + P;

subst([E_from_P], H_from_E_P)$
ratsimp(%)$
H_from_P:%;

solve(H_from_P,P)[1]$
ratsimp(%)$
P_from_H:%;

subst([P_from_E], H_from_E_P)$
ratsimp(%)$
H_from_E:%;

solve(H_from_E,E)[1]$
ratsimp(%)$
E_from_H:%;

/* conservative */
U_from_cons : U = [rho, rho * u, E];

/* flux */
/* F : transpose(matrix([rho*(u-xi), rho*u*(u-xi)+P, E*(u-xi)+P*u])); */
F = transpose(matrix([rho*u, rho*u*u+P, (E+P)*u]))$
F_from_prim_E : %;

subst(E_from_P, F_from_prim_E)$
F_from_prim : %;

csSq_from_gamma_P_rho : c[s]^2 = %gamma * P / rho;

solve(csSq_from_gamma_P_rho, P)[1]$
P_from_gamma_cs_rho : %;

/* flux jacobian */
genmatrix(lambda([i,j], diff(subst(F_from_prim,F)[i,1], subst(W_from_prim, W)[j])), 3, 3)$
subst([P_from_gamma_cs_rho],%)$
ratsimp(%)$
dF_dW_from_prim_cs : dF_dW = %;

genmatrix(lambda([i,j], diff(subst([E_from_P],subst(U_from_cons,U)[i]),subst(W_from_prim,W)[j])), 3, 3)$
ratsimp(%)$
dU_dW_from_prim : dU_dW = %;

invert(subst(dU_dW_from_prim,dU_dW))$
ratsimp(%)$
dW_dU_from_prim : dW_dU = %;

P_from_gamma_rho_eInt : P = (%gamma - 1) * rho * e[int];

hInt_from_gamma_eInt : h[int] = %gamma * e[int];

hTotal_from_gamma_eInt : h[total] = h[int] + u^2/2;

subst(dF_dW_from_prim_cs,dF_dW) . subst(dW_dU_from_prim,dW_dU)$
ratsimp(%)$
dF_dU = %;