% stałe tłumienie c2 zmienne c1
clc; close all; clear variables;
%% Dane
m1 = 24;        % [kg]      masa
m2 = 35;        % [kg]      masa
k1 = 1e5;       % [N/m]     sztywność sprężyny
k2 = 18e3;      % [N/m]     sztywność sprężyny
A  = 3.5e-3;    % [m]       amplituda wymuszenia
c1 = [];        % [kg/s]    tłumienie tłumika
c2 = 50;        % [kg/s]    tłumienie tłumika

tsim = 10;      % [s]       czas symulacji
tsam = 0.001;   % [s]       czas próbkowania
% zakres zmienności parametrów
f_v = 0.2:0.1:16;
c1_v = 50:200:1650;
% preinicjalizacja wektorów
T1 = zeros(length(c1_v),length(f_v));
T2 = zeros(length(c1_v),length(f_v));
%% Pętla główna obliczeń
for j = 1:length(c1_v)
    for i = 1:length(f_v)
        f = f_v(i);
        c1 = c1_v(j);
        A1 = 2*A*pi*f;
        sim('sym3_1');
        T1(j,i) = rms(z1)./rms(w);
        T2(j,i) = rms(z2)./rms(w);
    end
end
%%
figure(1)
plot(f_v,T1);
figure(2)
plot(f_v,T2);

%% Wykreślenie wyników
% legenda
for i = 1:length(c1_v)
        leg(i) = {['c1 = ' num2str(c1_v(i)) ' kg/s']};
end
for q = 1:2
    fig = figure(q);
    switch q
    case 1
        fig = figure(q);
        plot(f_v, T1,'LineWidth', 1.2);
        tit = title({'Współczynnik przenoszenia drgań Twz1 między wymuszeniem w a masą m1';...
                     'dla modelu liniowego przy stałym współczynniku c2 = 50 kg/s'});
        legend(leg);
        XLab = xlabel('f [Hz]');
        YLab = ylabel('T_{wz1} [-]');
        figname = 'lab3_1_T1';
    case 2
        fig = figure(q);
        plot(f_v, T2,'LineWidth', 1.2);
        tit = title({'Współczynnik przenoszenia drgań Twz2 między wymuszeniem w a masą m2';...
                     'dla modelu liniowego przy stałym współczynniku c2 = 50 kg/s'});
        legend(leg);
        XLab = xlabel('f [Hz]');
        YLab = ylabel('T_{wz2} [-]');
        figname = 'lab3_1_T2';
    end

    ax = gca;                       grid on; axis tight;
    set(fig,'Color','w');
    set(tit,'FontName','Times',     'FontSize',24);
    set(XLab,'FontName','Times',    'FontSize',22);
    set(YLab,'FontName','Times',    'FontSize',22);
    set(ax,'FontName','Times',      'FontSize',18);
    ax.GridLineStyle = '-';
    ax.GridColor = 'k';
    ax.GridAlpha = 1; % maximum line opacity
    set(gcf, 'Position', [100 100 1100 600]);
    set(gca,'LooseInset',get(gca,'TightInset'))
    print(gcf,figname,'-dpng','-r300');
end