% 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
%In this example, the analytical solution solution, a numerical 
%solution using Euler's method and the error between the two
%solutions are calculated and then plotted.
%
clc
clear

g=9.8   %m/s^2   gravitational constant
m=70     %kg      mass of falling sphere 65
c=12.5  %kg/s    drag coeff. of sphere
dt=2   %s       time step
nt=24   %        number of time steps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  calculate analytical solution v1  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

maxtime=(nt*dt)                                 %duration of experiment (seconds)
t1=0:dt:maxtime;                                %generate time array
v_analytic =(g*m/c)*(1-exp(-t1*c/m));           %calculate velocity analytically
v_terminal(1:length(v_analytic))=m*g/c;             %calculate terminal velocity

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%generate time array for Euler method %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% implement Euler's method for falling sphere
veuler(1)=0.                                            %initialize velocity
for i=2:nt+1,
    veuler(i)=veuler(i-1)+(g-veuler(i-1)*c/m)*(t2(i)-t2(i-1));
end

% calculate absolute error - difference between the analytical solution
% abd Euler Solution
for i=1:nt+1
    error(i) = abs(v_analytic(i)-veuler(i));
end


% generate plot and labels
plot(t1,v_terminal,'-',t1,v_analytic,t2,veuler,'*',t2,error*10)
title('Sphere speed as a function of Time')
xlabel('Time (s)')
ylabel('V (m/s)')
key2 = ['Euler Method dt =', num2str(dt)];          %convert dt into a character string
                                                    %and concatonate with text to be
                                                    %used in legend
legend('Terminal Velocity','Analytic Solution',key2,'Absolute Error*10',0)                
