plotting projectile with drag

26 views (last 30 days)
I have a script to plot the projectile motion path with drag. However, the plot comes out as with no drag. Using parameter v0=150,h=10,dt=0.01,theta=0.1745
I think the mistake is in my while loop, buti can't figure it out.
function [rx,ry,vx,vy]=solve_ode_euler(v0,h,dt,theta)
n=1000; %number of maximum iterations
rho=1.225;
cd=0.479;
m=5.5;
D=0.1;
R=D/2;
A=4*pi*R^4;
g=9.81;
% initializing values of distance and velocity
rx(1)=0;
ry(1)=h;
v(1)=v0;
vx=v0*cos(theta);
vy=v0*sin(theta);
vx(1)=vx;
vy(1)=vy;
k=cd*rho*A/(2*m);
t(1)=0;
i=1;
dt=0.01;
% while loop to solve projectile with Euler method
while i<n
v(i+1)=sqrt((vx(i)^2)+(vy(i)^2)); % velocity magnitude as a vector
rx(i+1)=rx(i)+vx(i)*dt;
vx(i+1)=vx(i)-(k*v(i)*vx(i))*dt;
ry(i+1)=ry(i)+vy(i)*dt;
vy(i+1)=vy(i)-g*dt-(k*v(i)*vy(i))*dt;
t(i+1)=t(i)+dt;
if ry(i+1)<0 %stops the projectile if reaches the ground
i=n;
else
i=i+1;
end
plot(rx,ry)
end

Accepted Answer

James Tursa
James Tursa on 13 Dec 2018
Edited: James Tursa on 13 Dec 2018
When doing physics problems, I would advise that you always include units and description in a comment off to the side to make sure it is easy to spot mismatches. E.g.,
rho=1.225; % (kg/m^3), density of atmosphere near Earth surface
cd=0.479; % (dimensionless), coefficient of drag
m=5.5; % (kg), mass
D=0.1; % (m), diameter? (you fill this in ...)
R=D/2; % (m), radius? (you fill this in ...)
A=4*pi*R^4; % (m^4), cross-section area? (you fill this in ...) WRONG?
g=9.81; % (m/s^s), acceleration due to gravity at Earth surface
and
rx(1)=0; % (m) initial x position
ry(1)=h; % (m) initial y position
v(1)=v0; % (m/s), initial velocity magnitude
Just looking at it, the A term, which I assume is Area, is very suspicious with that R^4 term giving what looks like m^4 for the units. I would expect this to be m^2 for the units. So maybe this should be an R^2 in your formula (only you know the shape of your object and what the correct formula for the cross-section area of your object is ... but the units must come out at m^2 to be valid for the downstream code). If it's a sphere then maybe you just need pi*R^2 on the rhs.

More Answers (2)

madhan ravi
madhan ravi on 13 Dec 2018

the cyclist
the cyclist on 13 Dec 2018
If I increase cd to a very large value, for example
cd=10000;
then I see the impact of drag -- the trajectory is no longer a parabola.
My guess is that you have a units problem, and some parameters is off by orders of magnitude.

Categories

Find more on Polar Plots in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!