	p, r, u, v;

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

	//==== Jacobian matrix transformation in x-direction:====
	// Ax = dFx / dU, Ax = [Sx-1][Lx][Sx], dUx ~= [Sx-1]*dWx, where dWx ~= [Sx]*dUx
	A[0][0] = 0;                          A[0][1] = 1;                        A[0][2] = 0;          A[0][3] = 0;
	A[1][0] = u*u*(K-3)/2 + v*v*(K-1)/2;  A[1][1] = (3-K)*u;                  A[1][2] = (1-K)*v;    A[1][3] = K-1;
	A[2][0] = -u*v;                       A[2][1] = v;                        A[2][2] = u;          A[2][3] = 0;
	A[3][0] = -K*E*u+(K-1)*u*(u*u+v*v);   A[3][1] = K*E+(1-K)*(3*u*u+v*v)/2;  A[3][2] = (1-K)*u*v;  A[3][3] = K*u;

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

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

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


	//==== Jacobian matrix transformation in Y-direction: ====
	// Ay = dFy / dU, Ay = [Sy-1][Ly][Sy], dUy ~= [Sy-1]*dWy, where dWy ~= [Sy]*dUy
	A[0][0] = 0;                          A[0][1] = 0;         A[0][2] = 1;                       A[0][3] = 0;
	A[1][0] = -u*v;                       A[1][1] = v;         A[1][2] = u;                       A[1][3] = 0;
	A[2][0] = v*v*(K-3)/2 + u*u*(K-1)/2;  A[2][1] = (1-K)*u;   A[2][2] = (3-K)*v;                 A[2][3] = K-1;
	A[3][0] = -K*E*v+(K-1)*v*(u*u+v*v);   A[3][1] = (1-K)*u*v; A[3][2] = K*E+(1-K)*(3*v*v+u*u)/2; A[3][3] = K*v;

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

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

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