function [cmp, nfreq, rc] = mfs_freqrespc(cmp, f, opts, args)

# usage: [cmp, nfreq, rc] = mfs_freqrespc(cmp, f, opts, args)
#
# Input  cmp        Structure with aeroelastic component
#        f(:)       List of excitation frequencies
#        opts       Structure with wake options
#        args       Parameters
#
# Output cmp        Structure with aeroelastic compnent
#        nfreq      Number of frequencies
#        rc         Return code: 0 = no errors
#                                1 = errors
#
# Parameters:
#
#   Name       Values    Default   Meaning
#  -----------------------------------------------------------------------
#   method     fsm         fsm     Force summation method
#              modal               Modal reduction
#              direct              Direct frequency response
#   loadcase    int        1       Loadcase number
#   nband       int        0       Number of additional frequencies in
#                                  halfpower bandwidth 
#
# ------------------------------------------------------------------------

# Copyright (c) 2018 by Johannes Wandinger

# Check arguments

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

# Initialize output

  nfreq = 0;
  rc    = 0;

# Default values of parameters

  meth =  4;
  nb   =  0;
  lc   =  1;

# Process parameters

  for k = 1 : 2 : length(args)
      switch args{k}
      case "method"
         switch args{k + 1};
         case "fsm"
            meth = 4;
         case "modal"
            meth = 2;
         case "direct"
            meth = 1;
         otherwise
            printf("*E* mfs_freqresp: unknown method %s\n", args{k + 1})
            rc = 1;
         end 
      case "loadcase"
         lc = args{k + 1};
         if (lc < 1)
            printf("*E* mfs_freqresp: loadcase must be positive\n");
            rc = 1;
         end
      case "nband"
         nb = args{k + 1};
         if (nb < 0)
            printf("*E* mfs_freqresp: nband must be positive\n");
            rc = 1;
         end
      otherwise
         printf("*E* mfs_freqresp: unknown parameter %s", args{k})
         rc = 1;
      end
  end

  if (rc) return; end

# Make sure excitation frequencies are ascending

  f = sort(f);

# Compute response

  if (meth == 1)
     [resps, respa, rc] = mfs_dfreq(cmp, f, lc, opts);
     if (rc) return; end
     nfreq = length(f);
  else
     [resps, respa, nfreq, rc] = mfs_mfreq(cmp, f, nb, lc, meth, opts);
     if (rc); return; end
     if (! isfield(cmp.aero, "modes"))
        Xa = cmp.splines.Snh * (cmp.splines.Shg * cmp.solid.modes.disp);
        modes = struct("nofmod", cmp.solid.modes.nofmod,
                       "rigmod", cmp.solid.stiff.ndofr,
                       "omega",  cmp.solid.modes.omega,
                       "freq",   cmp.solid.modes.freq,
                       "disp",   Xa);
        cmp.aero.modes = modes;
     end
  end

# End

  if (! isfield(cmp.solid, "freqresp"))
     cmp.solid.freqresp.nfa = lc;
     cmp.aero.freqresp.nfa  = lc;
  else
     cmp.solid.freqresp.nfa = max(lc, cmp.solid.freqresp.nfa);
     cmp.aero.freqresp.nfa  = max(lc, cmp.solid.freqresp.nfa);
  end
  cmp.solid.freqresp.resp{lc} = resps;
  cmp.aero.freqresp.resp{lc}  = respa;

end
