clear
clc

%% Dati da Lotus Suspension Analysis
position = [
40.000
35.000
30.000
25.000
20.000
15.000
10.000
5.000
0.000
-5.000
-10.000
-15.000
-20.000
-25.000
-30.000
-35.000
-40.000
-45.000
-50.000
-55.000
-60.000
-65.000
];

tire_forces = [
4595.7847
4358.7183
4123.7476
3890.7935
3659.7786
3430.6235
3203.2480
2977.5703
2753.5076
2530.9746
2309.8843
2090.1472
1871.6707
1654.3597
1438.1155
1222.8350
1008.4114
794.7330
581.6828
369.1380
156.9698
-54.9575
];

spring_forces = [
6751.6
6421.9
6092.7
5764.0
5435.7
5107.9
4780.5
4453.6
4127.2
3801.1
3475.5
3150.3
2825.5
2501.1
2177.1
1853.4
1529.9
1206.8
883.9
561.2
238.7
-83.6
];

spring_travel = [
26.8343
23.4638
20.0979
16.7366
13.3800
10.0280
6.6807
3.3381
0.0000
-3.3337
-6.6629
-9.9879
-13.3089
-16.6260
-19.9395
-23.2495
-26.5564
-29.8605
-33.1620
-36.4613
-39.7588
-43.0549
];

%% Interpolazioni
% Curva interpolata per la forza del pneumatico
tire_fit_x = min(position):0.1:max(position);
tire_fit_y = interp1(position, tire_forces, tire_fit_x, 'spline');

% Curva interpolata per la forza della molla
spring_fit_x = min(position):0.1:max(position);
spring_fit_y = interp1(position, spring_forces, spring_fit_x, 'spline');

%% Grafico della forza verticale sullo pneumatico in funzione dello
% spostamento della ruota
figure;
hold on
plot(position, tire_forces, '.', tire_fit_x, tire_fit_y, '-r', 'LineWidth', 1.5);

% Seleziona i punti sul grafico: 1 = bump 40mm, 17 = rebound 40mm
x1 = position(1);
x2 = position(17);
y1 = tire_forces(1);
y2 = tire_forces(17);

% Calcola la rigidezza a terra (ground stiffness)
g_s = (y2 - y1) / (x2 - x1);

% Traccia la retta blu tra i due punti -40 e 40
plot([-40, 40], [y2, y1], 'b--', 'LineWidth', 1);

% Posizione per il testo della costante della molla
text_x = mean([x1, x2]);
text_y = mean([y1, y2]);

% Visualizza il valore della costante della molla
text(text_x + 3, text_y + 3, ['k = ', num2str(g_s), ' N/mm'], 'FontSize', 15, 'Color', 'blue', 'VerticalAlignment', 'cap', 'HorizontalAlignment', 'left');

title('Tire Force relative to wheel position');
grid on;
grid minor;
xlabel('Wheel position [mm]');
ylabel('Force [N]');
text(position, tire_forces, string(tire_forces), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
hold off;

%% Grafico della forza sulla molla in funzione dello spostamento della ruota
figure;
plot(position, spring_forces, '.', spring_fit_x, spring_fit_y, 'LineWidth', 1.5);
title('Spring Force relative to wheel position');
grid on;
grid minor;
xlabel('Wheel position [mm]');
ylabel('Force [N]');
text(position, spring_forces, string(spring_forces), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');

%% Grafico della forza sulla molla su forza sulla ruota
figure;
plot(spring_forces, tire_forces, 'LineWidth', 1.5);
title('Spring Force relative to Tire Force');
grid on;
grid minor;
xlabel('Spring Force [N]');
ylabel('Tire Force [N]');

%% Grafico dello spring travel su spring force
figure;
plot(spring_travel, spring_forces, 'LineWidth', 1.5);
hold on;

% Seleziona due punti sul grafico
% 1 = bump 40mm, 9 = condizione statica (auto a terra)
x1 = spring_travel(1);
x2 = spring_travel(9);
y1 = spring_forces(1);
y2 = spring_forces(9);

% Calcola la costante della molla (k)
k = (y2 - y1) / (x2 - x1);

% Posizione per il testo della costante della molla
text_x = mean([x1, x2]);
text_y = mean([y1, y2]);

% Visualizza il valore della costante della molla
text(text_x+4, text_y+4, ['k = ', num2str(k), ' N/mm'], 'FontSize', 15, 'Color', 'blue', 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left');

title('Spring Travel vs Spring Force');
xlabel('Spring Travel [mm]');
ylabel('Spring Force [N]');
grid on;
hold off;

%% Grafico della rigidezza a terra (secante)

% Calcolo del delta della forza della molla rispetto alla condizione di position = 0
delta_spring_forces = tire_forces - tire_forces(9);

% Calcolo della rigidezza per ogni posizione della ruota
stiffness = delta_spring_forces ./ position;

% Curva interpolata per la rigidezza a terra
stiffness_fit_x = min(position):0.1:max(position);
stiffness_fit_y = interp1(position, stiffness, stiffness_fit_x, 'spline');

figure;
plot(position, stiffness, '-', stiffness_fit_x, stiffness_fit_y, 'LineWidth', 1.5);
text(position, stiffness, string(stiffness), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
grid on;
grid minor;
title('Vehicle ground stiffness');
xlabel('Wheel travel [mm]');
ylabel('Stiffness [N/mm^2]');