function dw = mfs_vlmdw(panels, controls, config, refpnt)

# usage: dw = mfs_vlmdw(panels, controls, config, refpnt)
#
# Input  panels    Structure with panel data
#        controls  Structure with control surface data
#        config    Structure with configuration data
#        refpnt(3) Coordinates of reference point for rotations
#
# Output dw(:, :)  Downwash matrix: rows correspond to panels
#
# The function computes the downwash matrix for the vortex-lattice
# method. The columns correspond to trim variables:
#
#  Col. 1:  Downwash due to camber and angle of incidence
#       2:  Angle of attack
#       3:  Side slip angle
#       4:  Pitch rate
#       5:  Yaw rate
#       6:  Roll rate
#      >6:  Control surfaces
#
# If refpnt is undefined, the columns for pitch, yaw and control rate
# are not built.
#
# ---------------------------------------------------------------------

# Check arguments

  if (nargin < 3 || nargin > 4 || nargout != 1)
     print_usage();
  end

# Initialize

  nofpan = panels.nofpan;
  ncntrl = config.ncntrl;

  if (nargin == 4)
     ncol = 6 + ncntrl;
  else
     ncol = 3 + ncntrl;
  end 

  nvec = panels.nvec;

# Process trim paramters

  dw = zeros(nofpan, ncol);

  dw(:, 1) = panels.slope';   % Camber
  dw(:, 2) = -nvec(3, :)';    % Angle of Attack
  dw(:, 3) = -nvec(2, :)';    % Sideslip angle

  if (nargin == 4)
     coorC = panels.C - refpnt';
     dw(:, 4) = (coorC(3, :) .* nvec(1, :))';    % Pitch
     dw(:, 4) -= (coorC(1, :) .* nvec(3, :))';
     dw(:, 5) = (coorC(1, :) .* nvec(2, :))';    % Yaw
     dw(:, 5) -= (coorC(2, :) .* nvec(1, :))';
     dw(:, 6) = (coorC(2, :) .* nvec(3, :))';    % Roll
     dw(:, 6) -= (coorC(3, :) .* nvec(2, :))';
     n = 7;
  else
     n = 4;
  end

# Process control surfaces

  if (config.ncntrl)
     for [angle, cntnam] = config.cntls
         cntl    = getfield(controls, cntnam);
         pids    = cntl.pids;
         factors = cntl.factors;
         for s = 1 : cntl.nofls
             p  = pids(s, 1) : pids(s, 2);
             np = length(p);
             dw(p, n) = -repmat(factors(s), np, 1);
         end
         n++;
     end
  end

end
