/**

Script to compute the eigensystem of the 10-moment equation
system. The heat-tensor is set to 0. Load into Maxima using the
batch("tenmom-eig.mac"); command.

Author: Ammar Hakim.
 */

/** Allocate space for matrix */
Atm : zeromatrix(10, 10);

/** Define constants to make it easy for creating matrix */
N : 1;
U1 : 2;
U2 : 3;
U3 : 4;
P11 : 5;
P12 : 6;
P13 : 7;
P22 : 8;
P23 : 9;
P33 : 10;

n : p0; /** it is better to work with rho (p0) than number density */

/** Positivity for system */
assume(n>0);
assume(p0>0);
assume(p11>0);

/** n */
Atm[N,N] : u1; Atm[N,U1] : n;

/** u1 */
Atm[U1,U1] : u1; Atm[U1,P11] : 1/p0;

/** u2 */
Atm[U2,U2] : u1; Atm[U2,P12] : 1/p0;

/** u3 */
Atm[U3,U3] : u1; Atm[U3,P13] : 1/p0;

/** p11 */
Atm[P11,U1] : 3*p11; Atm[P11,P11] : u1;

/** p12 */
Atm[P12,U1] : 2*p12; Atm[P12,U2] : p11; Atm[P12,P12] : u1;

/** p13 */
Atm[P13,U1] : 2*p13; Atm[P13,U3] : p11; Atm[P13,P13] : u1;

/** p22 */
Atm[P22,U1] : p22; Atm[P22,U2] : 2*p12; Atm[P22,P22] : u1;

/** p23 */
Atm[P23,U1] : p23; Atm[P23,U2] : p13; Atm[P23,P23] : u1; Atm[P23,U3] : p12;

/** p33 */
Atm[P33,U1] : p33; Atm[P33,U3] : 2*p13; Atm[P33,P33] : u1;

/** Now compute eigenvalues */
/**eigVals : eigenvalues(Atm);*/

/** Compute eigensystem. */
[vals, vects] : eigenvectors(Atm); 

/** Right eigenvectors. This nastiness is needed as the eigesystem is
returned as a list of nested lists. These nested lists are taken apart
based on eigenvalue multiplicity and put as columns in a matrix of
right-eigenvectors */

Rev : (Rev : matrix([]), for i from 1 thru length (vals[1]) 
        do (for j from 1 thru vals[2][i] 
          do ( (Rev : addcol(Rev, transpose(matrix(vects[i][j])))))), Rev); 

/** Massages right-eigenvectors */
Revm : matrix([]);
Revm : addcol( Revm, -sqrt(p11)/sqrt(p0)*col(Rev,1) );
Revm : addcol( Revm, -sqrt(p11)/sqrt(p0)*col(Rev,2) );
Revm : addcol( Revm, sqrt(p11)/sqrt(p0)*col(Rev,3) );
Revm : addcol( Revm, sqrt(p11)/sqrt(p0)*col(Rev,4) );

Revm : addcol( Revm, p0*p11*col(Rev,5) );
Revm : addcol( Revm, p0*p11*col(Rev,6) );

Revm : addcol( Revm, col(Rev,7) );
Revm : addcol( Revm, col(Rev,8) );
Revm : addcol( Revm, col(Rev,9) );
Revm : addcol( Revm, col(Rev,10) );

/** Compute left-eigenvectors */
Levm : fullratsimp( invert(Revm) ); /** Reduce as best you can */
/** Compute check: this should be an identiy matrix */
id : fullratsimp( Levm . Revm );

/** Now create the Jacobian of the transformation */
phiprime : zeromatrix(10, 10);

/** n */
phiprime[N,N] : 1;

/** u1 */
phiprime[U1,N] : u1; phiprime[U1,U1] : n;

/** u2 */
phiprime[U2,N] : u2; phiprime[U2,U2] : n;

/** u3 */
phiprime[U3,N] : u3; phiprime[U3,U3] : n;

/** p11 */
phiprime[P11,N] : u1*u1; phiprime[P11,U1] : 2*n*u1; phiprime[P11,P11] : 1;

/** p12 */
phiprime[P12,N] : u1*u2; phiprime[P12,U1] : n*u2; phiprime[P12,U2] : n*u1; phiprime[P12,P12] : 1; 

/** p13 */
phiprime[P13,N] : u1*u3; phiprime[P13,U1] : n*u3; phiprime[P13,U3] : n*u1; phiprime[P13,P13] : 1;

/** p22 */
phiprime[P22,N] : u2*u2; phiprime[P22,U2] : 2*n*u2; phiprime[P22,P22] : 1;

/** p23 */
phiprime[P23,N] : u2*u3; phiprime[P23,U2] : n*u3; phiprime[P23,U3] : n*u2; phiprime[P23,P23] : 1;

/** p33 */
phiprime[P33,N] : u3*u3; phiprime[P33,U3] : 2*n*u3; phiprime[P33,P33] : 1;

/** Now compute inverse transform */
phiprimeInv : fullratsimp(invert(phiprime));