clear all;
close all;

% Initialize some parameters
rhog = 0.9; % Density of air at 3000m
a    = 0.01; % Particle radius
mu   = 1.69E-5; % Viscosity of air at 3000m
g    = 9.8; % Gravity
mp   = 917*4*pi/3*a^3; % mass of ice particle 

% Initial conditions
u0 = 0.0;

% Set time length of integration, and number of steps
Tmax = 25;
N = 100;
dt = Tmax/N;

v(1) = u0;
% Start iterative loop
for n = 2:N+1,

% Calculate drag at n-1
  u = v(n-1);
  Re = rhog*u*(2*a)/mu;
  if (abs(Re) > 0),
    CD = 24/Re + 6/(1+sqrt(Re)) + 0.4;
    D = 0.5*rhog*u^2*pi*a^2*CD;
  else
    D = 0;
  end

% Calculate right-hand sides at n-1
  f = g - D/mp;

% Update using Forward Euler
  v(n) = v(n-1) + dt*f;

end

% Plot results
t = linspace(0,Tmax,N+1);
subplot(311);
plot(t,v);
ylabel('v (m/sec)');
title('Forward Euler method');

subplot(312);
Re = rhog*v*(2*a)/mu;
plot(t,Re);
ylabel('Re');

subplot(313);
CD = 24./Re + 6./(1+sqrt(Re)) + 0.4;
plot(t,CD);
xlabel('time (sec)');ylabel('C_D');