function [cmps, cmpa, rc] = mfs_flextrim(cmp, meth)

# usage: [cmps, cmpa, rc] = mfs_flextrim(cmp, meth)
#
# Input  cmp      Structure with component data
#        meth     Method: 1 = restrained
#                         2 = unrestrained
#
# Output cmps     Structure with data of solid component
#        cmpa     Structure wiht data of aerodynamic component
#        rc       Return code: 0 = no errors
#                              1 = errors
#
# The function performs a trim analysis of a flexible aircraft.
#
# --------------------------------------------------------------------

# Check arguments

  if (nargin != 2 || nargout != 3)
     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"};

# Trim types

  types = {"restrained", "unrestrained"};

# Initialize

  rc     = 0;
  cmps   = struct();
  cmpa   = struct();
  nconf  = cmp.aero.config.nconf;

# Check existence of required data

  if (! isfield(cmp.solid, "stiff"))
     printf("*E* mfs_trim: Stiffness matrix missing\n");
     rc = 1;
  end
  if (! isfield(cmp.solid, "mass"))
     printf("*E* mfs_trim: Mass matrix missing\n");
     rc = 1;
  end

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

  if (! isfield(cmp, "splines"))
     printf("*E* mfs_trim: Splines missing\n");
     rc = 1;
  end

  if (rc) return; end;

  ncntrl = cmp.aero.config.ncntrl;
  nofpan = cmp.aero.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)];

# Initialize Output

  cmps = cmp.solid;
  cmpa = cmp.aero;

# Preset trim parameters

  tpall = [tpnames, fieldnames(cmp.aero.config.cntls)'];
  zz    = zeros(1, nconf);
  ntp   = length(tpall);
  ntpa  = ntp - 6; % Number of aerodynamic trim parameters

  if (sum(cmp.aero.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.aero.config.tpdef, name{1});
      tpdef(n++, :) = bitget(val, 1 : nconf);
  end

# Check the number of trim variables that are defined

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

  if (rc) return; end

# List of local and dependent dofs

  dofl = cmp.solid.dofs.dofl;

  if ((ndofd = cmp.solid.dofs.ndofd))
     dofd = cmp.solid.dofs.dofd;
  end

# Build the rigid body modes and the rigid body mass matrix

  R = mfs_rigidmotion(cmp.solid.nodes, cmp.solid.dofs.mxdofpnt,
                      [0, 0, 0]);
  Mrr = R' * cmp.solid.mass.M * R;
  m   = Mrr(1);
  cm  = [Mrr(2, 6), Mrr(3, 4), Mrr(1, 5)] / m;

  R = mfs_rigidmotion(cmp.solid.nodes, cmp.solid.dofs.mxdofpnt, cm);

  Mr  = cmp.solid.mass.M * R;
  Mrr = R' * Mr;
  Mlr = mfs_matpartr(Mr, cmp.solid.dofs);

# Build the downwash matrix

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

# Build the matrix G to compute forces from vortices

  G = 2 * diag(cmp.aero.panels.width);

# Build the matrix of influence coefficients

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

# Compute matrices to get vortices from downwash

  Vxmat = C \ [dw, cmp.splines.D1h];  % [G0, GK, Gh]

# Compute building blocks of aerodynamic matrices

  Sh  = cmp.splines.Svh' * G;
  Vxh = Sh * Vxmat;
  q0h = Vxh(:, 1);
  qhK = Vxh(:, 2 : ntpa + 1);
  Qhh = Vxh(:, ntpa + 2 : end);

# Spline matrices to local and rigid body dofs

  Shg  = cmp.splines.Shg;
  Shl  = mfs_matpartc(Shg, cmp.solid.dofs);
  Shr  = Shg * R;

# Initialize displacements

  ndofg = cmp.solid.dofs.ndofg;
  disp  = zeros(ndofg, nconf);

  if (meth == 2)
     qr = zeros(6, nconf);
  end

# Perform trim analysis of each configuration

  tpval = zeros(ntp, nconf);

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

  [nh, nk] = size(qhK);
  nr       = columns(Mlr);

  RHS0  = [Mlr, Shl' * [qhK, q0h, Qhh]];
  Kll   = mfs_matpart(cmp.solid.stiff.K, cmp.solid.dofs, [1, 0]);
  RHS   = Kll \ RHS0;
  ixM   = 1 : nr;
  ixqlK = (nr + 1) : (nr + nk);
  ixql1 = (nr + nk + 1) : (nr + nk + 1);
  ixqlh = (nr + nk + 2) : (nr + nk + 1 + nh);
  Ihh   = eye(nh);

  for n = 1 : nconf

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

      U0lt  = [qdyn * RHS(:, ixql1), RHS(:, ixM), -qdyn * RHS(:, ixqlK)];
      Rhh   = Ihh - (qdyn * Shl) * RHS(:, ixqlh); 
      U0lt += RHS(:, ixqlh) * (Rhh \ (Shl * (qdyn * U0lt)));

      TF1 = [Mrr, -qdyn * Shr' * qhK];
      TF  = TF1 + qdyn * Shr' * Qhh * (Shl * U0lt(:, 2 : end));

      rhs0 = qdyn * Shr' * (q0h + Qhh * Shl * U0lt(:, 1));

      if (nlc)
         C      = zeros(nlc, ntp + 1);
         lincon = cmp.aero.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
         TF   = [TF; C(:, 1 : ntp)];
         rhs0 = [rhs0; C(:, end)];
      end

      TFl  = TF(:, ixtl);
      TFp  = TF(:, ixtp);
      rhs  = rhs0 - TFp * tpval(ixtp, n);
      tpval(ixtl, n) = TFl \ rhs;

      disp(dofl, n)  = U0lt(:, 1);
      disp(dofl, n) -= U0lt(:, 2 : end) * tpval(:, n);
      if (ndofd)
         disp(dofd, n) = cmp.solid.dofs.C(dofd, :) * disp(:, n);
      end

      if (meth == 2)

         qr(:, n) = Mrr \ (Mlr' * disp(dofl, n));
         disp(:, n) -= R * qr(:, n);

         rhs0 = qdyn * Shr' * (q0h + Qhh * Shg * disp(:, n));

         if (nlc)
            TF1  = [TF1; C(:, 1 : ntp)];
            rhs0 = [rhs0; C(:, end)];
         end

         Tl  = TF1(:, ixtl);
         Tp  = TF1(:, ixtp);
         rhs  = rhs0 - Tp * tpval(ixtp, n);
         tpval(ixtl, n) = Tl \ rhs;
 
      end

  end

# Vortex strengths

  tpvala = tpval(7 : ntp, :);
  Gamma  = Vxmat(:, 1) + Vxmat(:, 2 : ntpa + 1) * tpvala;
  Gamma += Vxmat(:, ntpa + 2 : end) * (Shg * disp);

# Store the results

  cmpa.trim = struct("type", types{meth}, "nconf", nconf, 
                     "ntp", ntp, "G", Gamma);

  cmpa.trim.tpnames = tpall;
  cmpa.trim.tpval   = tpval;
  cmpa.trim.isangle = isangle;
  cmpa.trim.disp    = cmp.splines.Snh * (Shg * disp);

  cmps.trim = struct("nconf", nconf, "disp", disp);

  if (meth == 2)
      cmpa.trim.qr = qr;
      cmps.trim.qr = qr;
  end

end
