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

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


% Condizioni iniziali
v0=10;
v0_x=v0*cos(alpha);
v0_y=v0*sin(alpha);

z0=[0, v0_x, hin,  v0_y]; %[x, vx, y, vy]


tspan = [0, 5]; % intervallo di tempo per l'integrazione

% Integrazione numerica
%[t, z]=ode45(dvdt, tspan, z0);
[t, z]=ode45(@(t,z) moto(t,z,m,g,A,rho_a, Cd), tspan, z0);

% Risultati finali
x = z(:,1);
y = z(:,3);

% Plot della traiettoria
plot(x, y)
hold on
plot(gittata, hfin, 'r*')
xlabel('Posizione lungo x [m]')
ylabel('Posizione lungo y [m]')
title('Traiettoria del lancio della palla')
grid on