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 

params = [rhog, a, mu, g, mp];

% Initial conditions
u0 = 0.0;

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

v(1) = u0;

% For first step, use forward Euler
v(2) = v(1) + dt*drop_rhs(v(1),params);

% Start iterative loop using midpoint method
for n = 3:N+1,

% Update using midpoint method
  v(n) = v(n-2) + 2*dt*drop_rhs(v(n-1),params);

end

% Plot results
t = linspace(0,Tmax,N+1);
subplot(311);
plot(t,v);
ylabel('v (m/sec)');
title('Midpoint 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');