assume(rho > 0, P_star > 0, R > 0, C_v > 0, %gamma_1 > 0, Cs > 0);

load("diag");

swap(eqn) := rhs(eqn) = lhs(eqn);

/*
rho <-> bar rho
v_i <-> tilde v_i
e_total <=> tilde E
P_star <-> P*
P_bar <-> bar P
T <-> tilde T

gamma-1 = gamma_1 = R/C_v
alpha = 1 - 2/3 C_v/R = 1 - 2/3 /(gamma-1)
*/;

gUs : genmatrix(lambda([i,j], if i <= j then gU[i,j] else gU[j,i]),3,3);

gLs : genmatrix(lambda([i,j], if i <= j then gL[i,j] else gL[j,i]),3,3);

vLs : [vLx, vLy, vLz];

vUs : [vUx, vUy, vUz];

vSq_def : vSq = vUs . gLs . vUs;

transpose(vLs) = gLs . vUs$
makelist(lhs(%)[i,1] = rhs(%)[i,1], i, 1, 3)$
vL_for_vU : %;

vL_for_vU$
makelist(swap(%[i]), i, 1, 3)$
vL_from_vU : %;

transpose(vUs) = gUs . vLs$
makelist(lhs(%)[i,1] = rhs(%)[i,1], i, 1, 3)$
vU_for_vL : %;

vU_for_vL$
makelist(swap(%[i]), i, 1, 3)$
vU_from_vL : %;

W : transpose(matrix([rho, vUx, vUy, vUz, P_star, k, omega]));

etotal_from_eint_v_k : e_total = e_int + 1/2 * vSq + k;

eint_def : e_int = C_v * T;

T_from_eint : solve(eint_def, T)[1];

T_def : T = P_bar / (rho * R);

Pbar_from_rho_T : solve(T_def, P_bar)[1];

Pbar_def : P_bar = P_star - 2/3 * rho * k;

Pstar_from_rho_k_Pbar : solve(Pbar_def, P_star)[1];

subst([Pbar_from_rho_T, T_from_eint], Pstar_from_rho_k_Pbar)$
ratsimp(%);

alpha = 1 - (2 * C_v) / (3 * R);
ratsimp(%);
alpha_def : %$

R_for_alpha : solve(alpha_def, R)[1];

gamma1_def : %gamma_1 = R / C_v;

R_for_gamma1 : solve(gamma1_def, R)[1];

subst([R_for_gamma1], R_for_alpha)$
ratsimp(% / C_v)$
gamma1_for_alpha : %;

solve(gamma1_for_alpha, alpha)[1];
alpha_for_gamma1 : %$

subst([eint_def, T_def, Pbar_def, R_for_gamma1], etotal_from_eint_v_k)$
ratsimp(%)$
etotal_from_W : %;

solve(etotal_from_W, P_star)[1]$
Pstar_from_W : %;

transpose(matrix([rho, rho * vUx, rho * vUy, rho * vUz, rho * e_total, rho * k, rho * omega]))$
U : %;

U$
subst([etotal_from_W], %)$
subst([vSq_def], %)$
jacobian( transpose(%)[1], transpose(W)[1] )$
subst([swap(vSq_def)], %)$
subst(  ratsimp(rho*vL_from_vU), ratsimp(%))$
dU_dW : %;

dU_dW$
invert(%)$ ratsimp(%)$
subst([vSq = vLs . vUs], %)$ ratsimp(%)$
subst( ratsimp(%gamma_1 / 2 * [vLs . vUs = vSq]), %)$
dW_dU : %;

makelist(
transpose(matrix([
            vUs[j]*rho,
            vUs[j]*rho*vUx + P_star*gUs[j,1],
            vUs[j]*rho*vUy + P_star*gUs[j,2],
            vUs[j]*rho*vUz + P_star*gUs[j,3],
            vUs[j]*(rho*e_total + P_star),
            vUs[j]*rho*k,
            vUs[j]*rho*omega
        ])),
    j,1,3)$
subst([etotal_from_W], %)$
ratsimp(%)$
F_def : %;

F_def$
subst([vSq_def], %)$
makelist(ratsimp(jacobian(transpose(%[i])[1], transpose(W)[1])), i,1,3)$

makelist( 
    subst(ratsimp(vUs[i] * (1/2 * swap(vSq_def) + (1 - 2/3/%gamma_1) * k) ), %[i]),
    i, 1, 3)$ ratsimp(%)$

makelist(
    subst(
        apply(append, makelist(makelist(ratsimp(vUs[k]*rho*swap(vL_for_vU[j])), j, 1, 3), k, 1, 3)),
        %[i]),
    i, 1, 3)$ ratsimp(%)$

makelist(
    subst(
        [ratsimp(swap(rho*(1/2*vSq_def + (1 - 2/3/%gamma_1)*k + vUs[i] * vL_for_vU[i]) + (1+1/%gamma_1)*P_star))],
        %[i]),
    i,1,3)$ ratsimp(%)$

dF_dW_def : %;

makelist(dF_dW_def[i] . dW_dU, i, 1, 3)$ ratsimp(%)$
makelist( subst([ ratsimp( -((vUs.vLs=vSq) - 1/2 * (%gamma_1 + 1) * vSq + k*(1-2/3/%gamma_1) + (1/%gamma_1+1)*P_star / rho) * vUs[i] ) ], %[i]), i,1,3)$ ratsimp(%)$
dF_dU_def : %;

makelist(ratsimp(dW_dU . dF_dW_def[i] - ident(7) * vUs[i]), i,1,3)$
subst([vSq=vLs.vUs], %)$ ratsimp(%)$
subst(vU_for_vL, %)$ ratsimp(%)$
acoustic_def : %;

results : makelist(eigenvectors(acoustic_def[i]), i, 1,3);

CsSq_def : Cs^2 =  (%gamma_1+1) * P_star / rho;

ratsimp(sqrt(CsSq_def))$
Cs_def : %;

solve(Cs_def, P_star)[1]$
Pstar_for_Cs : %;

perm : [
    matrix(
        [1,0,0,0,0,0,0],
        [0,0,0,0,0,0,1],
        [0,1,0,0,0,0,0],
        [0,0,1,0,0,0,0],
        [0,0,0,1,0,0,0],
        [0,0,0,0,1,0,0],
        [0,0,0,0,0,1,0]
    ),
    matrix(
        [1,0,0,0,0,0,0],
        [0,0,0,0,0,0,1],
        [0,0,1,0,0,0,0],
        [0,1,0,0,0,0,0],
        [0,0,0,1,0,0,0],
        [0,0,0,0,1,0,0],
        [0,0,0,0,0,1,0]
    ),
    matrix(
        [1,0,0,0,0,0,0],
        [0,0,0,0,0,0,1],
        [0,0,0,1,0,0,0],
        [0,1,0,0,0,0,0],
        [0,0,1,0,0,0,0],
        [0,0,0,0,1,0,0],
        [0,0,0,0,0,1,0]
    )]    ;

iperm : makelist( ratsimp(invert(perm[i])), i, 1,3);

makelist( diag(apply(append, makelist( makelist(results[k][1][1][j], i,1, results[k][1][2][j]), j,1,length(results[k][1][1])))), k, 1,3)$
makelist( subst([Pstar_for_Cs], %[i]), i, 1, 3)$ ratsimp(%)$
makelist( iperm[i] . %[i] . perm[i], i, 1, 3)$
acoustic_Lambda_def : %;

makelist( transpose(apply(append, makelist(apply(matrix, results[k][2][i]),i,1,3))), k, 1, 3)$
makelist( subst([Pstar_for_Cs], %[i]), i,1,3)$ ratsimp(%)$
makelist( subst([gamma1_for_alpha], %[i]), i,1,3)$ ratsimp(%)$
makelist(%[i] . perm[i], i, 1, 3)$
acoustic_evR : %;

acoustic_evR$
makelist(invert(%[i]), i,1,3)$ ratsimp(%)$
makelist(subst([alpha_for_gamma1], %[i]), i,1,3)$ ratsimp(%)$
acoustic_evL : %;

/* verify that the left and right eigenvectors and eigenvalues of the acoustic matrix reconstruct it correctly */
makelist(acoustic_evR[i] . acoustic_Lambda_def[i] . acoustic_evL[i] - acoustic_def[i], i,1,3)$ ratsimp(%)$
makelist( subst([Cs_def, alpha_for_gamma1], %[i]), i,1,3)$ ratsimp(%)$
%;

makelist(acoustic_Lambda_def[i] + ident(7) * vUs[i], i,1,3)$ ratsimp(%)$
Lambda_def : %;

perm2 : [
    diag([1,1,1/rho,1/rho,1/rho,1/rho,1]),
    diag([1,1/rho,1,1/rho,1/rho,1/rho,1]),
    diag([1,1/rho,1/rho,1,1/rho,1/rho,1])
];

makelist(dU_dW . acoustic_evR[i], i,1,3)$ ratsimp(%)$
makelist(subst([ratsimp(swap(6 * Cs * %gamma_1 * vU_for_vL[i]))], %[i]), i,1,3)$ ratsimp(%)$
evR : %;

makelist(acoustic_evL[i] . dW_dU, i,1,3)$ ratsimp(%)$
evL : %;

/* verify that the flux jacobian wrt cons.' eigenvectors & values constructed by the acoustic matrix eigenvectors & values matches the flux jacobian wrt cons.' constructed by the flux jacobian wrt prim.'s*/
makelist(evR[i] . Lambda_def[i] . evL[i] - dF_dU_def[i], i,1,3)$ ratsimp(%)$
makelist( subst([CsSq_def, vSq=vUs.vLs], %[i]), i,1,3)$ ratsimp(%)$
makelist( subst(vU_for_vL, %[i]), i,1,3)$ ratsimp(%)$
%;

Xs : transpose(makelist(X[i], i,0,6));

makelist( (dF_dU_def[i] . Xs), i,1,3);