Vizualizing Orientation Tensors with Matlab/Octave

Recently, I had to vizualize the standard fiber orientation tensor A=<pp> where p is the orientation of a single fiber as shown in the figure below, and the angular brackets denote an average over the distribution function. Although, I explicitly refer to the fiber orientation tensor in this blog, the discussion below is valid for any symmetric tensor, such as a stress tensor, for example.

For some simple cases, the tensor has a direct physical meaning:

For non-simple cases, it can be tricky to mentally visualize the tensor. However, if we calculate the eigenvalues and eigenvectors of the matrix A, visualizing the distribution can be done easily by drawing an ellipsoid. The ellipsoid points in a direction determined by the eigenvectors and the shape and size of the ellipsoid is determined by the eigenvalues.

I modified some code (VizOrient.m) that I found online. It works fine with any relatively recent version of Matlab and Octave 3.0 and higher. Below are five examples:

A=[0.45 0 0; 0 0.45 0; 0 0 0.1];
VizOrient(A);

A=[0.90 0 0; 0 0.05 0; 0 0 0.05];
VizOrient(A);

A=[0.33 0 0; 0 0.33 0; 0 0 0.33];
VizOrient(A);

A=[0.33 0.1 0; 0.1 0.33 0; 0 0 0.33];
VizOrient(A);

A=[0.45 0.2 0; 0.2 0.45 0; 0 0 0.1];
VizOrient(A);

The only problem with Octave is actually a problem with gnuplot. The statement "axis equal" is not strictly implemented, and hence some mental gymnastics could be required. Works great on Matlab, though.
 
References:
1. Chung and Kwon, Korea-Australia Rheology Journal, 14(4), 2002, 175-188

Code:
function VizOrient(A)
%
% take a 3x3 fiber orientation (symmetric) matrix A, and visualize the 
% resulting ellipsoid. using a program from Nima Moshtagh, and erasing the 2D
% lines - also correcting the meaning of a, b, and c axes (1/a)
%
%------------------------------------------------------------------------------
% Copyright (c) 2009, Nima Moshtagh
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
%
%    * Redistributions of source code must retain the above copyright 
%      notice, this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright 
%      notice, this list of conditions and the following disclaimer in 
%      the documentation and/or other materials provided with the distribution
%      
%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
%AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
%IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
%ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
%LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
%CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
%SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
%INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
%CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
%ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
%POSSIBILITY OF SUCH DAMAGE.
%------------------------------------------------------------------------------
N       = 20;
[U D V] = svd(A);
 
%----------------------------------
% generate the ellipsoid at (0,0,0)
%----------------------------------
a = D(1,1);
b = D(2,2);
c = D(3,3);
 
[X,Y,Z] = ellipsoid(0,0,0,a,b,c,N);
    
%-------------------------
% rotate and the ellipsoid
%-------------------------
 
XX = zeros(N+1,N+1);
YY = zeros(N+1,N+1);
ZZ = zeros(N+1,N+1);
 
for k = 1:length(X),
  for j = 1:length(X),
    point = [X(k,j) Y(k,j) Z(k,j)]';
    P = V * point;
    XX(k,j) = P(1);
    YY(k,j) = P(2);
    ZZ(k,j) = P(3);
  end
end
 
%-----------------
% Plot the ellipse
%-----------------
surf(XX,YY,ZZ);
xlabel('x');
ylabel('y');
zlabel('z');
axis equal
hidden on
 
end