source: issm/trunk-jpl/src/m/plot/plot_tensor_principalaxis.m

function plot_tensor_principalaxis(md,options,width,i,tensor,type,plot_options)
%PLOT_TENSOR_PRINCIPALAXIS - plot ytensor principal axis
%
%   Usage:
%      plot_tensor_principalaxis(md,options,width,i);
%
%   See also: PLOTMODEL

%prepare subplot
subplot(width,width,i); 

%process data and model
[x y z elements is2d isplanet]=processmesh(md,[],options);

if dimension(md.mesh)==2,
        eval(['Vx=tensor.principalaxis' type(end) '(:,1); Vy=tensor.principalaxis' type(end) '(:,2);'])
        eval(['value=tensor.principalvalue' type(end) ';']);
        [Vx datatype]=processdata(md,Vx,options);
        [Vy datatype]=processdata(md,Vy,options);
        [value datatype]=processdata(md,value,options);
else
        eval(['Vx=tensor.principalaxis' type(end) '(:,1); Vy=tensor.principalaxis' type(end) '(:,2); Vz=tensor.principalaxis' type(end) '(:,3);'])
        [Vx datatype]=processdata(md,Vx,options);
        [Vy datatype]=processdata(md,Vy,options);
        [Vz datatype]=processdata(md,Vz,options);
        [value datatype]=processdata(md,value,options);
end

%take the center of each element if ~isonnode
if datatype==1,
        x=mean(x(elements'))'; y=mean(y(elements'))'; z=mean(z(elements'))';
end

%plot quivers
if dimension(md.mesh)==2,

        %density
        if exist(options,'density')
                density=getfieldvalue(options,'density');
                x=x(1:density:end);
                y=y(1:density:end);
                Vx=Vx(1:density:end);
                Vy=Vy(1:density:end);
                value=value(1:density:end);
        end

        %scaling:
        delta=((min(x)-max(x))^2+(min(y)-max(y))^2)/numel(x);
        scale=0.5/max(sqrt((Vx.^2+Vy.^2)/delta));
        Vx=scale*Vx; Vy=scale*Vy;

        pos=find(value>=0);
        q1=quiver(x(pos),y(pos),Vx(pos),Vy(pos),'Color','r','ShowArrowHead','off','AutoScale','off');
        hold on
        pos=find(value<0);
        q2=quiver(x(pos),y(pos),Vx(pos),Vy(pos),'Color','b','ShowArrowHead','off','AutoScale','off');

else
        %density
        if exist(options,'density')
                density=getfieldvalue(options,'density');
                x=x(1:density:end);
                y=y(1:density:end);
                z=z(1:density:end);
                Vx=Vx(1:density:end);
                Vy=Vy(1:density:end);
                Vz=Vz(1:density:end);
                value=value(1:density:end);
        end

        %scaling:
        delta=((min(x)-max(x))^2+(min(y)-max(y))^2)/numel(x);
        scale=0.5/max(sqrt((Vx.^2+Vy.^2)/delta));
        Vx=scale*Vx; Vy=scale*Vy; Vz=scale*Vz;

        pos=find(value>=0);
        q1=quiver3(x(pos),y(pos),z(pos),Vx(pos),Vy(pos),Vz(pos),'Color','r','ShowArrowHead','off','AutoScale','off');
        hold on
        pos=find(value<0);
        q2=quiver3(x(pos),y(pos),z(pos),Vx(pos),Vy(pos),Vz(pos),'Color','b','ShowArrowHead','off','AutoScale','off');
end

%legend
if strcmpi(type(1:6),'strain')
        legend([q1 q2],'extension','compression')
elseif strcmpi(type(1:6),'stress')
        legend([q1 q2],'compression','traction')
end

%apply options
strings=strsplit_strict(type,'_');
string=strings{1};
options=addfielddefault(options,'title',[upper(string(1)) string(2:end) ' principal axis ' type(end)]);
options=addfielddefault(options,'colorbar',0);
applyoptions(md,[],options);