function rederr = mfs_reductionerror(fid, cmp, f, lcs)

# usage: rederr = mfs_reductionerror(fid, cmp, f, lcs)
#
# Input  fid          File handle
#        cmp          Structure with component
#        f(nf)        List of excitation frequencies (optional)
#        lcs(nlc)     List of load cases
#                     (optional, default is all loadcases)
# Output rederr(nlc)  Structure array with reduction errors
#                     (optional)
#
# Fields of structure array rederr:
#    f(nf)            Frequencies in acending order
#    abs(nf)          absolute error without static correction
#    rel(nf)          relativ error without static correction
#    eabs(nf)         absolute error with static correction
#    erel(nf)         relative error with static correction
#
# The function computes the modal strain energies and writes them to the 
# output file.
#
# If excitation frequencies are defined, then also the strain energy
# error is computed and upper bounds on the error are written to the 
# output file. If defined, the strain energy errors are returned in the
# output arguments.
#
# ------------------------------------------------------------------------

# Copyright (c) 2022 by Johannes Wandinger

  t0 = clock();

# Initialize

  rederr = struct();
  err  = 0;

# Check arguments

  if (nargin < 2 || nargin > 4 || nargout > 1)
     print_usage();
  end

  if (! is_valid_file_id(fid))
     error("mfs_reductionerror: first argument does not refer to an open file\n");
  end
  if (! isstruct(cmp))
     fclose(fid);
     error("mfs_reductionerror: second argument must be a structure\n");
  end
  isfreq = (nargin > 2) && (! isempty(f));
  if (nargout > 0 && ! isfreq)
     fclose(fid);
     error("mfs_reductionerror: frequencies needed to compute errors\n");
  end

# Check component type

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

# Check availability of normal modes

  if (! isfield(cmp, "modes"))
     fclose(fid);
     error("mfs_reductionerror: Normal modes not available\n");
  end

# Check availability of loads

  if (! isfield(cmp, "load"))
     fclose(fid);
     error("mfs_reductionerror: No loads defined\n");
  end

# Check loadcases

  nofldc = cmp.load.nofldc;
  if (nargin < 4)
     lcs = 1 : nofldc;
  else
     lcs = sort(lcs);
     if (lcs(end) > nofldc || lcs(1) < 1)
        fclose(fid);
        error("mfs_reductionerror: bad loadcase list\n");
     end
  end
  nlc = length(lcs);

# Get degree of freedom data

  ndofg = cmp.dofs.ndofg;
  ndofl = cmp.dofs.ndofl;
  dofl  = cmp.dofs.dofl;

# Get modal data needed

  nofmod = cmp.modes.nofmod;
  w      = cmp.modes.omega;
  fr     = cmp.modes.freq;
  fp     = fr(nofmod);
  X      = cmp.modes.disp(dofl, :);

# Check maximum frequency

  if (isfreq)
     f    = sort(f);
     fmax = f(end);
     err  = fmax > fp;
  end

# Get information on requested loadcases

  inflc  = cmp.load.inflc;
  infall = 0;
  for lc = lcs
      infall = bitor(infall, inflc(lc));
      chk = 0;
      for k = 1 : 4
          chk = chk + bitget(inflc(lc), k);
      end
      if (chk > 1)
         fclose(fid);
         error("mfs_reductionerror: loadcase %d has mixed types of excitation\n",
               lc);
      end
  end
  
  force  = bitget(infall, 1);
  enfmot = bitand(infall, 14);

  mfs_paths("add", "mfs_reductionerror.m", "util");

# Get stiffness matrix partitions

  if (enfmot)
     reqflags = [1, 1];
  else
     reqflags = [1, 0];
  end

  [Kll, Klp] = mfs_matpart(cmp.stiff.K, cmp.dofs, reqflags);

# Get data needed to process rigid body modes

  ndofr = cmp.stiff.ndofr;

  if (ndofr)
     dofe = cmp.stiff.dofe;
     Kee  = Kll(dofe, dofe);
     Xr   = X(:, 1 : ndofr);
     Mll  = mfs_matpart(cmp.mass.M, cmp.dofs, [1, 0]);
  end

# Get base motion

  if (enfmot)
     ndofp = cmp.dofs.ndofp;
     dofp  = cmp.dofs.dofp;
     up    = zeros(ndofp, nlc);
     ub    = zeros(ndofg, nlc);
     if (bitget(infall, 2))
        up = up + cmp.load.u(dofp, lcs);
     end
     if (bitget(infall, 3))
        up = up + cmp.load.v(dofp, lcs);
     end
     if (bitget(infall, 4))
        up = up + cmp.load.a(dofp, lcs);
     end
     if (ndofr)
        ubl = zeros(ndofl, nlc);
        ubl(dofe, :) = - Kee \ Klp(dofe, :) * up;
        XTM = Xr' * Mll;
        ubl += Xr * (XTM * ubl);
        ub(dofl, :) = ubl;
     else
        ub(dofl, :) = - Kll \ Klp * up;
     end
     ub(dofp, :) = up;
     if (cmp.dofs.ndofd)
        ub(cmp.dofs.dofd, :) = cmp.dofs.C(cmp.dofs.dofd, :) * ub;
     end
  end

# Compute static solution

  rhs = sparse(ndofl, nlc);
  if (force)
     rhs = mfs_matpartr(cmp.load.f(:, lcs), cmp.dofs);
  end
  if (enfmot)
     Mub = cmp.mass.M * ub;
     rhs = rhs + mfs_matpartr(Mub, cmp.dofs);
  end

  if (ndofr)
     qr   = Xr' * rhs;
     rhs  = rhs - Mll * (Xr * qr);
     us   = zeros(ndofl, nlc);
     us(dofe, :) = Kee \ rhs(dofe, :);
     n1 = ndofr + 1;
     X  = X(:, n1 : nofmod); w = w(n1 : nofmod); fr = fr(n1 : nofmod);
     nofmod -= ndofr;
  else
     us = Kll \ rhs;
  end

# Loop over loadcases

  fprintf(fid, "\n");
  for k = 1 : 8
      fprintf(fid, "----------");
  end
  fprintf(fid, "\n\nModal strain energies of component \"%s\"\n\n", ...
          inputname(2));

  for lc = 1 : nlc

     if (! inflc(lcs(lc))) continue; end

     fprintf(fid, "Loadcase %2d:\n\n", lcs(lc));

#    Compute strain energies

     Es  = us(:, lc)' * rhs(:, lc);
     Lx  = X' * rhs(:, lc) ./ w;
     EFm = Lx .^2 / Es;

#    Output strain energies

     fprintf(fid, "  mode   frequency      En/ES       Sum       1 - Sum\n\n");

     esum = 0;
     for m = 1 : nofmod
         esum  = esum + EFm(m);
         ediff = 1 - esum;
         fprintf(fid, "  %4d  %8.2f Hz  %8.5e  %8.6f  %8.5e\n", ...
                 m + ndofr, fr(m), EFm(m), esum, ediff);
     end

     fprintf(fid, "\n");

#    Compute strain energy errors

     if (isfreq && ! err)
        if (nargout)
           eta = f / fp;
        else
           eta = fmax / fp;
        end
        H  = 1 ./ (1 - eta.^2);
        rederr(lc).f    = f;
        rederr(lc).rel  = 0.5 * H.^2 * ediff;        
        rederr(lc).erel = 0.5 * (H.^2 - 1) * ediff;        
        rederr(lc).abs  = 0.5 * Es * rederr(lc).rel;
        rederr(lc).eabs = 0.5 * Es * rederr(lc).erel;
        if (inflc(lc) == 2)
           w    = (2 * pi * f).^4;
           rederr(lc).abs  = rederr(lc).abs .* w;
           rederr(lc).eabs = rederr(lc).eabs .* w;
        elseif (inflc(lc) == 4)
           w    = (2 * pi * f).^2;
           rederr(lc).abs  = rederr(lc).abs .* w;
           rederr(lc).eabs = rederr(lc).eabs .* w;
        end
        fprintf(fid, "  Upper bound on relative strain energy error");
        fprintf(fid, " (fmax = %7.2f Hz)\n", fmax);
        fprintf(fid, "     with static correction:    %10.4e\n",
                     rederr(lc).erel(end));
        fprintf(fid, "     without static correction: %10.4e\n\n",
                     rederr(lc).rel(end));
     end

  end

  mfs_paths("remove");

# Report an error if the error estimation is not valid

  if (err)
     fprintf(fid, "  Error estimation not possible: ");
     fprintf(fid, "fmax (%.2f Hz) > fp (%.2f Hz)\n", fmax, fp); 
     if (nargout > 0)
        text   = "mfs_reductionerror: error estimation not possible";
        mfs_errexit(text, fid);
     end
  end

# End

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

end
