	p, r, u;

	HI = K-1.;
	c = sqrt( K * p / r );
	E = p / ( HI * r) + u * u / 2;

	//==== Jacobian matrix transformation in -direction:====
	//  = dF / dU, A = [S-1][L][S], dU ~= [S-1]*dW, where dW ~= [S]*dU
	A[0][0] = 0;                  A[0][1] = 1;                    A[0][2] = 0;
	A[1][0] = u*u*(K-3)/2;        A[1][1] = (3-K)*u;              A[1][2] = K-1;
	A[2][0] = -K*E*u+(K-1)*u*u*u; A[2][1] = K*E+(1-K)*(3*u*u)/2;  A[2][2] = K*u;

	// Diagonal matrix of eigenvalues
	L[0][0] = u;  L[0][1] = 0;    L[0][2] = 0;
	L[1][0] = 0;  L[1][1] = u+c;  L[1][2] = 0;
	L[2][0] = 0;  L[2][1] = 0;    L[2][2] = u-c;

	// Transformation matrix
	S[0][0] = -c*c + (K-1)*u*u/2;  S[0][1] =  -(K-1)*u;  S[0][2] = (K-1);
	S[1][0] = -c*u + (K-1)*u*u/2;  S[1][1] = c-(K-1)*u;  S[1][2] = (K-1);
	S[2][0] =  c*u + (K-1)*u*u/2;  S[2][1] =-c-(K-1)*u;  S[2][2] = (K-1);

	// Inverted transformation matrix
	S_[0][0] = -1/(c*c);       S_[0][1] = 1/(2*c*c);                           S_[0][2] = 1/(2*c*c);
	S_[1][0] = -u/(c*c);       S_[1][1] = u/(2*c*c)+1/(2*c);                   S_[1][2] = u/(2*c*c)-1/(2*c);
	S_[2][0] = -(u*u)/(2*c*c); S_[2][1] = u*u/(4*c*c) + u/(2*c) + 1/(2*(K-1)); S_[2][2] = u*u/(4*c*c) - u/(2*c) + 1/(2*(K-1));

