clear all;
%Построение упругой линии балки при прямом изгибе методом начальных параметров.
%Шарнирно-опертая балка под действием центральной сосредоточенной силы.

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

%Начальные параметры для первого участка:
v_0 = 0;                % м, прогиб               
fi_0 = -P*a^2/(4*EI);   % рад,  угол поворота 
M_0 = 0;                % Нм, момент
Q_0 = P/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) - P;        % Н, сила
q_0 = 0;                % Н/м, интенсивность распределенной нагрузки

%Второй участок
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;

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





