MATLAB: Sample Code of Projectile Motion. Unknown constant that’s used.

projectile motion

Disclamer: I am using the following code to help guide me through a projectile motion problem.
function projectile(theta,v0,dt,tf)
t = 0:dt:tf;
vx = v0*cos(theta);
vy = v0*sin(theta);
vx(1) = vx;
vy(1) = vy;
v = sqrt((vx^2)+(vy^2));
v(1) = v;
x(1) = 0;
y(1) = 0;
B2divm = 0.00004;
g = 9.81;
for i = 1:length(t)-1
    v = sqrt((vx(i)^2)+(vy(i)^2));
    x(i+1) = x(i) + vx(i)*dt;
    vx(i+1) = vx(i) - (B2divm*v*vx(i))*dt;
    y(i+1) = y(i) + vy(i)*dt;
    vy(i+1) = vy(i) - g*dt - (B2divm*v*vy(i))*dt;
end
plot(x,y,'r')
end
I follow and reasonably understand the code, except for the use of the constant 'B2divm=0.0004'.
I strongly suspect that 'B2divm' above has something to do with Drag and air density, but I just don't know.
Anybody know what it means? How it's derived?
Thanks!!

Best Answer

You are probably right about that.
As you can see, the position states (X & Y) are updated by adding a delta position which is computed as V*dt, i.e.
x(i+1) = x(i) - vx(i)*dt;
y(i+1) = y(i) - vy(i)*dt;
You would expect a similar expression for the velocity update. 
vx(i+1) = vx(i) + ax(i)*dt
vy(i+1) = vy(i) + ay(i)*dt
The delta velocity is the acceleration times the time step,
deltaV = a*dt. So, a reaonable expression for the acceleratio of a projectile is the sum of (weight + drag force) divided by the mass. Note that the constant name is B2divm - maybe "B2 divided by mass". This suggests that "B2" is the sum of the forces acting on the projectile. I see that there is a g*dt in the Y axis, therefore "B2" represents only the aerodynamic drag force.
Drag Force = 0.5*(air density)*Velocity^2*(Reference Aera)*(Drag Coefficient)
Note that the Drag Force is computed using the total velocity (a scalar value). We want to determine the amount of the force in the X and Y directions, so the X component is (Drag Force) * cos(gamma) and
the Y component is (Drag Force)*sin(gamma) where gamma is the flight path angle (i.e. direction of the velocity vector).
We substitude cos(gamma) = Vx/V and sin(gamma) = Vy/V, now we have:
Fx = (Drag Force) * Vx/V
Fy = (Drag Force) * Vy/V
or
Fx = 0.5*(air density) * V^2 * (Reference Area) * (Drag Coefficient) * Vx / V
= 0.5*(air density) * (Reference area) * (Drag Coefficient) * V * Vx
Now the acceleration in the X direction due to air resistance = Fx / m, so
Ax(from drag) = {0.5 * (air density) * (Reference Area) * (drag Coefficient) / m } * V * Vx
For short trajectories, everything inside the curly brackets can be considered as constant, and this would be the "B2 divided by m" term, B2divm.
If the trajectory has any significant variation in altitude, the air density term must be computed inside the for loop for each pass.
Related Question