% Clear workspace and command window
clear
clc

%% Trasmissibilita: forze trasmesse al basamento

ground_stiffness = 44842.2; % rigidezza secante [N/m]

susp_mass = 291.75;         % massa sospesa ant sx [kg]

wn = sqrt(ground_stiffness/susp_mass);

a=0.1:0.01:3;
w = wn .* a;
r2 = w.^2. / wn^2;
r = w./wn;

figure;
hold on;

% Range valori zita
zita_range = 0.3:0.05:0.4;

for i = 1:length(zita_range)
    zita = zita_range(i);
    
    % equazione trasmissibilita
    t = sqrt((2 .* zita .* r).^2 + 1) ./ sqrt((1 - r2).^2 + (2 .* zita .* r).^2);
    
    plot(r2, t, 'LineWidth', 1.5);
end

%% Trova intersezione

for i = 1:length(t)
    if t(i) > 0.1 && t(i) < 1
        intersection_point = r2(i);
        break;
    end
end

text(intersection_point + 0.5, 1, ['\omega^2 / \omega_n^2 = ', num2str(intersection_point)], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontSize', 11);
text(intersection_point + 0.5, 1, 'τ = 1', 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left', 'FontSize', 11);

%% Plot
xlabel('\omega^2 / \omega_n^2');
ylabel('τ');
title('Transmissibility relative to. \omega^2 / \omega_n^2');
grid on;
legend(cellstr(num2str(zita_range', 'ζ = %0.2f')));