% DRAG EXAMPLE - Newton's 2nd Law
%ME 2040 - Matlab script to calculate the velocity of a falling
%sphere subject to aerodynamic forces as a function of time.
%The differential equation that is modeled here is:
% dv/dt=g-(c/m)v
%   v is the velocity of the sphere
%   t is time
%   c drag coef.
%   g gravitational accelleration cost.
%   m mass of sphere
%Two equivalent methods for handling arrays are shown.

g=9.81 %m/s^2   gravitational constant
m=68.1 %kg      mass of falling sphere
c=12.5 %kg/s    drag coeff. of sphere
dt=0.1 %s       time step
nt=60  %        number of time steps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  method without for loop     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
maxtime=(nt*dt)                     %duration of experiment (seconds)
t1=0:dt:maxtime;                    %generate time array - t1
v1=(g*m/c)*(1-exp(-t1*c/m));        %calculate velocity - v1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  method using for loop where array %
%  elements are explicitly defined   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t2(1)=0;                               %initialize time
for i=2:nt+1
    t2(i)=t2(i-1)+dt;                  %time array
end

for i=1:nt+1,
v(i)=(g*m/c)*(1-exp(-t2(i)*c/m));                    
end

% generate plot and label
plot(t2,v,t1,v1,'*')
title('Sphere speed as a function of Time')
xlabel('Time (s)')
ylabel('V (m/s)')

% write data to an ASCII file
fid=fopen('sphere_vel.dat','w');
fprintf(fid,'%7.3f  \t %7.3f \n',t1,v1);            %c-style formatting
fclose(fid);