clear
global alpha beta gamma delta g

alpha=0;
beta=1;
gamma=2.5e-07;
delta=0.0124;
g=9.81;

r0=0.0001;
v0=0;

figure(1)
tspan=[0 120];

% y = [v, r]
y0=[v0;r0];
[t,y]=ode45(@raindropf,tspan,y0);
plot(1000*y(:,2),y(:,1),'k')
r=y(:,2);
set(gca,'XTick',[.1 .1001 .1002 .1003])

figure(2)
subplot(1,2,1)
plot(t,1000*y(:,2),'k')
set(gca,'YTick',[.1 .1001 .1002 .1003])
ylabel('r (mm)')
xlabel('t (s)')
subplot(1,2,2)
plot(t,y(:,1),'k')
xlabel('t (s)')
ylabel('v (m/s)')

figure(1)
% Solution for alpha=0, beta=0
v=g/4/gamma*r.*(1-(r0./r).^4);
hold on
%plot(1000*r,v)


% Solution for alpha=0, beta=1
v=sqrt(2*g/7/gamma*r).*sqrt(1-(r0./r).^7);

%plot(1000*r,v,'r')
hold off

xlabel('r (mm)')
ylabel('v (m/s)')
%axis([0 1 0 110])