clc; clear all; close all;

L = 600*10^(-2);
E = 50000*10^(7);
I = 30000 * 10^(-8);
w0 = 2.5 * 10^(5);
dx = 0.01; % 원래는 10cm 이어야 하지만 전진차분법으로 좀더 정확히 계산하기 위해 요소를 1cm로 쪼갬
x = 0:dx:L;
y = (w0/(120*E*I*L))*(-x.^5+2*L^2*x.^3-L^4.*x);
s = length(y);


%% Slope
for i=1:s-1
slope(i) = (y(i+1)-y(i))./dx;
end
slope(s) = slope(s-1); % 전진차분법이므로 최종값은 직전값과 동일하게 어림잡아 계산

%% Moment
for i=1:s-1
Moment(i) = E*I*(slope(i+1)-slope(i))./dx;

 % 조건문 : 이전함수의 값이 직전값과 동일하게 시작될 경우 기울기는 0이 되므로 그것을 피하기 위해 그 이후부터는 같은 값을 대입
     if slope(i+1) == slope(i)  
        for t = i:s
            Moment(t) = Moment(t-1);
        end
    end
    
end

%% Shear Force
for i=1:s-1
Shear(i) = (Moment(i+1)-Moment(i))./dx;

    if Moment(i+1) == Moment(i)
        for t = i:s
            Shear(t) = Shear(t-1);
        end
    end
    
end

%% Loading
for i=1:s-1
Loading(i) = -(Shear(i+1)-Shear(i))./dx;

    if Shear(i+1) == Shear(i)
        for t = i:s
            Loading(t) = Loading(t-1);
        end
    end
    
end

subplot(5,1,1);
plot(x,y);
axis([0 L -6*10^(-3) 0 ]);
xlabel('beam location(m)');
ylabel('Displacement(m)');
title('Displacement');
grid on;

subplot(5,1,2);
plot(x,slope);
xlabel('beam location(m)');
ylabel('Slope(rad)');
title('Slope(rad)');
grid on;

subplot(5,1,3)
plot(x,Moment);
xlabel('beam location(m)');
ylabel('Moment(N*m)');
title('Moment(N*m)');
grid on;

subplot(5,1,4)
plot(x,Shear);
xlabel('beam location(m)');
ylabel('Shear(N)');
title('Shear');
grid on;

subplot(5,1,5)
plot(x,Loading);
xlabel('beam location(m)');
ylabel('loading(N/m)');
title('loading(N/m)');
grid on;