function [resp, rc] = mfs_mtrans(cmp, dt, t, load, meth, ic, nwmp)

# usage  [resp, rc] = mfs_mtrans(cmp, dt, t, load, meth, ic, nwmp)
#
# Input  cmp       Structure with component data
#        dt        Time step
#        t(:)      Array with time
#        load(:)   Structure array defining loads
#        meth      Method: 2 = modal reduction
#                          3 = enhanced modal reduction
#        ic        Array with initial conditions
#        nwmp(2)   Array with Newmark parameters alpha and beta
#
# Output resp      Structure with transient response
#        rc        Return code: 0 means no errors
#
# This function performs a modal transient response analysis.
#
# ------------------------------------------------------------------------

# Copyright (c) 2021 by Johannes Wandinger

# Check arguments

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

# Initialize

  resp  = struct();
  rc    = 0;
  ntime = length(t);

# Check availability of data needed

  if (! isfield(cmp, "modes"))
     printf("*E* mfs_transresp: Normal modes not found\n");
     rc = 1;
  end
  if (! isempty(load) && ! isfield(cmp, "load"))
     printf("*E* mfs_transresp: Loads not found\n");
     rc = 1;
  end
  if (! isfield(cmp, "damping"))
     printf("*E* mfs_transresp: No damping defined\n");
     rc = 1;
  end

  if (rc) return; end

# Get load patterns and load functions

  if (! isempty(load))
     [pf, pu, ixf, ixu, rc] = mfs_trload(cmp, load);
     if (rc) return; end
     nf = length(ixf);
     nu = length(ixu);
     np = nf + nu;
     if (! isfield(load, "params"));
        load(1).params = [];
     end
  else
     nf = 0; nu = 0; np = 0;
     meth = 2;
  end

# Get information on degrees of freedom

  ndofg = cmp.dofs.ndofg;
  ndofl = cmp.dofs.ndofl;
  ndofp = cmp.dofs.ndofp;
  ndofd = cmp.dofs.ndofd;
  ndofr = cmp.stiff.ndofr;

  dofl = cmp.dofs.dofl;
  if (ndofp)
     dofp = cmp.dofs.dofp;
  end
  if (ndofd)
     dofd = cmp.dofs.dofd;
  end
  if (ndofr)
     dofe = cmp.stiff.dofe;
  end

# Get modal data

  w  = cmp.modes.omega;
  X  = cmp.modes.disp;
  nx = length(w);

# Partition stiffness and mass matrices

  if (nu)
     flags = [1, 1];
  else
     flags = [1, 0];
  end
  [Kll, Klp] = mfs_matpart(cmp.stiff.K, cmp.dofs, flags);
  Mll        = mfs_matpart(cmp.mass.M, cmp.dofs, [1, 0]);
  if (ndofr)
     Kee = Kll(dofe, dofe);
     Xr  = X(dofl, 1 : ndofr);
     XTM = Xr' * Mll;
  end 

# Get base motion

  if (nu)

     rhs = -Klp * pu;

     if (ndofr)
        rhs -= XTM' * (Xr' * rhs);
        ubl = zeros(ndofl, np, "double");
        ubl(dofe, :) = Kee \ rhs(dofe, :);
        ubl -= Xr * (XTM * ubl);
     else
        ubl = Kll \ rhs;
     end

     ub = zeros(ndofg, nu, "double");
     ub(dofl, :) = ubl; ub(dofp, :) = pu;

     if (ndofd)
        ub(dofd, :) = cmp.dofs.C(dofd, :) * ub;
     end

  end

# Get static mode shapes

  if (meth == 3)

     Vl = zeros(ndofl, np);
     Xs = zeros(ndofg, np);

     if (nf)
        Vl(:, ixf) = mfs_matpartr(pf, cmp.dofs);
     end

     if (nu)
        Mub = cmp.mass.M * ub;
        Vl(:, ixu) = -mfs_matpartr(Mub, cmp.dofs);
     end

     if (ndofr)
        Vl -= XTM' * (Xr' * Vl);
        Vl(dofe, :) = Kee \ Vl(dofe, :);
     else
        Vl = Kll \ Vl;
     end

     Vl -= X(dofl, :) * (X(dofl, :)' * (Mll * Vl));
     k = Vl' * Kll * Vl; m = Vl' * Mll * Vl;
     [qs, ws] = eig(k, m);
     Xs(dofl, :) = Vl * qs;
     if (ndofd)
        Xs(dofd, :) = cmp.dofs.C(dofd, :) * Xs;
     end

     w  = [w; sqrt(diag(ws))];
     X  = [X, Xs];
     nx = nx + np;

  end

# Get inertia loads due to base motion

  if (nu)
     Mgb = cmp.mass.M * ub;
     Mqb = X' * Mgb;
  end

# Get damping

  damping = cmp.damping;
  switch damping.type
  case "ratios"
     nd = length(damping.data);
     if (nd >= nx)
        D = damping.data(1 : nx)';
     else
        Dl = damping.data(nd);
        D  = [damping.data'; Dl(ones(nx - nd, 1))];
     end
     if (nu)
        Dqb = 2 * diag(2 .* D) * X' * cmp.mass.M * ub;
     end
  case "Rayleigh"
     if (ndofr)
        elmod     = ndofr + 1 : modes.nofmod;
        D         = zeros(modes.nofmod, 1);
        D(elmod)  = 0.5 * damping.data(1) * w(elmod);
        D(elmod) += 0.5 * damping.data(2) ./ w(elmod);
     else
        D = 0.5 * (damping.data(1) * w + damping.data(2) ./ w);
     end
     if (nu)
        Kgb  = cmp.stiff.K * ub;
        Dqb  = damping.data(1) * X' * Kgb;
        Dqb += damping.data(2) * Mqb;
     end
  end

# Get modal stiffness and damping

  dqq = 2 * w .* D;
  kqq = w .* w;

# Initialize excitation

  lq = zeros(nx, ntime, "double");

# Get contribution of loads to excitation

  if (nf)
     fn = zeros(nf, ntime, "double");
     for n = 1 : nf
         ix = ixf(n);
         func = str2func(load(ix).func);
         y    = func(t, load(ix).params);
         fn(n, :) = y(1, :);
     end
     lx = X' * pf;
     lq = lx * fn;
  end

# Get contribution of prescribed motion to excitation

  if (nu)
     fn   = zeros(nu, ntime, "double");
     fnd  = zeros(nu, ntime, "double");
     fndd = zeros(nu, ntime, "double");
     for n = 1 : nu
         ix = ixu(n);
         func = str2func(load(ix).func);
         y    = func(t, load(ix).params);
         fn(n, :)   = y(1, :);
         fnd(n, :)  = y(2, :);
         fndd(n, :) = y(3, :);
     end
     lq -= Mqb * fndd + Dqb * fnd;
  end

# Initialize response

  q   = zeros(nx, ntime);
  qd  = zeros(nx, ntime);
  qdd = zeros(nx, ntime);

# Get initial conditions

  qdd(:, 1) = lq(:, 1);

  if (! isempty(ic))
     if (columns(ic) == 1)
        if (nu)
           ic -= ub * fn(:, 1);
        end
        q(:, 1)    = X' * (cmp.mass.M * ic);   
        qdd(:, 1) -= kqq .* q(:, 1);
     else
        if (nu)
           ic(:, 1) -= ub * fn(:, 1);
           ic(:, 2) -= ub * fnd(:, 1);
        end
        qi = X' * (cmp.mass.M * ic);
        q(:, 1)    = qi(:, 1);
        qd(:, 1)   = qi(:, 2);
        qdd(:, 1) -= kqq .* q(:, 1) + dqq .* qd(:, 1);
     end
  end

# Compute response

  alpha = nwmp(1); beta = nwmp(2);

  a2 = 1 / (alpha * dt); a0 = a2 / dt; a1 = beta * a2; 
  a3 = 1 / (2 * alpha) - 1;
  a4 = beta / alpha - 1; a5 = 0.5 * dt * (a4 - 1);
  a6 = dt * (1 - beta); a7 = beta * dt;

  S = 1 ./ (kqq + a1 * dqq + a0);

  for n = 2 : ntime
      m = n - 1;
      q(:, n) = lq(:, n) + a0 * q(:, m) + a2 * qd(:, m) + a3 * qdd(:, m);
      q(:, n) += dqq .* (a1 * q(:, m) + a4 * qd(:, m) + a5 * qdd(:, m));
      q(:, n) = S .* q(:, n);
      qdd(:, n) = a0 * (q(:, n) - q(:, m)) - a2 * qd(:, m) - a3 * qdd(:, m);
      qd(:, n) = qd(:, m) + a6 * qdd(:, m) + a7 * qdd(:, n);
  end

# Build transient response structure

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

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

  if (meth == 3)
     resp.Xs = Xs;
  end
               
end
