clear all;
%Сопротивление материалов, ДПМ, семестр 4.
%Построение упругой линии статически-неопределимой балки при прямом изгибе методом начальных параметров.
%Три участка длиной a, при z=0, z=a, z=2a - шарнирно-подвижные опоры,
%правый конец - шарнирно-неподвижная опора. Нагрузка: левый край -
%положительный сосредоточенный момент qa2, средний участок под
%распределенной нагрузкой q.

%Исходные данные:
a = 1;                  % м, длина первого и второго участка балки 
I = 7.080*10^-5;        % м^4, момент инерции сечения
E = 2*10^11;            % Па, модуль Юнга для стали
EI = E*I;               % жесткость при изгибе
q = 1000;               % H/м, распределенная нагрузка

%Система граничных условий (неизвестные - [fi0  Ay  By  Cy]):
A = [a      a^3/(6*EI)      0               0;
    2*a     (4/3)*(a^3/EI)  a^3/(6*EI)      0;
    3*a     (9/2)*(a^3/EI)  (4/3)*(a^3/EI)  a^3/(6*EI);
    0       3*a              2*a            a];
b = [-q*a^4/(2*EI);  -47*q*a^4/(24*EI);  -93*q*a^4/(24*EI); (1/2)*q*a^2];

x = A^-1*b;

%Начальные параметры для первого участка:
v_0 = 0;                % м, прогиб               
fi_0 = x(1);   % рад,  угол поворота 
M_0 = q*a^2;                % Нм, момент
Q_0 = x(2);              % Н, сила
q_0 = 0;                % Н/м, интенсивность распределенной нагрузки

%Первый участок
N1 = 20;                    % число точек  
z1 = linspace(0, a, N1);     % координата 

for i=1:N1
    
    % координата
    z(i) = z1(i);
    
    %поперечная сила
    Q(i) = Q_0 + q_0*z(i);  
    
    %изгибающий момент
    M(i) = M_0 + Q_0*z(i) + q_0*(z(i))^2/2; 
    
    %угол поворота сечения
    fi(i)= fi_0 + M_0*z(i)/EI + Q_0*(z(i))^2/(2*EI) + q_0*(z(i))^3/(6*EI);  
    
    %прогиб
    v(i) = v_0 + fi_0*z(i) + M_0*z(i)*z(i)/(2*EI) + Q_0*(z(i))^3/(6*EI) + q_0*(z(i))^4/(24*EI);  
    
end;

%Начальные параметры для второго участка:
v_0 = v(N1);            % м, прогиб               
fi_0 = fi(N1);          % рад,  угол поворота 
M_0 = M(N1);            % Нм, момент
Q_0 = Q(N1) + x(3);     % Н, сила
q_0 = -q;                % Н/м, интенсивность распределенной нагрузки

%Второй участок
N2 = 20;                     % число точек  
z2 = linspace(a, 2*a, N2);   % координата 

for i=1:N2
    
    ii=i + N1;
    
    %глобальная координата
    z(ii) = z2(i);
    
    %поперечная сила
    Q(ii) = Q_0 + q_0*z(i);  
    
    %изгибающий момент
    M(ii) = M_0 + Q_0*z(i) + q_0*(z(i))^2/2; 
    
    %угол поворота сечения
    fi(ii)= fi_0 + M_0*z(i)/EI + Q_0*(z(i))^2/(2*EI) + q_0*(z(i))^3/(6*EI);  
    
    %прогиб
    v(ii) = v_0 + fi_0*z(i) + M_0*z(i)*z(i)/(2*EI) + Q_0*(z(i))^3/(6*EI) + q_0*(z(i))^4/(24*EI);  
    
end;

%Начальные параметры для третьего участка:
v_0 = v(N1+N2);            % м, прогиб               
fi_0 = fi(N1+N2);          % рад,  угол поворота 
M_0 = M(N1+N2);            % Нм, момент
Q_0 = Q(N1+N2) + x(4);     % Н, сила
q_0 = 0;                   % Н/м, интенсивность распределенной нагрузки

%Третий участок
N3 = 20;                     % число точек  
z3 = linspace(2*a, 3*a, N3);   % координата 

for i=1:N3
    
    ii=i + N1 + N2;
    
    %глобальная координата
    z(ii) = z3(i);
    
    %поперечная сила
    Q(ii) = Q_0 + q_0*z(i);  
    
    %изгибающий момент
    M(ii) = M_0 + Q_0*z(i) + q_0*(z(i))^2/2; 
    
    %угол поворота сечения
    fi(ii)= fi_0 + M_0*z(i)/EI + Q_0*(z(i))^2/(2*EI) + q_0*(z(i))^3/(6*EI);  
    
    %прогиб
    v(ii) = v_0 + fi_0*z(i) + M_0*z(i)*z(i)/(2*EI) + Q_0*(z(i))^3/(6*EI) + q_0*(z(i))^4/(24*EI);  
    
end;


subplot(4,1,1), area(z, Q, 'FaceColor', [0.75 1 0.75]), xlabel('z, m'); ylabel('Q, H'), title([' Поперечная сила Q ']), grid on,
subplot(4,1,2), area(z, M, 'FaceColor', [0.75 1 0.75]); xlabel('z, m'); ylabel('M, Hm'), title([' Изгибающий момент М ']),grid on,
subplot(4,1,3), area(z, fi, 'FaceColor', [0.75 1 0.75]); xlabel('z, m'); ylabel('\phi, rad'), title([' Угол поворота \phi ']),grid on,
subplot(4,1,4), area(z, v, 'FaceColor', [0.75 1 0.75]); xlabel('z, m'); ylabel('v, m'), title([' Прогиб v ']),grid on





