Why is my projectile motion code only working at certain input angles.

1 view (last 30 days)
This projectile motion is supposed to project the distance of a baseball with drag consdered. It is not working at some input angles but it is at others. For some reason i can get a good output reading when my input is 51 m/s and 15 deg. Ive gotten a few other random ggod readings but i'm not sure why it's like that so any help would be greatly appreciated.
clear
close all
V = input('input exit velocity [m/s]: ');
launch_angle = input('input launch angle [deg]: ');
Vx = V * cos(launch_angle);
Vy = V * sin(launch_angle);
%intial conditions
G = 9.80665; %m/s^2 Acceleration due to Gravity
DC = 0.4; %Drag coefficient
Area = 0.00426; %m^2 cross section area of ball
Mass = 0.149; %aKg mass of ball
x(1) = 0; %intial x postion
y(1) = 1.2; %inital y postion
xf(1) = 0; %inital xf postion
yf(1) = 0; %intial yf postion
AP = 1.2; %kg/m^3 Air Density @ Sea Level
D = AP*DC*Area/2; %constant for drag calculations
t(1) = 0; %sets intial time
dt = 0.01; %s set the intervals at which time will be evalutated
i = 1; %sets counter/index
while min(y)> -0.01;
t = t + dt;
i = i + 1;
xf(i) = xf(i-1)+ Vx.*dt;
AxD = - ( D / Mass ) * V * Vx
AyD = -G - ( D / Mass ) * V * Vy
Vx_new = Vx + AxD * dt;
Vy_new = Vy + AyD * dt;
x(i) = x(i-1) + Vx_new * dt + 0.5 * AxD * dt^2;
y(i) = y(i-1) + Vy_new * dt + 0.5 * AyD * dt^2;
Vx = Vx_new;
Vy = Vy_new;
end;
proj_distance = -xf(i).*3.28;
disp(proj_distance)
  3 Comments
Nicola Pugliese
Nicola Pugliese on 14 Apr 2022
I think I want the last reading because I only want the output to be the final distance. And yes that definitely had to be fixed with my sin and cos but the outputs are still off for some reason.
Voss
Voss on 14 Apr 2022
That makes sense.
Maybe see if you can get it to look right with no drag (D = 0) first.

Sign in to comment.

Answers (2)

Mathieu NOE
Mathieu NOE on 14 Apr 2022
hello
made a for loop on angles to test the code (for one given V velocity) and it worked as soon as I put the "clear x y" line in the code
clear
close all
% V = input('input exit velocity [m/s]: ');
% launch_angle = input('input launch angle [deg]: ');
V = 10;
figure(1), hold on
launch_angle = (5:5:85);
for ci = 1:numel(launch_angle)
clear x y % here <=
Vx = V * cosd(launch_angle(ci));
Vy = V * sind(launch_angle(ci));
%intial conditions
G = 9.80665; %m/s^2 Acceleration due to Gravity
DC = 0.4; %Drag coefficient
Area = 0.00426; %m^2 cross section area of ball
Mass = 0.149; %aKg mass of ball
x(1) = 0; %intial x postion
y(1) = 1.2; %inital y postion
xf(1) = 0; %inital xf postion
yf(1) = 0; %intial yf postion
AP = 1.2; %kg/m^3 Air Density @ Sea Level
D = AP*DC*Area/2; %constant for drag calculations
t(1) = 0; %sets intial time
dt = 0.01; %s set the intervals at which time will be evalutated
i = 1; %sets counter/index
while min(y)> -0.01
t = t + dt;
i = i + 1;
xf(i) = xf(i-1)+ Vx.*dt;
AxD = - ( D / Mass ) * V * Vx;
AyD = -G - ( D / Mass ) * V * Vy;
Vx_new = Vx + AxD * dt;
Vy_new = Vy + AyD * dt;
x(i) = x(i-1) + Vx_new * dt + 0.5 * AxD * dt^2;
y(i) = y(i-1) + Vy_new * dt + 0.5 * AyD * dt^2;
Vx = Vx_new;
Vy = Vy_new;
end;
proj_distance = -xf(i).*3.28;
disp(proj_distance)
legstr{ci} = (['Angle :' num2str(launch_angle(ci)) ]);
plot(x,y)
end
legend(legstr);
hold off
  1 Comment
Nicola Pugliese
Nicola Pugliese on 15 Apr 2022
Thank you this works much better but my results are still a bit inacurrate. Thats most likely a mathematical issue though.

Sign in to comment.


James Tursa
James Tursa on 15 Apr 2022
Edited: James Tursa on 15 Apr 2022
Drag depends on current velocity, not initial velocity. So you need to recalculate V at each step. E.g.,
Vx = Vx_new;
Vy = Vy_new;
V = norm([Vx Vy]);
Also it looks like you are double banging the acceleration into the position solution. The acceleration gets into velocity, which you then integrate into delta position via the term. But you also integrate it directly into position via the term. So it is getting into the position solution twice. For now, I would suggest doing simple Euler integration for everything. E.g.,
x(i) = x(i-1) + Vx * dt;
y(i) = y(i-1) + Vy * dt;
Vx = Vx + AxD * dt;
Vy = Vy + AyD * dt;
V = norm([Vx Vy]);
No need for the Vx_new and Vy_new variables.
Once you get the simple Euler scheme working, you can explore more advanced schemes such as Modified Euler, RK4, or even ode45( ).
  2 Comments
Nicola Pugliese
Nicola Pugliese on 15 Apr 2022
This works better but the results that I am getting are still very innacurate. Could it be something wrong with the code mathematically? This is where I am at right now with it.
clear
close all
Vmph = input('input exit velocity [mph]: ');
V = Vmph./2.237; %convert mph to m/s
launch_angle = input('input launch angle [deg]: ');
figure(1), hold on
for ci = 1:numel(launch_angle)
clear x y
Vx = V * cosd(launch_angle(ci));
Vy = V * sind(launch_angle(ci));
%intial conditions
G = 9.80665; %m/s^2 Acceleration due to Gravity
DC = 0.4; %Drag coefficient
Area = 0.00426; %m^2 cross section area of ball
Mass = 0.149; %aKg mass of ball
x(1) = 0; %intial x postion
y(1) = 1.2; %inital y postion
xf(1) = 0; %inital xf postion
yf(1) = 0; %intial yf postion
AP = 1.2; %kg/m^3 Air Density @ Sea Level
D = AP*DC*Area/2; %constant for drag calculations
t(1) = 0; %sets intial time
dt = 0.01; %s set the intervals
i = 1; %sets counter/index
while min(y)> -0.01
t = t + dt;
i = i + 1;
xf(i) = xf(i-1)+ Vx.*dt;
AxD = - ( D / Mass ) * V * Vx;
AyD = -G - ( D / Mass ) * V * Vy;
x(i) = x(i-1) + Vx * dt;
y(i) = y(i-1) + Vy * dt;
Vx = Vx + AxD * dt;
Vy = Vy + AyD * dt;
V = norm([Vx Vy]);
end;
xfinal = xf(i).*3.28;
disp(xfinal)
legstr{ci} = (['Angle :' num2str(launch_angle(ci)) ]);
plot(x,y)
end
legend(legstr);
hold off
Mathieu NOE
Mathieu NOE on 19 Apr 2022
Edited: Mathieu NOE on 19 Apr 2022
hello
what makes you believe there is something wrong mathematically ? the longest range is obtained with an initial angle of 45° which is the correct answer.
I looked again at the drag force equation and the projection on x and y axes are correct (IMO)
beside that , minor upgrade of the code as constants don't need to be reinitialized at every loop
clear
close all
% V = input('input exit velocity [m/s]: ');
% launch_angle = input('input launch angle [deg]: ');
V = 10;
figure(1), hold on
launch_angle = (5:5:85);
% constants
G = 9.80665; %m/s^2 Acceleration due to Gravity
Cd = 0.4; %Drag coefficient
Area = 0.00426; %m^2 cross section area of ball
Mass = 0.149; %aKg mass of ball
rho_air = 1.2; %kg/m^3 Air Density @ Sea Level
D = rho_air*Cd*Area/2; %constant for drag calculations
dt = 0.01; %s set the intervals at which time will be evalutated
for ci = 1:numel(launch_angle)
clear x y
Vx = V * cosd(launch_angle(ci));
Vy = V * sind(launch_angle(ci));
%intial conditions
x(1) = 0; %intial x postion
y(1) = 1.2; %inital y postion
t(1) = 0; %sets intial time
i = 1; %sets counter/index
while min(y)> -0.01
t = t + dt;
i = i + 1;
AxD = - ( D / Mass ) * V * Vx;
AyD = -G - ( D / Mass ) * V * Vy;
Vx = Vx + AxD * dt;
Vy = Vy + AyD * dt;
V = sqrt(Vx.^2 + Vy.^2);
x(i) = x(i-1) + Vx * dt;
y(i) = y(i-1) + Vy * dt;
end
proj_distance = x(i).*3.28;
disp(proj_distance)
legstr{ci} = (['Angle :' num2str(launch_angle(ci)) ]);
plot(x,y)
end
legend(legstr);
hold off

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!