function [cmp, rc] = mfs_freq2timex(cmp, spectra, rcase, dt, bc)

# usage: [cmp, rc] = mfs_freq2timex(cmp, spectra, rcase, dt, bc)
#
# Input  cmp         Structure with component
#        spectra     Structure with spectra
#        rcase       Result case identifier
#        dt          Time step (optional)
#        bc          Base correction (optional)
#                    0 = no correction (default)
#                    1 = subtract a constant
#                    2 = subtact a linear function
#
# Output cmp         Structure with component
#        rc          Return code: 0 = no errors
#                                 1 = errors
#
# Fields of structure spectra
#
#   df          Frequency step
#   lc(:)       Identifiers of load cases with transfer functions
#   spec(:, :)  Spectra assigned to load cases: rows correspond to
#               load cases, columns to frequencies
#
# The function transforms results from a frequency response analysis
# to the time domain. The results are stored as results of a transient
# response analysis.
#
# ---------------------------------------------------------------------

# Copyright (c) 2022 by Johannes Wandinger

  required_fields = {"df", "lc", "spec"};

  if (nargin != 5 || nargout != 2)
     print_usage();
  end

# Initialize

  rc    = 0;
  ndofg = cmp.dofs.ndofg;

# Check spectra

  defined_fields = fieldnames(spectra);

  reqdef = ismember(required_fields, defined_fields);
  defleg = ismember(defined_fields, required_fields);

  undef = find(! reqdef); lnundef = length(undef);
  illeg = find(! defleg); lnilleg = length(illeg);

  if (lnundef)
     rc = 1;
     for n = 1 : lnundef
         printf("*E* mfs_freq2time: spectra: required field %s undefined\n",
                required_fields{undef(n)});
     end
  end

  if (lnilleg)
     rc = 1;
     for n = 1 : lnilleg
         printf("*E* mfs_freq2time: spectra: illegal field %s found\n",
                defined_fields{illeg(n)});
     end
  end

  if (rc) return; end 

  ldcs  = spectra.lc;
  noflc = length(ldcs);
  nspec = rows(spectra.spec);

  if (noflc != nspec)
     printf("*E* mfs_freq2time: number of load cases (%2.0d) ", noflc);
     printf("does not match number of spectra (%2.0d)\n",
            nspec); 
     rc = 1;
     return
  end
  
# Check existence of load cases referenced, collect information on method,
# maximum frequency and base motion

  meths = zeros(1, nspec);
  fm    = zeros(1, nspec);
  nbase = 0;   % Number of base motions
  nsms  = 0;   % Number of static modes shapes

  for n = 1 : nspec
      lc = ldcs(n);
      if (lc < 1 || lc > cmp.freqresp.nfa ||
          isempty(cmp.freqresp.resp{lc}))
         printf("*E* mfs_freq2time: load case %d referenced by", lc);
         printf(" spectrum %d does not exist\n", n);
         rc = 1;
      else
         meths(n) = cmp.freqresp.resp{lc}.method;
         fm(n)    = cmp.freqresp.resp{lc}.freq(end);
         nbase += isfield(cmp.freqresp.resp{lc}, "base");
         nsms  += meths(n) == 3;
      end
  end

  if (rc) return; end

# Check compatibility of methods

  meth = meths(1);

  if (meth == 1 && sum(ismember(2 : 4, meths)))
     printf("*E* mfs_freq2time: direct frequency response cannot be");
     printf(" mixed with other methods\n");
     rc = 1;
  end

  if (meth == 4 && sum(ismember(1 : 3, meths)))
     printf("*E* mfs_freq2time: force summation method cannot be");
     printf(" mixed with other methods\n");
     rc = 1;
  end

  meth = max(meths);

  if (rc) return; end

# Get highest frequency up to which all transfer functions are defined

  fmax = min(fm);

# Get frequencies corresponding to spectra up to fmax

  f = 0 : spectra.df : fmax; nf = length(f);

# Extract spectrum values needed

  spec = spectra.spec(:, 1 : nf);

# Get time step and number of time steps

  T = 1 / spectra.df;

  if (isempty(dt))
     ntime = 2 * nf - 1;
     dt    = T / ntime;
  else
     ntime = round(T / dt);
  end
  tsteps = (0 : ntime - 1) * dt;

# Get w to compute derivatives

  w = 2 * pi * i * f;
  if (meth == 3)
     wq = w .* w;
  end

# Branch according to method and perform IFFTs

  switch meth

  case 1        % Direct frequency response analysis

     U = zeros(ndofg, nf, "double");

     for n = 1 : nspec
         lc = ldcs(n);
         U += mfs_addspec(f, spec(n, :), cmp.freqresp.resp{lc}.freq,
                          cmp.freqresp.resp{lc}.U);
     end

     U(:, 2 : end) *= 2;

     UD = w .* U; UDD = w .* UD;

     u = real(ifft(U, ntime, 2)) / dt;
     v = real(ifft(UD, ntime, 2)) / dt;
     a = real(ifft(UDD, ntime, 2)) / dt;

     if (bc)
        u -= u(:, 1);
        if (bc == 2)
           u -= v(:, 1) * tsteps;
           v -= v(:, 1);
        end
     end
         
     resp = struct("method", 1, "ntime", ntime, "tsteps", tsteps,
                   "nback", ntime, "tback", tsteps, 
                   "u", u, "v", v, "a", a);
 
  case {2, 3}   % Modal frequency response analysis

     nofmod = cmp.modes.nofmod;

     Q = zeros(nofmod, nf, "double");

     if (nbase)
        base = zeros(ndofg, nbase, "double");
        fn   = zeros(nbase, ntime, "double");
        fnd  = zeros(nbase, ntime, "double");
        fndd = zeros(nbase, ntime, "double");
     end
     if (nsms)
        Xs   = zeros(ndofg, nsms, "double");
        qs   = zeros(nsms, ntime, "double");
        qsd  = zeros(nsms, ntime, "double");
        qsdd = zeros(nsms, ntime, "double");
     end

     ixb = 0; ixsm = 0;

     for n = 1 : nspec

         lc = ldcs(n);
         Q += mfs_addspec(f, spec(n, :), cmp.freqresp.resp{lc}.freq,
                          cmp.freqresp.resp{lc}.Q);

         if (isfield(cmp.freqresp.resp{lc}, "base"))
            base(:, ++ixb) = cmp.freqresp.resp{lc}.base;
            type = cmp.freqresp.resp{lc}.type; 
            d = zeros(3, nf);
            for m = 1 : 3
                k = m - type + 1;
                if (k != 0)
                   d(m, 2 : nf) = w(2 : nf).^k;
                else
                   d(m, :) = ones(1, nf); 
                end
            end
            S = spec(n, :); S(2 : nf) *= 2;
            fn(ixb, :)   = real(ifft(S .* d(1, :), ntime, 2)) / dt;
            fnd(ixb, :)  = real(ifft(S .* d(2, :), ntime, 2)) / dt;
            fndd(ixb, :) = real(ifft(S .* d(3, :), ntime, 2)) / dt;
         end

         if (meths(n) == 3)
            Xs(:, ++ixsm) = cmp.freqresp.resp{lc}.mct;
            if (cmp.freqresp.resp{lc}.type == 1)
               S = spec(n, :); S(2 : nf) *= 2;
               qs(ixsm, :)   = real(ifft(S, ntime, 2)) / dt;
               qsd(ixsm, :)  = real(ifft(S .* w, ntime, 2)) / dt;
               qsdd(ixsm, :) = real(ifft(S .* wq, ntime, 2)) / dt;
            else
               qs(ixsm, :)   = fn(ixb, :);
               qsd(ixsm, :)  = fnd(ixb, :);
               qsdd(ixsm, :) = fndd(ixb, :);
            end
         end

     end

     Q(:, 2 : end) *= 2;
     QD = w .* Q; QDD = w .* QD;

     q   = real(ifft(Q, ntime, 2)) / dt;
     qd  = real(ifft(QD, ntime, 2)) / dt;
     qdd = real(ifft(QDD, ntime, 2)) / dt;

     if (nsms)
        q = [q ; qs]; qd = [qd ; qsd]; qdd = [qdd ; qsdd];
     end

     if (bc)
        rigmod = 1 : cmp.stiff.ndofr;
        q(rigmod, :) -= q(rigmod, 1);
        if (nbase)
           fn -= fn(:, 1);
        end
        if (bc == 2)
           q(rigmod, :)  -= qd(rigmod , 1) * tsteps;
           qd(rigmod, :) -= qd(rigmod, 1);
           if (nbase)
              fn  -= fnd(:, 1) * tsteps;
              fnd -= fnd(:, 1);
           end
        end
     end

     resp = struct("method", meth, "ntime", ntime, "tsteps", tsteps,
                   "nback", 0, "q", q, "qd", qd, "qdd", qdd);

     if (nbase)
        resp.base = base;
        resp.f = fn; resp.fd = fnd; resp.fdd = fndd;
     end

     if (nsms)
        resp.Xs = Xs;
     end

  case 4        % Force summation method

     nofmod = cmp.modes.nofmod;

     Q  = zeros(nofmod, nf, "double");
     Ue = zeros(ndofg, nf, "double");

     for n = 1 : nspec
         lc = ldcs(n);
         Q  += mfs_addspec(f, spec(n, :), cmp.freqresp.resp{lc}.freq,
                           cmp.freqresp.resp{lc}.Q);
         Ue += mfs_addspec(f, spec(n, :), cmp.freqresp.resp{lc}.freq,
                           cmp.freqresp.resp{lc}.Ue);
     end

     Q(:, 2 : end) *= 2; Ue(:, 2 : end) *= 2;

     QD  = w .* Q; QDD = w .* QD;

     q   = real(ifft(Q, ntime, 2)) / dt;
     qd  = real(ifft(QD, ntime, 2)) / dt;
     qdd = real(ifft(QDD, ntime, 2)) / dt;

     UeD = w .* Ue; UeDD = w .* UeD;

     ue = real(ifft(Ue, ntime, 2)) / dt;
     ve = real(ifft(UeD, ntime, 2)) / dt;
     ae = real(ifft(UeDD, ntime, 2)) / dt;

     if (bc)
        rigmod = 1 : cmp.stiff.ndofr;
        q(rigmod, :)  -= q(rigmod, 1);
        if (bc == 2)
           q(rigmod, :)  -= qd(rigmod, 1) * tsteps;
           qd(rigmod, :) -= qd(rigmod, 1);
        end
     end
         
     resp = struct("method", 4, "ntime", ntime, "tsteps", tsteps,
                   "nback", 0, "q", q, "qd", qd, "qdd", qdd,
                   "ue", ue, "ve", ve, "ae", ae);

  end

# Store results in component

  if(! isfield(cmp, "transresp"))
     cmp.transresp.nta = rcase;
  else
     cmp.transresp.nta = max(rcase, cmp.transresp.nta);
  end
  cmp.transresp.resp{rcase} = resp;

end

function Y = mfs_addspec(fspec, spec, freq, tf)

# This function interpolates the transfer function matrix to the
# frequencies of the spectrum and multiplies it with the spectrum.
# If the zero frequency is not included, a value of zero is
# added.
#
# Input  fspec(:) Frequencies correspondig to spectrum
#        spec(:)  Spectrum
#        freq(:)  Frequencies corresponding to transfer matrix
#        tf(:, :) Transfer matrix
#
# Output Y(:, :)  Values of transfer matrix multplied by
#                 corresponding values of spectrum
#
# ---------------------------------------------------------------------

# Add zero frequency if necessary and transpose transfer matrix
# (interp1 interpolates along columns)

  if (freq(1) > 0)
     freq           = [0, freq];
     nf             = length(freq);
     tft(2 : nf, :) = tf.';
  else
     tft = tf.';
  end

# Interpolate

  Y = interp1(freq, tft, fspec).' .* spec;

end
