Why is my projectile motion code only working at certain input angles.
Show older comments
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
Voss
on 14 Apr 2022
Use cosd and sind when the angle is in degrees:
launch_angle = input('input launch angle [deg]: ');
Vx = V * cosd(launch_angle);
Vy = V * sind(launch_angle);
And I think you want to convert all xf at the end, not only the last element of xf:
proj_distance = -xf.*3.28;
(This is not an exhaustive list of things that might need to be fixed.)
Nicola Pugliese
on 14 Apr 2022
Voss
on 14 Apr 2022
That makes sense.
Maybe see if you can get it to look right with no drag (D = 0) first.
Answers (2)
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
on 15 Apr 2022
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.,
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
on 15 Apr 2022
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
Categories
Find more on Creating, Deleting, and Querying Graphics Objects in Help Center and File Exchange
on 14 Apr 2022
on 19 Apr 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
