MATLAB: Plotting projectile with drag

dragplotprojectile

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

Best Answer

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.
Related Question