function dw = mfs_vlmdwh(panels, controls, kred, motion, cref)

# usage: dw = mfs_vlmdwh(panels, controls, kred, motion, cref)
#
# Input  panels    Structure with panel data
#        controls  Structure with control surface data
#        kred      Reduced frequency
#        motion    Structure with definition of motion
#        cref      Reference chord length
#
# Output dw(:)     Downwash matrix
#
# The function computes the downwash matrix for the harmonic 
# vortex-lattice method for one reuced frequency.
#
# ---------------------------------------------------------------------

# Check arguments

  if (nargin != 5 || nargout != 1)
     print_usage();
  end

# Initialize

  nofpan = panels.nofpan;
  refpnt = motion.refpnt;

  if ((ncntrl = controls.ncntrl))
     fldnam = fieldnames(controls);
     cntnam = fldnam(2 : end);
  else
     cntnam = [];
  end

  dw = complex(zeros(nofpan, 1));

# Process motion definitions

  for [val, key] = motion

      switch key

      case "heave"
         dw += (2 * i * kred * val / cref) * ones(nofpan, 1);

      case "pitch"
         v   = (panels.C(1, :)' - refpnt(1)) .* panels.nvec(3, :)';
         dw -= val * (1 + 2 * i * kred * v / cref);

      case "roll"
         val = 2 * i * kred * val / cref;
         v   = (panels.C(2, :)' - refpnt(2)) .* panels.nvec(3, :)';
         v  -= (panels.C(3, :)' - refpnt(3)) .* panels.nvec(2, :)';
         dw += val * v;

      case cntnam
         ctl = controls.(key);
         for l = 1 : ctl.nofls
             pids = ctl.pids(l, 1) : ctl.pids(l, 2);
             eyL = [panels.nvec(3, pids(1)), -panels.nvec(2, pids(1))];
             dy  = eyL * (ctl.P2(l, 2 : 3) - ctl.P1(l, 2 : 3))';
             s   = (eyL / dy) * (panels.C(2 : 3, pids) - ctl.P1(l, 2 : 3)');
             xLE  = (1 - s) * ctl.P1(l, 1) + s * ctl.P2(l, 1);
             valf = val * ctl.factors(l);
             dw(pids) -= valf * (1 + 2 * i * kred * ...
                         (panels.C(1, pids) - xLE)' / cref);
         end

      otherwise

      end

  end

end
