Simulating projectile trajectories with (simple) air resistance.
Update: I originally wrote this script after learning that a lot of undergraduate students - as I was, at the time - thought that, because over the three years the B.Sc. lasts - they had always neglected air resistance, so it must be negligible, right? Of course it wasn't - it just indicated a model that was more difficult to solve. This is a simulation of a very simple air resistance on a sphere; you see how much the trajectories are *not* parabolas.
These are simply the calculated trajectories when you include simple air resistance.
The simple air resistance is just 1/2 rho v^2 A c_d. Everything in there is estimated,
with the exception of the velocity v, which is calculated.
You can also see why, at one point, projectile trajectories were thought to be triangular. In common day terms, it was thought that the energy given made an artificial movement, after which the energy would slowly be used, then the natural movement would kick in and it'd fall straight down. Yes, that's what they once thought.
Matlab script is included.
These are simply the calculated trajectories when you include simple air resistance.
The simple air resistance is just 1/2 rho v^2 A c_d. Everything in there is estimated,
with the exception of the velocity v, which is calculated.
You can also see why, at one point, projectile trajectories were thought to be triangular. In common day terms, it was thought that the energy given made an artificial movement, after which the energy would slowly be used, then the natural movement would kick in and it'd fall straight down. Yes, that's what they once thought.
Matlab script is included.
projectilemotion.m
close all;
format compact;
clear;
cd = 0.43; %Dimensionless
rho = 1.225; % kg/m/m/m
vStart = 340;
vAngle = pi/4;
m = 0.340e-1; %kg
g = 9.81; % m/s/s/kg
%[x,y] = projectilefunction(cd,rho,vStart,vAngle,m,g);
%plot(x,y);
%xlabel('Horizontal distance, metres');
%ylabel('Vertical distance, metres');
vLin = linspace(1,340,25);
aLin = linspace(0, pi/2, 25);
vNum = 1;
aNum = 1;
s = 1;
%Attempt to save gif
[xpos,ypos, conv] = projectilefunction(cd,rho,vStart,vAngle,m,g);
plot(xpos,ypos);
title(sprintf('Trajectory for vStart=%2.2e, angle=%2.2e pi', vStart, vAngle/pi));
xlim([0 25]);
ylim([0 25])
xlabel('Horizontal Distance (m)');
ylabel('Vertical Distance (m)');
figure(1);
frame = getframe(figure(1));
[im,map] = rgb2ind(frame.cdata,256,'nodither');
im(1,1,1,25*25) = 0;
%end attempt
for vStart = vLin
for vAngle = aLin
[xpos,ypos, vs, conv] = projectilefunction(cd,rho,vStart,vAngle,m,g);
fprintf('Max distance %2.2e, max height %2.2e \n', max(xpos), max(ypos));
plot(xpos,ypos);
title(sprintf('Trajectory for vStart=%2.2e, angle=%2.2e pi', vStart, vAngle/pi));
%legend(sprintf('Trajectory for vStart=%2.2e, angle=%2.2e pi', vStart, vAngle/pi));
xlim([0 1]);
ylim([0 1])
xlabel('Horizontal Distance (m)');
ylabel('Speed size (m/s)');
Frames(s) = getframe;
s = s + 1;
end
end
movie(Frames, 1);
function [ x,y, vs, conv ] = projectilefunction( cd,rho,vStart,vAngle,m,g )
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
%cd = 0.43; %Dimensionless
%rho = 1.225; % kg/m/m/m
%vStart = 10;
%vAngle = pi/4;
%m = 1; %kg
%g = 9.81; % m/s/s/kg
position = zeros(2,1e3);
dt = 1e-3; %s
i = 1;
v = [vStart*cos(vAngle); vStart*sin(vAngle)];
gravity = [0; -m*g];
grounded = 0;
vs = zeros(1,1e3);
while(grounded == 0)
drag = -0.5*rho*v'*v*cd*v/sqrt(v'*v);
if( v(1,1) == 0 && v(2,1) == 0)
break
end
v = v + (gravity + drag)*dt/m;
if(i > 1)
position(1,i) = position(1,i-1) + v(1,1)*dt;
position(2,i) = position(2,i-1) + v(2,1)*dt;
else
position(1,i) = v(1,1)*dt;
position(2,i) = v(2,1)*dt;
end
if position(2,i) <= 0
grounded = 1;
else
i = i + 1;
end
end
x = position(1, 1:i);
y = position(2, 1:i);
conv = i;
vs = v'*v;
end
