흙의 압밀비배수 3축압축 시험데이터를 통해 모어원을 그리고 최종적으로 파괴포락선을 그리는 건데요.<br>요 밑에 -밑줄 - 아래에 있는게 제가 작성한건데 잘 되고 값도 얼추 손으로 구한거랑도 맞고요.<br>제가 쓰려고 만들었는데 생각해보니깐 후배들 그래프 작도가 엉망이라 지적하기 보단 걍 이걸 활용해라 하고 던져주려고 합니다.<br>근데 이게 중간에 solve 기능? 함수? 암튼 그걸로 좌표값을 찾다보니 매트랩에서는 간단한데<br>이걸 옥타브로 하려고 보니 symbolic설치->python 설치->sympy 설치를 해야하는데...<br>이걸 후배들이 할 수 있을 것 같지가 않아요.(옥타브용으로 수정은 했는데...)<br>그래서 프로그래머 기사님들 중 아시는 분께 여쭤봅니다.<br>     1. 심볼릭-solve 기능을 사용하지 않고 좌표를 구할 수 있는 방법?<br>또는 2. 원의 공통외접선을 구하는 다른 방법이 있는지?(전 도해적으로 구하는 방법으로 구했습니다.)<br><br>--------------------------------------밑줄--------------------------------------<br>minPS=[10, 20, 30]; % 최소주응력<br>maxPS=[25, 47.3, 67.5]; % 최대주응력<br>u=[2, 4, 8]; % 과잉간극수압<br><br>leng=length(minPS); % 데이터 수<br>n=100; % 모어원의 x좌표 수<br>rad=linspace(0,2*pi,n); % 원의 내부각(360도)을 100등분<br><br>    PSdif=maxPS-minPS; % 주응력차(전응력=유효응력)<br>    PSsum=maxPS+minPS; % 주응력합(전응력)<br>    PSsum1=(maxPS-u)+(minPS-u); % 주응력합(유효응력)<br>    t=PSdif/2; % 시료의 최대 전단강도(전응력=유효응력)<br>    tx=maxPS-t; % 시료의 최대 전단강도에서 연직응력(전응력)<br>    tx1=(maxPS-u)-t; % 시료의 최대 전단강도에서 연직응력(유효응력)<br><br>% 모어원의 x,y 좌표 계산<br>for i=1:leng<br>        if i==1<br>            x1=t(i)*cos(rad)+tx(i);<br>            xe1=t(i)*cos(rad)+tx1(i);<br>            y1=t(i)*sin(rad);<br>    elseif i==2<br>            x2=t(i)*cos(rad)+tx(i);<br>            xe2=t(i)*cos(rad)+tx1(i);<br>            y2=t(i)*sin(rad);<br>    else<br>            x3=t(i)*cos(rad)+tx(i);<br>            xe3=t(i)*cos(rad)+tx1(i);<br>            y3=t(i)*sin(rad);<br>    end<br>end<br><br>% 파괴포락선 계산<br>for i=1:leng-1<br>slope(i)=(t(i+1)-t(i))/(tx(i+1)-tx(i)); % i번째 시료와 i+1번째 시료의 기울기1<br>eslope(i)=(t(i+1)-t(i))/(tx1(i+1)-tx1(i));<br><br>intercept(i)=t(i)-slope(i)*tx(i);<br>tempx=intercept(i)/slope(i)*-1; % 기울기1과 절편값1으로부터 y=0인 x좌표1 계산<br>eintercept(i)=t(i)-eslope(i)*tx1(i);<br>tempx1=eintercept(i)/eslope(i)*-1;<br><br>r(i)=(tx(i)-tempx)/2; % i번째 시료의 최대전단강도에서 x좌표와 x좌표1 사이를 지름으로 하는 원의 반지름 계산<br>rx(i)=tx(i)-r(i); % i번째 시료의 최대전단강도에서 x좌표와 x좌표1 사이를 지름으로 하는 원의 x좌표2 계산<br>er(i)=(tx1(i)-tempx1)/2;<br>erx(i)=tx1(i)-er(i);<br><br><br><br>syms xx yy xx1 yy1<br>sq=solve((xx-tx(i))^2+(yy)^2-t(i)^2,(xx-rx(i))^2+(yy)^2-r(i)^2);<br>% i번째 시료의 모어원과 i번째 시료의 최대전단강도에서 x좌표와 x좌표1 사이를 지름으로 하는 원의 교점 계산<br>sq1=solve((xx1-tx1(i))^2+(yy1)^2-t(i)^2,(xx1-erx(i))^2+(yy1)^2-er(i)^2);<br>ipx=double(sq.xx); % 교점의 x좌표<br>ipy=double(sq.yy); % 교점의 y좌표<br>ipxe=double(sq1.xx1);<br>ipy2=double(sq1.yy1);<br><br>slope2(i)=ipy(2)/(ipx(1)-tempx); % 교점의 x,y좌표로부터 새로운 기울기2 계산<br>intercept2(i)=slope2(i)*tempx*-1; % 기울기2로부터 절편값2 계산<br>eslope2(i)=ipy2(2)/(ipxe(1)-tempx1);<br>eintetcept2(i)=eslope2(i)*tempx1*-1;<br>end<br><br>x4=[0:1:70]; % 파괴포락선의 출력을 위한 x좌표 생성<br>avgslope=(slope2(1)+slope2(2))/2; % 기울기2의 평균값 계산<br>avgeslope=(eslope2(1)+eslope2(2))/2;<br>avgintercept=(intercept2(1)+intercept2(2))/2; % 절편값2의 평균값 계산<br>avgeintercept=(eintetcept2(1)+eintetcept2(2))/2;<br>k=(180/pi)*atan(avgslope); % 파괴포락선의 내부마찰각<br>k1=(180/pi)*atan(avgeslope);<br><br><br>%결과값 출력<br>hold on<br>plot(x1,y1,'b-')<br>plot(x2,y2,'b-')<br>plot(x3,y3,'b-')<br>plot(xe1,y1,'b--')<br>plot(xe2,y2,'b--')<br>plot(xe3,y3,'b--')<br>plot(x4,avgslope*x4+avgintercept,'r-')<br>plot(x4,avgeslope*x4+avgeintercept,'r--')<br>    <br>axis([0,70,0,70])<br>grid on<br>set(gca,'fontsize',14)<br>title('CU-Test') % 그래프 제목 출력<br>xlabel('Normal stress [tf/m^2]') % x축 래이블 출력<br>ylabel('Shear stress [tf/m^2]') % y축 래이블 출력<br><br>abox=[0,0,43,43,0; 55,70,70,55,55]; % 출력 값 박스 생성을 위한 좌표 값 <br>plot(abox(1,:), abox(2,:), 'b-'); % 출력 값 박스 작도(파란 사각형)<br><br>for i=1:4 % 박스 내에 c,Φ 값 출력을 위한 for루프<br>     <br>   if i==1<br>       s=['c = ',num2str(avgintercept),'tf/m^2']; % text함수 사용을 위해 숫자를 문자로 변환하여 S변수에 저장      <br>   elseif i==2<br>       s=['Φ = ',num2str(k),'deg'];<br>   elseif i==3<br>       s=['c` = ',num2str(avgeintercept),'tf/m^2']; % text함수 사용을 위해 숫자를 문자로 변환하여 S변수에 저장      <br>   else<br>       s=['Φ` = ',num2str(k1),'deg'];  <br>   end<br>   <br>   if i<=2<br>        text(2,70-5*i,s,'fontsize',12); % 그래프 상에 결과 값 출력<br>   else<br>        text(22,70-5*(i-2),s,'fontsize',12);<br>   end<br>end<br><br>[결과물]<br><div style="text-align:left;"><img width="560" height="420" src="http://thimg.todayhumor.co.kr/upfile/201608/1471373754d6eae016e715438e833e044becdf86e0__w560__h420__f39201__Ym201608.png" alt="cu-test.png" style="border:medium none;"></div><br>