r = ripples( H, T, h, D, rhow, rhos )  


g = 9.81;        % gravitational accel (m/s)
w = 2*pi ./T;    % wave angular frequency (radians/s) 

kh = qkhf(w, h)
do = H/(sinh(kh)
A = 2*d0         % wave orbital amplitude = 2*d0
ub = A.*w;       % wave velocity amplitude A*w (m/s)

s = (rhow/rhos)
psi = (ub).^2 ./(s-1)*g*D