% This is some code that demonstrates linear and quadratic drag in a
% variety of scenarios.
%
%% Startup Block
set(0,'DefaultAxesFontSize',20);
set(0,'DefaultLineLineWidth',3);
%% Linear Drag, Horizontal Part
t=[0:.01:10];
m=1;
b=0.3;
vox=10;
vold=vox*ones(size(t));
v=vox*exp(-b.*t./m);
figure
subplot(2,1,1)
plot(t,vold,'k')
hold on
plot(t,v,'r')
xlabel('Time')
ylabel('Horizontal Velocity')
legend('Classical','Linear Air Resistance','Location','East')
axis([min(t) max(t) 0 1.1*vox])
grid on
xold=vold.*t;
x=(m*vox./b)*(1-exp(-b.*t./m));
subplot(2,1,2)
plot(t,xold,'k')
hold on
plot(t,x,'r')
xlabel('Time')
ylabel('Horizontal Position')
legend('Classical','Linear Air Resistance','Location','NorthWest')
axis([min(t) max(t) 0 1.1*max(xold)])
plot(t,m*vox./b*ones(size(t)),'r:')
grid on
%% Linear Drag, Vertical Part
t=[0:.01:10];
m=1;
b=0.3;
voy=-30;
g=9.81;
vold=voy+g.*t;
v=exp(-b.*t./m).*(voy-m*g./b)+m*g./b;
figure
subplot(2,1,1)
plot(t,-vold,'k')
hold on
plot(t,-v,'r')
xlabel('Time')
ylabel('Vertical Velocity')
legend('Classical','Linear Air Resistance','Location','NorthEast')
grid on
yold=voy.*t+(1/2)*g.*t.*t;
y=(m/b)*((1-exp(-b.*t./m))*(voy-m*g/b)+g.*t);
subplot(2,1,2)
plot(t,-yold,'k')
hold on
plot(t,-y,'r')
xlabel('Time')
ylabel('Vertical Position')
legend('Classical','Linear Air Resistance','Location','SouthWest')
grid on
%% Linear Trajectory
x=[0:.1:200];
vox=20;
voy=30;
g=9.81;
m=1;
b=0.35;
yold=voy.*x./vox-(g/2)*(x.*x./(vox*vox));
y=((voy+m*g/b)./vox).*x+((m*g/b)*(m/b)*log(1-b.*x/(vox*m)));
newmax=min(find(y<0));
newmax2=min(find(yold<0));
figure
plot(x(1:newmax2-1),yold(1:newmax2-1),'k')
hold on
plot(x(1:newmax-1),y(1:newmax-1),'r')
xlabel('Horizontal Displacement')
ylabel('Vertical Displacement')
legend('Classical','Linear Air Resistance','Location','NorthWest')
grid on
%% Quadratic -- Horizontal Motion
vox=10;
m=1;
c=0.2;
t=[0:.01:5];
vold=vox*ones(size(t));
v=(m*vox./(m+c.*t.*vox));
figure
subplot(2,1,1)
plot(t,vold,'k')
hold on
plot(t,v,'r')
xlabel('Time')
ylabel('Horizontal Velocity')
title('Horizontal-Only Motion')
legend('Classical','Quadratic Air Resistance','Location','East')
axis([min(t) max(t) 0 1.1*vox])
grid on
xold=vold.*t;
x=(vold.*m./(c.*vox)).*log(1+t./(m./(c*vox)));
subplot(2,1,2)
plot(t,xold,'k')
hold on
plot(t,x,'r')
xlabel('Time')
ylabel('Horizontal Position')
legend('Classical','Quadratic Air Resistance','Location','NorthWest')
axis([min(t) max(t) 0 1.1*max(xold)])
grid on
%% Quadratic -- Vertical Motion (free fall)
voy=0;
m=1;
c=0.05;
t=[0:.01:8];
g=9.81;
vold=voy+g.*t;
v=sqrt(m*g/c)*tanh(g.*t./(sqrt(m*g/c)));
figure
subplot(2,1,1)
plot(t,-vold,'k')
hold on
plot(t,-v,'r')
xlabel('Time')
ylabel('Vertical Velocity')
title('Vertical-Only Motion')
legend('Classical','Quadratic Air Resistance','Location','SouthWest')
axis([min(t) max(t) -40 voy*1.1])
grid on
yold=voy.*t+(1/2)*g.*t.*t;
y=(m./c)*log(cosh(g.*t./(sqrt(m.*g./c))));
subplot(2,1,2)
plot(t,-yold,'k')
hold on
plot(t,-y,'r')
xlabel('Time')
ylabel('Vertical Position')
legend('Classical','Quadratic Air Resistance','Location','SouthWest')
axis([0 max(t) -80 1])
grid on