function [D, Dqb] = mfs_mdamp(damping, stiff, mass, modes, base, dofl)

# usage  [D, Dqb] = mfs_mdamp(damping, stiff, mass, modes, enfmot)
#
# Input  damping    Structure with damping data
#        stiff      Structure with stiffness matrix
#        mass       Structure with mass matrix
#        modes      Structure with normal modes 
#        base(:, :) Matrix with base motion due to enforced motion
#        dofl(:)    Indices of local degrees of freedom 
#
# Output D(:)       Modal damping ratios
#        Dqb(:, :)  Modal damping loads from base motion
#
# The function computes the damping matrices needed in a modal response
# analysis.
#
# ------------------------------------------------------------------------

  enfmot = ! isempty(base);
  if (! enfmot)
     Dqb = [];
  end

 switch damping.type
  case "ratios"
     nx = modes.nofmod;
     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 (enfmot)
        w   = modes.omega;
        X   = modes.disp;
        Dqb = 2 * diag(w .* D) * X(dofl, :)' * Mll * base(dofl, :);
     end
  case "Rayleigh"
     ndofr = stiff.ndofr;
     w     = modes.omega;
     X     = modes.disp;
     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 (enfmot)
        Kqg = X' * cmp.stiff.K;
        Mqg = X' * cmp.mass.M;
        Dqb = damping.data(1) * Kqg * base;
        Dqb = Dqb + damping.data(2) * Mqg * base;
     end
  end

end
