function [C, CPcoor, CPnvec] = mfs_vlmmat(ls, panels, symy, ys)

# usage: [C, CPcoor, CPnvec] = mfs_vlmmat(ls, panels, symy, ys)
#
# Input  ls           Structure with lifting surface data
#        panels       Structure with panel data
#        symy         Flag indicating symmetry with respect to plane
#                     y = ys
#        ys           y-coordinate of symmetry plane
#
# Output C(:, :)      Matrix C
#        CPcoor(:, 3) Control point coordinates
#        CPnvec(:, 2) Control point normals (ny and nz)
#
# The function computes the aerodynamic C-matrix of the vortex-lattice
# method. Rows correspond to control points and columns to vortices.
#
# ---------------------------------------------------------------------

# Copyright(c) 2020 by Johannes Wandinger

# Check arguments

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

# Initialize

  nofls  = ls.nofls;
  nofpan = panels.nofpan;
  fourpi = 4 * pi;
  C      = zeros(nofpan, nofpan);  
  rC     = panels.C';

# Symmetry

  if (symy)
     nloop = 2;
  else
     nloop = 1;
  end

# Collect control point coordinates and normal vectors

  CPcoor = panels.C'; CPnvec = panels.nvec(2 : 3, :)';

# Loop over the lifting surfaces
# ------------------------------

  for l = 1 : nofls

      nx = ls.surfs(l).nx; ny = ls.surfs(l).ny;
      m1 = ls.surfs(l).p1; m2 = ls.surfs(l).pend;

      iv  = m1 : m2;                % Vortices
      ixA = 1 : nx * ny;            % Points A
      ixB = nx + 1 : ixA(end) + nx; % Points B

      for loop = 1 : nloop  % Symmetry loop

# Collect vortex point coordinates

          if (loop == 1)
             VPcoor = [panels.A(:, iv), panels.B(:, m2 - nx + 1 : m2)];
             f = fourpi;
          else
             f = -fourpi;
             VPcoor(2, :) = ys - VPcoor(2, :);
          end

# Compute vectors between points

          xAB = VPcoor(1, ixB) - VPcoor(1, ixA);
          yAB = VPcoor(2, ixB) - VPcoor(2, ixA);
          zAB = VPcoor(3, ixB) - VPcoor(3, ixA);

          xVP = CPcoor(:, 1) - VPcoor(1, :);
          yVP = CPcoor(:, 2) - VPcoor(2, :);
          zVP = CPcoor(:, 3) - VPcoor(3, :);

# Compute cross product of AP x BP

          vx = yVP(:, ixA) .* zVP(:, ixB) - zVP(:, ixA) .* yVP(:, ixB);
          vy = zVP(:, ixA) .* xVP(:, ixB) - xVP(:, ixA) .* zVP(:, ixB);
          vz = xVP(:, ixA) .* yVP(:, ixB) - yVP(:, ixA) .* xVP(:, ixB);

# Compute cross product ex x VP

          wy = -zVP; wz = yVP;

# Normalize vectors VP

          nVP = sqrt(xVP .* xVP + yVP .* yVP + zVP .* zVP);
          xVP = xVP ./ nVP; yVP = yVP ./ nVP; zVP = zVP ./ nVP;

# Compute contribution of segment AB

          nv = vx .* vx + vy .* vy + vz .* vz;

          tiny = eps * nVP(:, ixA) .* nVP(:, ixB);
          nv(nv < tiny) = 1;

          nv = f * nv;

          rx = xVP(:, ixA) - xVP(:, ixB);
          ry = yVP(:, ixA) - yVP(:, ixB);
          rz = zVP(:, ixA) - zVP(:, ixB);

          s  = (xAB .* rx + yAB .* ry + zAB .* rz) ./ nv;

          C(:, iv) += CPnvec(:, 1) .* (vy .* s);
          C(:, iv) += CPnvec(:, 2) .* (vz .* s);

# Compute contribution of semi-infinite segments

          nw = f * (wy .* wy + wz .* wz);
          s  = (1 + xVP) ./ nw;

          C(:, iv) -= CPnvec(:, 1) .* (wy(:, ixA) .* s(:, ixA));
          C(:, iv) -= CPnvec(:, 2) .* (wz(:, ixA) .* s(:, ixA));
          C(:, iv) += CPnvec(:, 1) .* (wy(:, ixB) .* s(:, ixB));
          C(:, iv) += CPnvec(:, 2) .* (wz(:, ixB) .* s(:, ixB));

      end % Symmetry loop

  end

end
