Untitled

% Dati del problema
rho_p = 250; %[kg/m^3]
d = 0.05; %[m]
alpha = pi/12; %[rad]
hin = 2; %[m]
gittata = 18; %[m]
hfin = 1.5; %[m]
g = 9.81; %[m/s^2]
Cd = 0.5; %[-]
rho_a = 1.2; %[kg/m^3]
V_p = pi/6 * d^3; % volume sfera
m = rho_p * V_p; % massa sfera
A = pi * (d/2)^2; % area sezione trasversale della sfera

% Funzione per il moto
function dzdt = moto(t,z,m,g,A,rho_a,Cd)
    % z = [x, vx, y, vy]
    v = sqrt(z(2)^2 + z(4)^2); %velocità totale
    
    dzdt = zeros(4,1);
    dzdt(1) = z(2); %dx/dt = vx
    dzdt(2) = -(1/(2*m))*rho_a*A*Cd*v*z(2); %dvx/dt
    dzdt(3) = z(4); %dy/dt = vy
    dzdt(4) = -g -(1/(2*m))*rho_a*A*Cd*v*z(4); %dvy/dt
end

% Funzione per il sistema di equazioni da risolvere
function err = target_system(x, params)
    % x(1) = tempo di volo
    % x(2) = velocità iniziale
    
    % Estrai parametri
    alpha = params.alpha;
    hin = params.hin;
    m = params.m;
    g = params.g;
    A = params.A;
    rho_a = params.rho_a;
    Cd = params.Cd;
    gittata = params.gittata;
    hfin = params.hfin;
    
    % Calcola componenti velocità iniziale
    v0 = x(2);
    v0_x = v0*cos(alpha);
    v0_y = v0*sin(alpha);
    
    % Condizioni iniziali
    z0 = [0, v0_x, hin, v0_y];
    
    % Integra il moto fino al tempo di volo
    [~, z] = ode45(@(t,z) moto(t,z,m,g,A,rho_a,Cd), [0 x(1)], z0);
    
    % Calcola errori (distanza dal target)
    err = zeros(2,1);
    err(1) = z(end,1) - gittata;  % errore in x
    err(2) = z(end,3) - hfin;     % errore in y
end

% Raggruppa i parametri
params = struct('alpha', alpha, 'hin', hin, 'm', m, 'g', g, 'A', A, ...
    'rho_a', rho_a, 'Cd', Cd, 'gittata', gittata, 'hfin', hfin);

% Risolvi il sistema con fsolve
x0 = [2, 10];  % stima iniziale [tempo di volo, velocità iniziale]
options = optimoptions('fsolve', 'Display', 'iter');
[x_sol, fval] = fsolve(@(x) target_system(x, params), x0, options);

% Estrai soluzioni
tempo_volo = x_sol(1);
v0 = x_sol(2);

% Calcola traiettoria con i valori trovati
v0_x = v0*cos(alpha);
v0_y = v0*sin(alpha);
z0 = [0, v0_x, hin, v0_y];
[t, z] = ode45(@(t,z) moto(t,z,m,g,A,rho_a,Cd), [0 tempo_volo], z0);

% Plot risultati
figure
plot(z(:,1), z(:,3), 'b-', 'LineWidth', 1.5)
hold on
plot(gittata, hfin, 'r*', 'MarkerSize', 10)
plot(0, hin, 'go', 'MarkerSize', 10)
xlabel('Posizione lungo x [m]')
ylabel('Posizione lungo y [m]')
title(['Traiettoria con v0 = ' num2str(v0, '%.2f') ' m/s, tempo = ' num2str(tempo_volo, '%.2f') ' s'])
legend('Traiettoria', 'Target', 'Punto di partenza')
grid on

% Stampa risultati
fprintf('Velocità iniziale: %.2f m/s\n', v0)
fprintf('Tempo di volo: %.2f s\n', tempo_volo)