function [C, S, wake] = mfs_vlmhpre(lsf, panels, pancols, symy, ys, cref, opts)

# usage: [C, S, wake] = mfs_vlmhpre(lsf, panels, pancols, symy, ys, cref, opts)
#
# Input  lsf      Structure with lifting surface data
#        panels   Structure with panel data
#        pancols  Structure with panel column data
#        symy     Flag indicating symmetry with respect to plane y = ys
#        ys       y-coordinate of symmetry plane
#        cref     Reference chord length
#        opts     Structure with options
#
# Output C(:, :)  Matrix of influence coefficients of steady vortex
#                 lattice method
#        S(:, :)  Summation matrix
#        wake(:)  Structure array with wake information
#                 .nrows             Number of rows in matrix F
#                 .ncols             Number of columns in matrix F
#                 .rows(2)           Indices of first and last row in
#                                    matrix F
#                 .cols(2)           Indices of first and last column in
#                                    matrix F
#                 .s(nx)             Integration grid
#                 .FS(nrows, ntep)   Factors of integral IS
#                 .FF1(nrows, ncols) Factors of integral IF1
#                 .FF2(nrows, ncols) Factors of integral IF2
#                 .Y(:, nx)          Function values of IS, IF1 and IF2
#                 .ixY(:)            (1) Index of last entry of IS
#                                    (2) Index of last entry of IF1
#
# Fields of structure opts:
#
#    msg    File handle for messages (if zero, no messages will be issued)
#    nx     Number of intervals per chord length cref
#    m1     Factor to get size of first interval from mimimum control
#           point distance
#    lenw   Length of wake behind last control point in multiples of cref
#
# The function computes the matrix C of the steady vortex lattice method,
# the summation matrix S and the functions to be integrated along the wake.
#
# -------------------------------------------------------------------------

# Copyright(c) 2024 by Johannes Wandinger

# Check arguments

  if (nargin != 7 || nargout != 3)
     print_usage();
  endif

# Initialize

  wake = struct();
  msg  = opts.msg;

  if (msg)
     fprintf(msg, "Information on wakes\n\n");
  endif

  fourpi = 4 * pi * [1, -1];

  nofpan = panels.nofpan;   % Number of panels or control points
  npcol  = pancols.npcol;   % Number of panel columns

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

# Compute matrix C of steady VLM

  [C, P, nvec] = mfs_vlmmat(lsf, panels, symy, ys);

# Build summation matrix

  S = sparse(npcol, nofpan);
  for k = 1 : npcol
      pids = pancols.panids{k};
      S(k, pids) = 1;
  endfor

# Loop over panel column groups and lifting surfaces

  nws = 1;   % counts the wake segments

  for k = 1 : pancols.ng % Panel colum groups

      % Panel column indices of current wake segment

      col1  = pancols.gcx(1, k); 
      col2  = pancols.gcx(2, k);
      ncols = col2 - col1 + 1;

      % Number of trailing edge points Q (D or E)

      ntep  = ncols + 1;

      % Coordinates of trailing edge points Q (D or E)
      % (ntep columns)

      Q  = squeeze(pancols.tepcor(:, 1, col1 : col2));
      Q  = [Q, pancols.tepcor(:, 2, col2)];

      for loop = 1 : nloop  % Symmetry loop

          if (msg)
             if (symy)
                if (loop == 1)
                   fprintf(msg, "  Trailing edge %4.0f (original)\n", k);
                else
                   fprintf(msg, "  Trailing edge %4.0f (mirrored)\n", k);
                endif
             else
                fprintf(msg, "  Trailing edge %4.0f\n", k);
             endif
          endif

          fps = fourpi(loop);

          if (loop == 2)
             Q(2, :) = ys - Q(2, :);
          endif

          % Vectors DE (ncols columns)

          DE = Q(:, 2 : end) - Q(:, 1 : end-1);  % Vectors DE

          for l = 1 : lsf.nofls % Lifting surfaces

              if (msg)
                 fprintf(msg, "    Lifting surface %4.0f: ", lsf.ids(l));
              endif

              if (loop == 1)
                 te = lsf.surfs(l).gid == k;
              else
                 te = 0;
              endif

#           Wake column indices in matrix F

              wake(nws).cols  = [col1, col2]; 
              wake(nws).ncols = ncols;

#           Control point indices

              m1 = lsf.surfs(l).p1; m2 = lsf.surfs(l).pend;
              nrows = m2 - m1 + 1;
              wake(nws).rows = [m1, m2]; wake(nws).nrows = nrows;

#           Number of functions

              nsemi = ntep  * nrows; % semi-infinite vortex lines
              nfini = ncols * nrows; % finite vortex lines

#           Control point distances from trailing edge points Q

              XQP = P(m1 : m2, 1) - Q(1, :);   % rows = control points
              YQP = P(m1 : m2, 2) - Q(2, :);   % cols = vortex lines
              ZQP = P(m1 : m2, 3) - Q(3, :);

#           Integration grid

              [s, npx] = wakegrid(XQP, lsf.surfs(l).nx, te, cref, opts);
              cs   = cref * s;

#           Function Y1: IS

              % ex x rQP

              C1y = -ZQP; C1z = YQP;

              % Component normal to lifting surface

              C1n = nvec(m1 : m2, 1) .* C1y + nvec(m1 : m2, 2) .* C1z;

              % Factors FS

              wake(nws).FS = C1n ./ (fps * (C1y .* C1y + C1z .* C1z));

              % Coordinates in one column

              xQP = reshape(P(m1 : m2, 1) - Q(1, :), nsemi, 1);
              yQP = reshape(P(m1 : m2, 2) - Q(2, :), nsemi, 1);
              zQP = reshape(P(m1 : m2, 3) - Q(3, :), nsemi, 1);

              % x-coordinates of vortex points V (A or B)

              xVP = xQP - cs;

              % Norms of rVP

              nQP = 1 ./ sqrt(xVP .* xVP + yQP .* yQP + zQP .* zQP);

              % Wake function values

              Y1  = 1 + xVP .* nQP;

              clear xVP;

#           Functions Y2 and Y3: IF1 and IF2

              % rDP x rEP

              C2x = YQP(:, 1 : ncols) .* ZQP(:, 2 : ntep) ...
                  - ZQP(:, 1 : ncols) .* YQP(:, 2 : ntep);
              C2y = ZQP(:, 1 : ncols) .* XQP(:, 2 : ntep) ...
                  - XQP(:, 1 : ncols) .* ZQP(:, 2 : ntep);
              C2z = XQP(:, 1 : ncols) .* YQP(:, 2 : ntep) ...
                  - YQP(:, 1 : ncols) .* XQP(:, 2 : ntep);

          % ex x rDE

              C3y = -DE(3, :); C3z = DE(2, :);  

          % Factors FF1 and FF2

              wake(nws).FF1 = (C2y .* nvec(m1 : m2, 1) ...
                             + C2z .* nvec(m1 : m2, 2)) / fps;
              wake(nws).FF2 = (C3y / fps) .* nvec(m1 : m2, 1) ...
                            + (C3z / fps) .* nvec(m1 : m2, 2);

          % Common denominator of both integrals

              a0 = reshape(C2x .* C2x + C2y .* C2y + C2z .* C2z, nfini, 1);
              a1 = reshape(C2y .* C3y + C2z .* C3z, nfini, 1);
              a2 = repmat(C3y .* C3y + C3z .* C3z, nrows, 1);
              a2 = reshape(a2, nfini, 1);

              a  = a0 + 2 * a1 * cs + a2 * (cs .* cs);

              clear C2x; clear C2y; clear C2z; clear C3y; clear C3z;

              tiny = eps * max(abs(a0));
              a(a < tiny) = 1;
              a = 1 ./ a;

              clear a0; clear a1; clear a2;

              % Scalar product of DE with DP

              cx = 1 : ncols;
              B1 = DE(1, :) .* XQP(:, cx) ...
                 + DE(2, :) .* YQP(:, cx) ...
                 + DE(3, :) .* ZQP(:, cx);
              cx = 2 : ntep;
              b1 = reshape(B1, nfini, 1);
              clear B1;

              % Scalar product of DE with EP

              B2 = DE(1, :) .* XQP(:, cx) ...
                 + DE(2, :) .* YQP(:, cx) ...
                 + DE(3, :) .* ZQP(:, cx);
              b2 = reshape(B2, nfini, 1);
              clear B2;

              % DE multiplied by cref * s

              bs = DE(1, :)' .* cs;

              % Expanded bs (each row repeated nrows time)

              bsx = zeros(nfini, npx);          
              ir1 = 1; ir2 = nrows;
              for m = 1 : ncols
                  bsx(ir1 : ir2, :) = repmat(bs(m, :), nrows, 1);
                  ir1 += nrows; ir2 += nrows;
              endfor

              clear bs;

              % Common nominator function b

              iD = 1 : nfini; iE = nrows + 1 : nsemi;
              b  = (b1 - bsx) .* nQP(iD, :);
              b -= (b2 - bsx) .* nQP(iE, :);

              % Functions Y2 and Y3

              Y2 = b .* a;

              Y3 = (cs .* b) .* a;

              clear a; clear b;

              % Store values in structure wake

              wake(nws).s   = s;
              wake(nws).ixY = [nsemi, nsemi + nfini];
              wake(nws).Y   = [Y1; Y2; Y3];

              nws++;

          endfor  % End of lifting surface loop

      endfor % End of symmetry loop

  endfor % End of panel column group loop

endfunction

function [s, ns] = wakegrid(XQP, npr, te, cref, opts)

# usage: [s, ns] = wakegrid(XQP, npr, te, cref, opts)
#
# Input  XQP(nc, nv)   Control point distances from trailing
#                      edge points
#        npr           Number of panel rows (in x-direction)
#        te            is 1 if current lifting surface is in front of
#                      trailing edge processed, otherwise 0
#        cref          Reference chord length
#        opts          Structure with wake options
#
# Output s(ns)         Wake grid
#
# Wake options: 
#
#    msg    File handle for messages (if zero, no messages will be issued)
#    nx     Number of intervals per chord length cref
#    m1     Factor to get size of first interval from mimimum control
#           point distance
#    lenw   Length of wake behind last control point in multiples of cref
#
# The function generates the wake grid s. The x-coordinates of the
# edge points of the horseshoe vortices of the wake are given by
# x = s * cref.
#
# -------------------------------------------------------------------------

# Initialize

  msg      = opts.msg;
  [nP, nQ] = size(XQP);

  if (te)

# Lifting surface is in front of trailing edge

     % Data of regular grid

     lenw = opts.lenw;
     nlin = lenw * opts.nx + 1; % number of intervals

     % Near field grid

     dx1  = 0.5 * (XQP(npr, 1) + XQP(npr, 2));
     dx2  = 0.5 * (XQP(nP, nQ) + XQP(nP, nQ - 1));
     xmin = min(abs([dx1, dx2]));
     ds1  = xmin / (opts.m1 *cref);
     dsN  = 1 / opts.nx;

     if (ds1 < dsN)
        d   = dsN / ds1;
        r   = (1 - ds1) / (1 - dsN);
        nxn = 1 + ceil(log(d) / log(r));
        dsn = zeros(1, nxn); dsn(1) = ds1;
        for k = 2 : nxn
            dsn(k) = r * dsn(k-1);
        endfor
        s = [0, cumsum(dsn)];
        nlin -= opts.nx;
        s1 = s(end);          % starting value of regular grid
        ns1 = nxn + 1;        % index of first point of regular grid
        ns2 = nxn + nlin;     % index of last point of regular grid
     else
        s1   = 0;             % starting value of regular grid
        ns1  = 1; ns2 = nlin; % index of first and last point
     endif

  else

# Lifting surface not in front of trailing edge processed

     % Adjust wake length

     xmax = max(max(XQP([npr, nP], [1, nQ]))); 
     lenw = ceil(xmax/cref) + opts.lenw;

     % Data of regular grid

     nlin = lenw * opts.nx + 1; % number of intervals
     s1   = 0;                  % starting value
     ns1  = 1; ns2 = nlin;      % index of first and last point

  endif

# Regular grid

  s(ns1 : ns2) = linspace(s1, lenw, nlin);
  ns = length(s);

# Messages

  if (msg)
     fprintf(msg, "nx   = %4.0f, te = %1.0f\n", ns, te);
     if (te)
        fprintf(msg, "                          ");
        fprintf(msg, "ds1  = %10.3e, dsN = %10.3e", ds1, dsN);
        if (ds1 < dsN)
           fprintf(msg, ", nxn  = %4.0f\n", nxn);
        else
           fprintf(msg, "\n");
        endif
     else
        fprintf(msg, "                          ");
        fprintf(msg, "xmax = %10.3e, lenw = %4.0f c\n", xmax, lenw);
     endif
  endif

endfunction
