%  Generate some random unit vectors
%  and store them in the rows of a matrix n

numpts=2000;
r1 =-1+2*rand(numpts,1);
r2 =-1+2*rand(numpts,1);
r3 =-1+2*rand(numpts,1);
n = [r1 r2 r3];

for ic = 1:numpts
    fac = norm(n(ic,:));
    n(ic,:)=n(ic,:)/fac;
end

% Here is a state of stress

sigma = [ 10 3 0;  3 -5 0; 0 0 2]

%  for each plane given by a row in n, find
%  the normal and the shear stress

for ic = 1:numpts
    nvec = n(ic,:);
    P = sigma * nvec';
    Pmag = norm(P);
    normal = nvec * P;
    nmag = norm(normal);
    shearup =sqrt(Pmag^2-nmag^2);
    sheardown=-shearup;
    sn(ic)=normal;
    ss1(ic)=shearup;
    ss2(ic)=sheardown;
end

plot(sn,ss1,'.',sn,ss2,'.')
axis([ -10 20 -15 15])
axis square
grid