function [ustar_cr,theta_cr,d_star]=shields(d50,s)
% d50 in microns
%function [ustar_cr,theta_cr,d_star]=shields(d50,s)

% calculate shields curve mobility parameter
if nargin==1, s=2.65; end
nu=0.011; % cm**2/s
g=981.; % cm/s**2
d50=d50*1e-4; % microns ==> cm
d_star=d50*((s-1)*g/nu^2)^(1/3);

theta_cr=nan*ones(size(d_star));
ustar_cr=nan*ones(size(d_star));

ind=find(d_star <= 4.);
if(~isempty(length(ind))),
  theta_cr(ind)=.24*d_star(ind).^(-1);
end
ind=find(4 <= d_star & d_star <= 10.);
if(~isempty(length(ind))),
  theta_cr(ind)=.14*d_star(ind).^(-.64);
end
ind=find(10 <= d_star & d_star <= 20.);
if(~isempty(length(ind))),
  theta_cr(ind)=.04*d_star(ind).^(-.10);
end
ind=find(20 <= d_star & d_star <= 150.);
if(~isempty(length(ind))),
  theta_cr(ind)=.013*d_star(ind).^(.29);
end
ind=find(150 <= d_star );
if(~isempty(length(ind))),
  theta_cr(ind)=.055;
end
ustar_cr=(theta_cr.*d50*(s-1)*g).^(0.5);
