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

# usage  [resp, rc] = mfs_dtrans(cmp, dt, t, load, ic, nwmp)
#
# Input  cmp       Structure with component data
#        dt        Time step
#        t(:)      Array with time
#        load(:)   Structure array defining loads
#        ic        Array with initial conditions
#        nwmp(3)   Array with Newmark parameters alpha and beta and
#                  factor gamma for time step to compute initial
#                  accelerations
#
# Output resp      Structure with transient response
#        rc        Return code: 0 means no errors
#
# This function performs a direct transient response analysis using the
# Newmark method.
#
# ------------------------------------------------------------------------

# Copyright (c) 2021 by Johannes Wandinger

# Check arguments

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

# Initialize

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

# Check availability of data needed

  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;
  elseif (! strcmp(cmp.damping.type, "Rayleigh"))
     printf("*E* mfs_transresp: Damping type %s ignored in direct response\n",
            cmp.damping.type);
     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;
  end

# Get information on degrees of freedom

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

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

# 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, Mlp] = mfs_matpart(cmp.mass.M, cmp.dofs, flags);

# Get damping

  Dll = Kll * cmp.damping.data(1) + Mll * cmp.damping.data(2);
  if (nu)
     Dlp = Klp * cmp.damping.data(1) + Mlp * cmp.damping.data(2);
  end

# Initialize excitation

  ll = zeros(ndofl, 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
     ll = mfs_matpartr(pf, cmp.dofs) * fn;
  end

# Get contribution of prescribed motion to excitation

  if (nu)
     gn   = zeros(nu, ntime, "double");
     gnd  = zeros(nu, ntime, "double");
     gndd = zeros(nu, ntime, "double");
     for n = 1 : nu
         ix = ixu(n);
         func = str2func(load(ix).func);
         y    = func(t, load(ix).params);
         gn(n, :)   = y(1, :);
         gnd(n, :)  = y(2, :);
         gndd(n, :) = y(3, :);
     end
     up = pu * gn; vp = pu * gnd; ap = pu * gndd;
     ll -= Klp * up + Dlp * vp + Mlp * ap;
  end

# Initialize response

  ul = zeros(ndofl, ntime, "double");
  vl = zeros(ndofl, ntime, "double");
  al = zeros(ndofl, ntime, "double");

# Get initial displacements and velocities

  if (! isempty(ic))
     if (columns(ic) == 1)
        ul(:, 1) = ic(dofl, 1);
     else
        ul(:, 1) = ic(dofl, 1);
        vl(:, 1) = ic(dofl, 2);
     end
  end

# Get initial accelerations

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

  dtm = gamma * dt;
  a1 = 2 / dtm; a0 = a1 / dtm;

  Kdll = Kll + a0 * Mll + a1 * Dll;
  ul2  = ll(:, 1) + gamma * (ll(:, 2) - ll(:, 1));
  ul2 += Mll * (a0 * ul(:, 1) + a1 * vl(:, 1));
  ul2 += Dll * (a1 * ul(:, 1) + vl(:, 1));
  ul2  = Kdll \ ul2;
  al2  = a0 * (ul2 - ul(:, 1)) - a1 * vl(:, 1);
  vl2  = vl(:, 1) + al2 * dtm;

  b1 = 1 / ((1 - beta) * dtm); b2 = beta / (1 - beta);
  al(:, 1) = b1 * (vl2 - vl(:, 1)) - b2 * al2;

# Compute response

  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;

  Kdll = Kll + a0 * Mll + a1 * Dll;
  [L, U, p, q] = lu(Kdll, "vector");

  for n = 2 : ntime
      m = n - 1;
      ul(:, n) = ll(:, n);
      ul(:, n) += Mll * (a0 * ul(:, m) + a2 * vl(:, m) + a3 * al(:, m));
      ul(:, n) += Dll * (a1 * ul(:, m) + a4 * vl(:, m) + a5 * al(:, m));
      ulx = L \ ul(p, n);
      ulx = U \ ulx;
      ul(q, n) = ulx;
      al(:, n) = a0 * (ul(:, n) - ul(:, m)) - a2 * vl(:, m) - a3 * al(:, m);
      vl(:, n) = vl(:, m) + a6 * al(:, m) + a7 * al(:, n);
  end

  clear ll;

  u = zeros(ndofg, ntime);
  u(dofl, :) = ul;
  if (nu)
     u(dofp, :) = up;
  end
  if (ndofd)
     u(dofd, :) = cmp.dofs.C(dofd, :) * u;
  end
  clear ul; clear up;

  v = zeros(ndofg, ntime);
  v(dofl, :) = vl;
  if (nu)
     v(dofp, :) = vp;
  end
  if (ndofd)
     v(dofd, :) = cmp.dofs.C(dofd, :) * v;
  end
  clear vl; clear vp;

  a = zeros(ndofg, ntime);
  a(dofl, :) = al;
  if (nu)
     a(dofp, :) = ap;
  end
  if (ndofd)
     a(dofd, :) = cmp.dofs.C(dofd, :) * a;
  end
  clear al; clear ap;

# Build transient response structure

  resp = struct("method", 1, "ntime", ntime, "tsteps", t,
                "nback", ntime, "tback", t, "u", u, "v", v, "a", a);

end
