function [q_b,phi_b,theta_p]=mpm(d50,ustar_c,ustar,s);
%function [q_b,phi_b,theta_p]=mpm(d50,ustar_c,ustar,[s]);
%  Meyer Peter Muller formula
%    d50 in cm
%    ustar_c in cm/s
%    ustar in cm/s
%
%   q_b dimensional bedload (cm^2/s)
if(nargin==3),s=2.65;end;

g=981;    % cm/s^2
d=d50;    % cm

phi_b=8* ((ustar.^2-ustar_c^2)/((s-1)*g*d)).^(1.5);
q_b=phi_b*d*sqrt((s-1)*g*d);    % cm^2/s
