function [cmp, rc] = mfs_rigtrim(cmp)

# usage: [cmp, rc] = mfs_rigtrim(cmp)
#
# Input  cmp     Structure with component data
#
# Output cmp     Structure with component data
#                (Trim results added)
#        rc      Return code: 0 = no errors
#                             1 = errors
#
# The function performs a trim analysis of a rigid aircraft. The
# values of the trim parameters are written to the output file.
#
# --------------------------------------------------------------------

# Check arguments

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

# Trim parameter names (without control surfaces)
# Columns of downwash matrix must be in this order!

  tpnames = {"ax", "ay", "az", "racce", "pacce", "yacce", ...
             "alpha", "beta", "pitch", "yaw", "roll"};

# Initialize

  rc     = 0;
  tp     = struct();
  nconf  = cmp.config.nconf;

# Check existence of required data

  if (! nconf)
     printf("*E* mfs_trim: No configurations defined\n");
     rc = 1;
  end

  if (! isfield(cmp, "rbmp"))
     printf("*E* mfs_trim: No rigid body mass properties defined\n");
     rc = 1;
  end

  if (rc) return; end

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

  if (! ncntrl)
     printf("*E* mfs_trim: No control surfaces defined\n");
     rc = 1; return;
  end

# Mark trim parameters that are angles

  isangle = [0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0];
  isangle = [isangle, ones(1, ncntrl)];

# Preset trim parameters

  tpall = [tpnames, fieldnames(cmp.config.cntls)'];
  ntp   = length(tpall);
  ntpa  = ntp - 6;

  if (sum(cmp.config.nlinc))
     tpr = [tpall, "rhs"];
  end

# Identify the trim variables that are defined

  tpdef = zeros(ntp, nconf);
  n     = 1;
  for name = tpall
      val = getfield(cmp.config.tpdef, name{1});
      tpdef(n++, :) = bitget(val, 1 : nconf);
  end

# Check the number of trim variables

  for n = 1 : nconf
      ntpl = ntp - sum(tpdef(:, n)) - cmp.config.nlinc(n);
      if (ntpl != 6)
         printf("*E* Configuration %d:\n", n);
         printf("    Number of trim parameters to be determined ");
         printf("must be 6 but is %d\n", ntpl);
         rc = 1;
      end
  end

  if (rc) return; end

# Build the rigid body mass matrix

  MRR = zeros(6, 6);
  MRR(1 : 3, 1 : 3) = cmp.rbmp.m * eye(3, 3);
  MRR(4 : 6, 4 : 6) = cmp.rbmp.JS;

# Build the load integrator matrix to compute load resultants from
# vortex strengths

  X = zeros(6, nofpan);

  rG  = 0.5 * (cmp.panels.A + cmp.panels.B) - cmp.rbmp.cm';
  rAB = cmp.panels.B - cmp.panels.A;
  fn  = 2 * [-rAB(3, :); rAB(2, :)];  % [Fy, Fz]

  X(2, :) = fn(1, :);
  X(3, :) = fn(2, :);
  X(4, :) = -rG(3, :) .* fn(1, :) + rG(2, :) .* fn(2, :);
  X(5, :) = -rG(1, :) .* fn(2, :);
  X(6, :) = rG(1, :) .* fn(1, :); 

# Build the downwash matrix

  dw = mfs_vlmdw(cmp.panels, cmp.controls, cmp.config, cmp.rbmp.cm);

# Build the matrix of influence coefficients

  C = mfs_vlmmat(cmp.ls, cmp.panels, 0, 0);

# Compute matrix to compute load resultants from configuration
# parameters (first column is contribution of camber)

  G = C \ dw;
  Q = X * G;

# Perform trim analysis of each configuration

  tpval = zeros(ntp, nconf);

  n = 1;
  for name = tpnames
      tpval(n++, :) = getfield(cmp.config, name{1});
  end
  for [val, key] = cmp.config.cntls
      tpval(n++, :) = val;
  end

  for n = 1 : nconf

      nlc  = cmp.config.nlinc(n);
      ixtl = find(tpdef(:, n) == 0);
      ixtp = find(tpdef(:, n));
      qdyn = cmp.config.qdyn(n);

      T    = [MRR, -qdyn * Q(:, 2 : ntpa + 1)];
      rhs  = qdyn * Q(:, 1);

      if (nlc)
         C      = zeros(nlc, ntp + 1);
         lincon = cmp.config.lincon{n};
         for l = 1 : nlc
             lc = lincon{l};
             for [val, key] = lc
                 ix = find(strcmp(key, tpr));
                 C(l, ix) = val;
             end
         end
         T   = [T; C(:, 1 : ntp)];
         rhs = [rhs; C(:, end)];
      end

      Tl  = T(:, ixtl);
      Tp  = T(:, ixtp);
      rhs = rhs - Tp * tpval(ixtp, n);

      tpval(ixtl, n) = Tl \ rhs;

  end

# Vortex strengths

  tpvala = tpval(7 : ntp, :);
  Gtrim  = G * [ones(1, nconf); tpvala];

# Store the results

  cmp.trim = struct("type", "rigid", "nconf", nconf, "ntp", ntp, 
                    "G", Gtrim);

  cmp.trim.tpnames = tpall;
  cmp.trim.tpval   = tpval;
  cmp.trim.isangle = isangle;

end
