%《MATLAB及其在理工課程中的應用指南》第七章例7－2－4程序
% 復雜應力狀態的分析，摩爾圓
% 西安電子科技大學出版社出版  陳懷琛編著 1999年10月初版，2004年11月第二版
%
Sx=input('σx(Mpa) =(書上取20)  '); % 輸入兩方向的正應力和剪應力,S=σ,T=τ
Sy=input('σy(Mpa) =(書上取0)  '); 
Txy=input('τxy(Mpa) =(書上取5)  ');
a=linspace(0,pi,37);		        % 將應力圓上的圓周角分為36份
Sa=(Sx+Sy)/2; Sd=(Sx-Sy)/2;
sigma=Sa+Sd*cos(2*a)-Txy*sin(2*a);	% 應力圓方程
tau = Sd*sin(2*a)+Txy*cos(2*a);
plot(sigma,tau,Sx,Txy,'*'),		    % 繪圖并繪出基準點，
set(gcf,'color','w')                % 置圖形背景色為白色
axis equal, v=axis;		            % 使x和y軸取等比例,找坐標軸的邊界v
line([v(1),v(2)],[0,0]),	        % 畫出x,y坐標,
line([0,0],[v(3),v(4)])						
hold, plot(Sa,0,'x')				% 標出應力圓的圓心
Smax=max(sigma),Smin=min(sigma),Tmax=max(tau)
% 以下是求任意斜截面a角上σ及τ的程序段
h=input('若不求應力,鍵入0.若求應力,鍵入1   ');
while h~=0;
   a=input('給出斜截面方向角a= [弧度](書上取pi/3或-pi/5) ')
   sigma=Sa+Sd*cos(2*a)-Txy*sin(2*a)
tau=Sd*sin(2*a)+Txy*cos(2*a)
plot(sigma,tau,'or')				% 畫出待求的應力狀態點
h=input('若不繼續求應力,鍵入0.若還要求應力,鍵入1');
end,hold off