%script to visualize stress tensor and traction vector in physical space
%JRA, March 6, 2000

%Arguments:
sigma1 = 100;
sigma3 = 50;
%theta = 65;

for theta=[0:10:180]

%Calculations:
thetarads=theta*pi/180;
[sigman, tau] = mohr(sigma1,sigma3,theta);

%plotting in physical space
figure(1)
clf 
boxscale = 2;
%%Plot box
%%%outline:
xbox = [-1 1 1 -1 -1];
xbox = boxscale.*xbox;
ybox = [-1 -1 1 1 -1];
ybox = boxscale.*ybox;
plot(xbox, ybox, 'k-');
axis([-2*boxscale 2*boxscale -2*boxscale 2*boxscale]);
title('Physical Space')
hold on

%%Plot plane
initialx=[0 0];
initialy=[1 -1];
initialy=boxscale.*initialy;
rotx=initialx*cos(thetarads)-initialy*sin(thetarads);
roty=initialx*sin(thetarads)+initialy*cos(thetarads);
plot(rotx,roty,'k-')

%%Plot pole
initialx=[0 0];
initialy=[0 -1];
initialx=0.25*boxscale.*initialx;
rotx=initialx*cos(thetarads+pi/2)-initialy*sin(thetarads+pi/2);
roty=initialx*cos(thetarads+pi/2)+initialy*cos(thetarads+pi/2);
plot(rotx,roty,'k-')

%For stresses and tractions, assume that sigma1 is equal to boxscale
%Plot stresses
%%sigma1
sigma1rightx = [boxscale 2*boxscale];
sigma1leftx =  [-boxscale -2*boxscale];
plot(sigma1rightx, [0 0], 'b-')
plot(sigma1leftx, [0 0], 'b-')
%%sigma3
sigma3topy = [boxscale boxscale+(sigma3/sigma1)*boxscale];
sigma3bottomy =  [-boxscale -boxscale-(sigma3/sigma1)*boxscale];
plot([0 0], sigma3topy, 'b-')
plot([0 0], sigma3bottomy, 'b-')

%Plot tractions
%%sigman
initialx=[0 0];
initialy=[0 -1];
initialy=(sigman/sigma1)*boxscale.*initialy;
rotx=initialx*cos(thetarads+pi/2)-initialy*sin(thetarads+pi/2);
roty=initialx*cos(thetarads+pi/2)+initialy*cos(thetarads+pi/2);
hsign =plot(rotx,roty,'r-')
set(hsign,'LineWidth',2)

%%tau
initialx=[0 0];
initialy=[0 1];
initialy=(tau/sigma1)*boxscale.*initialy;
rotx=initialx*cos(thetarads)-initialy*sin(thetarads);
roty=initialx*cos(thetarads)+initialy*cos(thetarads);
htau = plot(rotx,roty,'r-')
set(htau,'LineWidth',2)

%label stresses and tractions
inc = 0.25*boxscale;
xinc =2*boxscale;
s=sprintf('sigma1 = %.1f\n',sigma1)
text(xinc, boxscale, s);
s=sprintf('sigma3 = %.1f\n',sigma3)
text(xinc, boxscale-inc, s);
s=sprintf('theta = %.1f\n',theta)
text(xinc, boxscale-2*inc, s);
s=sprintf('sigman = %.1f\n',sigman)
text(xinc, boxscale-6*inc, s);
s=sprintf('tau = %.1f\n',tau)
text(xinc, boxscale-7*inc, s);


axis equal
axis off
drawnow


%Plot in Mohr space
figure(2)
clf
sigmann = [];
tauall = [];
for thetaforplot = 0:5:360
     [sigman,tau] = mohr(sigma1,sigma3,thetaforplot);
     sigmann = [sigmann sigman];
     tauall = [tauall tau];
end

plot(sigmann, tauall,'k-')
axis equal
axis([0 sigma1*2 -sigma3*2 sigma3*2])
xlabel('sigman');ylabel('tau'); title('Mohr Space');
hold on
plot([0 sigma1*2], [0 0], 'k-')

%Plot the traction pair
[sigman,tau] = mohr(sigma1,sigma3,theta);
plot(sigman,tau, 'ro')
drawnow


end