% Calculate the droplet velocity including films with the Brinkman model
% Created by Mathias Nagel - EPFL - LFMI November 2013 lfmi.epfl.ch/simcomics
function U  = dropvel(  Ca, lambda, rh )
% > U is the ratio of droplet velocity over carrier fluid mean velocity
% > Ca is the capillary number Ca=mu V/gamma, with the carrier fluid mean
% velocity
% > lambda is the viscosity ratio, lambda=mu_drop/mu_carrier, it will be
% invalid for lambda>>1
% > rh is the geometrical aspect ratio, rh=R/H, radius over channel height

    k= rh*sqrt(12);
    I0 = besseli(0,k);
    I1 = besseli(1,k);
    K0 = besselk(0,k);
    K1 = besselk(1,k);
    H = 4*K1./K0 + 4*K1.*(I0-2*I1./k)./K0.*K1.*k./(2*K1.*I1-K1.*I0.*k...
        -lambda.*I1.*K0.*k+K0.*(lambda-1).*(2*I0-4*I1./k));

    A = (1.5245./((1+lambda).*k+H)).^3./Ca-3*(2*k+H)./((1+lambda).*k+H);
    B = 3*((2*k+H)./((1+lambda).*k+H)).^2;
    C = ((2*k+H)./((1+lambda).*k+H)).^3;

    U = 1/6*(-2*A + (2*2^(1/3)*(A.^2 - 3 *B))./(27*C - 2*A.^3+9*A.*B+...
        3*sqrt(3)*sqrt(27*C.^2- 4*A.^3.*C + 18*A.*B.*C - A.^2.*B.^2+ ...
        4*B.^3) ).^(1/3) +2^(2/3)*(27*C - 2*A.^3 + 9*A.*B +3*sqrt(3)*...
        sqrt(27*C.^2 - 4*A.^3.*C+18*A.*B.*C-A.^2.*B.^2+4* B.^3)).^(1/3) );
    U = real(U);    % imaginary part is round-off error close to zero
end