흙의 압밀비배수 3축압축 시험데이터를 통해 모어원을 그리고 최종적으로 파괴포락선을 그리는 건데요.
요 밑에 -밑줄 - 아래에 있는게 제가 작성한건데 잘 되고 값도 얼추 손으로 구한거랑도 맞고요.
제가 쓰려고 만들었는데 생각해보니깐 후배들 그래프 작도가 엉망이라 지적하기 보단 걍 이걸 활용해라 하고 던져주려고 합니다.
근데 이게 중간에 solve 기능? 함수? 암튼 그걸로 좌표값을 찾다보니 매트랩에서는 간단한데
이걸 옥타브로 하려고 보니 symbolic설치->python 설치->sympy 설치를 해야하는데...
이걸 후배들이 할 수 있을 것 같지가 않아요.(옥타브용으로 수정은 했는데...)
그래서 프로그래머 기사님들 중 아시는 분께 여쭤봅니다.
     1. 심볼릭-solve 기능을 사용하지 않고 좌표를 구할 수 있는 방법?
또는 2. 원의 공통외접선을 구하는 다른 방법이 있는지?(전 도해적으로 구하는 방법으로 구했습니다.)

--------------------------------------밑줄--------------------------------------
minPS=[10, 20, 30]; % 최소주응력
maxPS=[25, 47.3, 67.5]; % 최대주응력
u=[2, 4, 8]; % 과잉간극수압

leng=length(minPS); % 데이터 수
n=100; % 모어원의 x좌표 수
rad=linspace(0,2*pi,n); % 원의 내부각(360도)을 100등분

    PSdif=maxPS-minPS; % 주응력차(전응력=유효응력)
    PSsum=maxPS+minPS; % 주응력합(전응력)
    PSsum1=(maxPS-u)+(minPS-u); % 주응력합(유효응력)
    t=PSdif/2; % 시료의 최대 전단강도(전응력=유효응력)
    tx=maxPS-t; % 시료의 최대 전단강도에서 연직응력(전응력)
    tx1=(maxPS-u)-t; % 시료의 최대 전단강도에서 연직응력(유효응력)

% 모어원의 x,y 좌표 계산
for i=1:leng
        if i==1
            x1=t(i)*cos(rad)+tx(i);
            xe1=t(i)*cos(rad)+tx1(i);
            y1=t(i)*sin(rad);
    elseif i==2
            x2=t(i)*cos(rad)+tx(i);
            xe2=t(i)*cos(rad)+tx1(i);
            y2=t(i)*sin(rad);
    else
            x3=t(i)*cos(rad)+tx(i);
            xe3=t(i)*cos(rad)+tx1(i);
            y3=t(i)*sin(rad);
    end
end

% 파괴포락선 계산
for i=1:leng-1
slope(i)=(t(i+1)-t(i))/(tx(i+1)-tx(i)); % i번째 시료와 i+1번째 시료의 기울기1
eslope(i)=(t(i+1)-t(i))/(tx1(i+1)-tx1(i));

intercept(i)=t(i)-slope(i)*tx(i);
tempx=intercept(i)/slope(i)*-1; % 기울기1과 절편값1으로부터 y=0인 x좌표1 계산
eintercept(i)=t(i)-eslope(i)*tx1(i);
tempx1=eintercept(i)/eslope(i)*-1;

r(i)=(tx(i)-tempx)/2; % i번째 시료의 최대전단강도에서 x좌표와 x좌표1 사이를 지름으로 하는 원의 반지름 계산
rx(i)=tx(i)-r(i); % i번째 시료의 최대전단강도에서 x좌표와 x좌표1 사이를 지름으로 하는 원의 x좌표2 계산
er(i)=(tx1(i)-tempx1)/2;
erx(i)=tx1(i)-er(i);



syms xx yy xx1 yy1
sq=solve((xx-tx(i))^2+(yy)^2-t(i)^2,(xx-rx(i))^2+(yy)^2-r(i)^2);
% i번째 시료의 모어원과 i번째 시료의 최대전단강도에서 x좌표와 x좌표1 사이를 지름으로 하는 원의 교점 계산
sq1=solve((xx1-tx1(i))^2+(yy1)^2-t(i)^2,(xx1-erx(i))^2+(yy1)^2-er(i)^2);
ipx=double(sq.xx); % 교점의 x좌표
ipy=double(sq.yy); % 교점의 y좌표
ipxe=double(sq1.xx1);
ipy2=double(sq1.yy1);

slope2(i)=ipy(2)/(ipx(1)-tempx); % 교점의 x,y좌표로부터 새로운 기울기2 계산
intercept2(i)=slope2(i)*tempx*-1; % 기울기2로부터 절편값2 계산
eslope2(i)=ipy2(2)/(ipxe(1)-tempx1);
eintetcept2(i)=eslope2(i)*tempx1*-1;
end

x4=[0:1:70]; % 파괴포락선의 출력을 위한 x좌표 생성
avgslope=(slope2(1)+slope2(2))/2; % 기울기2의 평균값 계산
avgeslope=(eslope2(1)+eslope2(2))/2;
avgintercept=(intercept2(1)+intercept2(2))/2; % 절편값2의 평균값 계산
avgeintercept=(eintetcept2(1)+eintetcept2(2))/2;
k=(180/pi)*atan(avgslope); % 파괴포락선의 내부마찰각
k1=(180/pi)*atan(avgeslope);


%결과값 출력
hold on
plot(x1,y1,'b-')
plot(x2,y2,'b-')
plot(x3,y3,'b-')
plot(xe1,y1,'b--')
plot(xe2,y2,'b--')
plot(xe3,y3,'b--')
plot(x4,avgslope*x4+avgintercept,'r-')
plot(x4,avgeslope*x4+avgeintercept,'r--')
    
axis([0,70,0,70])
grid on
set(gca,'fontsize',14)
title('CU-Test') % 그래프 제목 출력
xlabel('Normal stress [tf/m^2]') % x축 래이블 출력
ylabel('Shear stress [tf/m^2]') % y축 래이블 출력

abox=[0,0,43,43,0; 55,70,70,55,55]; % 출력 값 박스 생성을 위한 좌표 값
plot(abox(1,:), abox(2,:), 'b-'); % 출력 값 박스 작도(파란 사각형)

for i=1:4 % 박스 내에 c,Φ 값 출력을 위한 for루프
     
   if i==1
       s=['c = ',num2str(avgintercept),'tf/m^2']; % text함수 사용을 위해 숫자를 문자로 변환하여 S변수에 저장      
   elseif i==2
       s=['Φ = ',num2str(k),'deg'];
   elseif i==3
       s=['c` = ',num2str(avgeintercept),'tf/m^2']; % text함수 사용을 위해 숫자를 문자로 변환하여 S변수에 저장      
   else
       s=['Φ` = ',num2str(k1),'deg'];  
   end
   
   if i<=2
        text(2,70-5*i,s,'fontsize',12); % 그래프 상에 결과 값 출력
   else
        text(22,70-5*(i-2),s,'fontsize',12);
   end
end

[결과물]