function varargout = mfs_freqresp(cmp, fk, varargin)

# usage: varargout = mfs_freqresp(cmp, fk, varargin)
#
#  [cmp, f, ES, WR] = mfs_freqresp(cmp, f, varargin)
#  cmp              = mfs_freqresp(cmp, kred, opts)
#  cmp              = mfs_freqresp(cmp, f, varargin)
#  cmp              = mfs_freqresp(cmp, f, opts, varargin)
#
# The first format is for solid components, the second for aerodynamic
# components and the third and fourth for aeroelastic components.
#
# Input  cmp        Structure with component
#        f(:)       List of excitation frequencies
#        kred(:)    List of reduced frequencies
#        varargin   Parameters
#        opts       Structure with wake options (optional)
#
# Output cmp        Structure with component
#                   (Frequency response added)
#        f(:)       List of excitation frequencies (optional)
#        ES(:)      Strain energies (optional)
#        WR(:)      Work of residual load (optional)
#
# Parameters in solid mechanics:
#
#   Name       Values    Default   Meaning
#  -----------------------------------------------------------------------
#   method     enhanced  enhanced  Enhanced modal reduction
#              fms       fsm       Force summation (Aeroelasticity)
#              modal               Modal reduction
#              direct              Direct frequency response
#   nband       int        5       Excitation frequencies per
#                                  halfpower bandwidth
#   loadcase    int        1       Loadcase number
#
# Wake options:
#
#    Name    Default     Meaning
#  ------------------------------------------------------------------
#    msg        0        File handle of message file
#                        (if 0, no messages are printed)
#    nx        10        Number of intervals per chord length cref
#    m1         4        Factor to divide mimimum control point distance
#                        from trailing edge to get size of first interval
#    lenw      20        Length of wake behind last control point in
#                        multiples of cref
#    tmp        0        Indicates whether a temporary file should be
#                        used (1) or not (0)
#    ktol     1e-4       Tolerance for reduced frequencies
#
# ------------------------------------------------------------------------

# Copyright (c) 2020 by Johannes Wandinger

  t0 = clock();

# Default options for aerodynamics

  optsaero = struct("msg", 0, "nx", 10, "m1", 4, "lenw", 20,
                    "tmp", 0, "ktol", 1e-4);

# Check arguments

  if (nargin < 2 || nargout < 1)
     print_usage();
  end
  if (! isstruct(cmp))
     error("mfs_freqresp: %s must be a structure\n", inputname(1));
  end
  if (ismatrix(fk))
     if (! ismember(1, size(fk)))
        error("mfs_freqresp: %s must be 1-dimensional\n", inputname(2));
     end
     fk = sort(fk(:)');
  else
     error("mfs_freqresp: %s must be an array\n", inputname(2));
  end

# Branch on component type

  switch cmp.type

  case "solid"
     if (nargout > 4)
        print_usage();
     end
     time_stamp = ignore_function_time_stamp("all");
     mfs_paths("add", "mfs_freqresp.m", "solid", "util");
     [cmp, f, ES, WR, rc] = mfs_freqrespx(cmp, fk, varargin);
     nfreq = length(f);
     varargout{1} = cmp;
     varargout{2} = f;
     varargout{3} = ES;
     varargout{4} = WR;
     mfs_paths("remove");
     ignore_function_time_stamp(time_stamp);

  case "aero"
     if (nargout != 1)
        print_usage();
     end
     if (nargin == 2)
        opts = optsaero;
     else
        opts = mfs_defopts(varargin{1}, optsaero);
     end
     time_stamp = ignore_function_time_stamp("all");
     mfs_paths("add", "mfs_freqresp.m", "aero", "util");
     [cmp, rc] = mfs_freqrespx(cmp, fk, opts);
     nfreq = length(fk);
     varargout{1} = cmp;
     mfs_paths("remove");
     ignore_function_time_stamp(time_stamp);

  case "aeroelastic"
     if (nargout < 1 || nargout != 1)
        print_usage();
     end
     if (nargin > 2 && isstruct(varargin{1}))
        opts = mfs_defopts(varargin{1}, optsaero);
        args = varargin(2 : end);
     else
        opts = optsaero;
        args = varargin;
     end
     time_stamp = ignore_function_time_stamp("all");
     mfs_paths("add", "mfs_freqresp.m", "aeroelastic", "aero", "util");
     [cmp, nfreq, rc] = mfs_freqrespc(cmp, fk, opts, args);
     varargout{1} = cmp;
     mfs_paths("remove");
     ignore_function_time_stamp(time_stamp);

  otherwise
     error("mfs_freqresp: model type \"%s\" not supported\n",
           cmp.type);

  end

  if (rc)
     error("mfs_freqresp ended with errors\n");
  end

  elapsed_time = etime(clock(), t0);
  printf("%10.4f seconds needed to process %5d frequencies of component %s\n", ...
         elapsed_time, nfreq, inputname(1));
  fflush(stdout);

end
