PROGRAM: DRAGF.M
function [drf, trp,md]=dragf(vx,H,yc)
% This m-file computes the Drag force component
% in dimensionless form up to a value of y<=yc
% here
% yc is the (not dimensionless) crest elevation
% vx is the matrix of x-velocities, the columns are
%       phases and the rows are distances measured
%       vertically from the bottom
% H is the wave Height
%
% The numerical integration is performed using the trapezodial
%  rule
[ny,nth]=size(vx);
 trp1=tril(ones(ny,ny));
 trp=2*trp1-diag(ones(ny,1))-[[ones(ny,1)] zeros(ny,ny-1) ];
drf=1/2*trp*(vx.*abs(vx)/H);
dy=yc/200;
drf=dy*drf;
%now get the drag moment
md=1/2*trp*(o9nes(201,1)*[0:200]*(vx.*abs(vx))/(H*H));