% This function roughly models the entrance region velocity profile by
% simply linearizing the profile shape over the entrance region length from
% Vav to Va*(1-r.^2./R^2)

function V = Vfluid (r, R, z, Va, L_e)

%r: position in cylindrical coord
%R: Radius of the tube
%z: position in the tube
%Va: Vmax in the velocity profile
%L_e: Entrance Region Length

if z <= 0
    V=Va.*sin(pi/3).*ones(size(r));
elseif z > L_e
    V=Va .* (1-r.^2./R^2);
else
    Vav=Va./2;
    q=(L_e-z)./L_e;
    V=(Va-Vav.*q).*(1-r.^2./R^2) + Vav.*q ;
end