function res = mfs_modecont(varargin)

# usage: res = mfs_modecont(varargin)
#
#  res = mfs_modecont(fid, cmp, rid, freq, ldc, thresh)
#  res = mfs_modecont(cmp, rid, freq, ldc, thresh)
#
# Input  fid          File handle (optional)
#        cmp          Structure with component
#        rid(:, 2)    Response identifiers
#        freq(:)      List of frequencies
#        ldc(:)       List of load cases (optional)
#        thresh(2)    Threshold (1) for printed output
#                               (2) for structure polar
#
#
# Output res(:, :)    Structure array with results (optional)
#                     (columns correspond to load cases)
#
# Fields of structure res:
#   rid(2)            Nodal point and degree of freedom identifier
#   rfreq(nofmod)     Eigenfrequencies
#   freq(:)           Frequencies
#   resp(:)           Total response
#   mcf(nofmod, :)    Modal contribution factors
#   polar(nfreq)      Structure array with data to create polar plots
#
# Fields of structure polar:
#   title             Title containing load case number, nodal point,
#                     degree of freedom identifier and frequency
#   theta(2, :)       Phase angles of contribution factors
#   rho(2, :)         Magnitudes of contribution factors
#   legend{:}         Legend
#
# --------------------------------------------------------------------

# Copyright (c) 2023 by Johannes Wandinger

  t0 = clock();

# Initialize

  res = struct();
  prt = 0; nxa = 1;
  rc  = 0;
  rd2dg = 180 / pi;

  thresh = [0.01, 0.1];

# Checks
# ------

# Arguments

  nargs = length(varargin);
  if (nargs < 3 || nargs > 6 || nargout > 1)
     print_usage();
  end

  out = nargout > 0;

  if (is_valid_file_id(varargin{1}))
     fid = varargin{1};
     prt = 1;
     nxa = 2;
  end

  cmp = varargin{nxa};
  if (! isstruct(cmp))
     error("mfs_modecont: argument %s is not a component\n",
           inputname(nxa));
  end
  ixcmp = nxa++;

  rid = varargin{nxa++};
  if (columns(rid) != 2)
     error("mfs_modecont: wrong definition of response identifiers %s\n",
           inputname(nxa-1));
  end

  freq = varargin{nxa++};

  if (nxa <= nargs)
     ldc    = varargin{nxa};
     nofldc = length(ldc);
  else
     ldc    = 1;
     nofldc = 1;
  endif

  if (++nxa == nargs)
     th  = varargin{nxa};
     nth = length(th);
     if (nth == 2)
        thresh = th;
     elseif (nt == 1)
        thresh = [th, th]
     else
        error("mfs_modecont: wrong dimensions of argument %s\n",
              inputname(nxa));
     end
  end

# Component type

  if (! strcmp(cmp.type, "solid"))
     error("mfs_modecont: component type must be \"solid\"\n");
  end

# Existence of frequency reponse results

  if (! isfield(cmp, "freqresp"))
     error("mfs_modecont: no frequency response data found\n");
  end

# Loadcases

  if (min(ldc) < 1 || max(ldc) > cmp.freqresp.nfa)
     error("mfs_modecont: loadcase %d or loadcase %d does not exist\n",
           min(ldc), max(ldc));
  end

# Computations
# ------------
  
  mfs_paths("add", "mfs_modecont.m", "solid", "util");

# Get degree of freedom indices

  [ixdof, rc, noddof] = mfs_finddof(cmp.nodes.ids, cmp.dofs.mxdofpnt,
                                    rid(:, 1), rid(:, 2), cmp.nset);
  if (rc)
     mfs_paths("remove");
     error("mfs_modecont aborted due to previous errors\n");
  end

  nofresp = length(ixdof);
  noddof(:, 1) = cmp.nodes.ids(noddof(:, 1));

  modes  = cmp.modes.disp;
  rfreq  = cmp.modes.freq;
  ndofr  = cmp.stiff.ndofr;

  if (out)
     res(nofresp, nofldc) = struct();
     [res.rfreq] = deal(rfreq);
  end

  if (prt)
     fprintf(fid, "\n");
     for k = 1 : 8
         fprintf(fid, "----------");
     end
     fprintf(fid, "\n\nModal contribution factors of component \"%s\"\n\n",
             inputname(ixcmp));
  end

# Loop over loadcases

  for l = 1 : nofldc

      lc   = ldc(l);
      resp = cmp.freqresp.resp{lc};
      rfq  = rfreq;

      if (resp.method == 1)
         printf("*E* mfs_modecont: modal contribution factors not ");
         printf("available\n")
         printf("                  for direct frequency response ");
         printf("analysis\n");
         printf("                  load case %d skipped\n", lc);
         rc = 1;
         continue
      end

      if (prt)
         fprintf(fid, "Loadcase %2d:\n\n", lc);
      end

#     Get frequencies

      f     = sort(freq);
      nfreq = length(f);

      ftest = [resp.freq, 10 * resp.freq(end)];
      ixl   = lookup(ftest, f); ixu = ixl + 1;
      df    = [f - ftest(ixl); ftest(ixu) - f];
      [m, im] = min(df);
      ix  = [ixl; ixu];
      ixf = zeros(1, nfreq);
      for k = 1 : nfreq
          ixf(k) = ix(im(k), k);
      end
      f = ftest(ixf);

#     Compute contribution factors

      dspa = mfs_getdispf(resp, modes, ixdof, ndofr, 1);
      dspf = dspa(:, ixf);

      Q = resp.Q(:, ixf);

      for n = 1 : nofresp
          mcf{n} = modes(ixdof(n), :)' .* Q;
      end

      if (resp.method >= 3 || resp.type > 1)
         for n = 1 : nofresp
             sc     = dspf(n, :) - sum(mcf{n});
             mcf{n} = [mcf{n}; sc];
         end
         rfq = [rfq; Inf];
      end

#     Print contribution factors

      if (prt)

         nofmod = length(rfq);

         for n = 1 : nofresp

             fprintf(fid, "  Node %6d, Dof %1d:\n\n", noddof(n, :));

             mcfac = mcf{n};

             for k = 1 : nfreq

                 fprintf(fid, "    Frequency = %7.2f Hz: Resp. = (%10.3e, %10.3e)\n\n",
                         f(k), real(dspf(n, k)), imag(dspf(n, k)));
                 fprintf(fid, "                               absolute     ");
                 fprintf(fid, "         relative  \n");
                 fprintf(fid, "      mode   frequency     real       imag");
                 fprintf(fid, "        real     phase\n");

                 relmcf = mcfac(:, k) / dspf(n, k);
                 [rs, ix] = sort(abs(real(relmcf)), "descend");
                 mout = length(find(rs > thresh(1)));

                 for m = 1 : mout
                     mm = ix(m);
                     fprintf(fid, "     %4d  %8.2f Hz  %10.3e %10.3e",
                             mm, rfq(mm), real(mcfac(mm, k)), imag(mcfac(mm, k)));
                     if (abs(relmcf(mm) < 10))
                        fprintf(fid, "    %6.3f  %7.2f\n",
                                real(relmcf(mm)), rd2dg * arg(relmcf(mm)));
                     else
                        fprintf(fid, "     ****      ****\n");
                     end
                 end

                 fprintf(fid, "\n");

             end

         end

      end

#     Add results to output structure

      if (out)

         [res(:, l).freq] = deal(f);

         for n = 1 : nofresp

             res(n, l).rid  = noddof(n, :);
             res(n, l).resp = dspf(n, :);
             res(n, l).mcf  = mcf{n};

             for k = 1 : nfreq

                 p.title = sprintf("LC %d Node %d Dof %d: f = %7.2f Hz",
                                   l, noddof(n, :), f(k));

                 mcfac = mcf{n}(:, k);
                 rcfac = mcfac / dspf(n, k);
                 [s, ix] = sort(abs(real(rcfac)), "descend");
                 mout = length(find(s > thresh(2)));
                 ix   = ix(1 : mout++);

                 theta = zeros(2, mout); rho = theta;
                 theta(2, 1)       = arg(dspf(n, k));
                 theta(2, 2 : end) = arg(mcfac(ix));
                 rho(2, 1)         = abs(dspf(n, k));
                 rho(2, 2 : end)   = abs(mcfac(ix));

                 leg = cell(mout, 1);
                 leg{1} = "Total";
                 for m = 2 : mout
                     modno = ix(m - 1);
                     if (rfq(modno) < Inf)
                        leg{m} = sprintf("Mode %d", modno);
                     else
                        leg{m} = "Static";
                     end
                 end

                 p.theta = theta; p.rho = rho; p.legend = leg;

                 res(n, l).polar(k) = p;

             end

         end

      end

  end

  mfs_paths("remove");

  if (rc)
     error("mfs_modecont aborted due to previous errors\n");
  end

  elapsed_time = etime(clock(), t0);
  printf("%10.4f seconds needed to compute modal contributions of component %s\n", ...
         elapsed_time, inputname(2));
  fflush(stdout);

end
