function [disp, w] = mfs_getdispf(resp, modes, ixdof, ndofr, absdisp)

# usage: [disp, w] = mfs_getdispf(resp, modes, ixdof, ndofr, absdisp)
#
# Input  resp         Structure with frequency response results
#        modes(:, :)  Matrix with normal modes
#        ixdof(:)     Array with degree of freedom indices
#        ndofr        Number of rigid body modes
#        absdisp      = 1: compute absolute displacements
#                     = 0: compute relative displacements
#
# Output disp(:, :)   Array with displacements
#                     (columns correspond to frequencies)
#        w(:)         Array with circular frequencies (optional)
#
# The function returns the frequency-dependent displacements of selected
# degrees of freedom.
#
# -------------------------------------------------------------------------

# Check arguments

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

# Get circular frequency

  w = 2 * pi * resp.freq;

# Get displacements depending on method and load type

  switch resp.method

  case 1  % Direct Frequency Response

     disp = resp.U(ixdof, :);

  case 4  % Force Summation Method

     disp = resp.Ue(ixdof, :);
     if (absdisp && ndofr)
        disp += modes(ixdof, 1 : ndofr) * resp.Q(1 : ndofr, :); 
     end

  otherwise % Modal Frequency response

     disp = modes(ixdof, :) * resp.Q;

     switch resp.type

     case 1 % Force or moment

        if (resp.method == 3)
           disp += resp.mct(ixdof);
        end

     case 2 % Prescribed displacement

        disp += resp.base(ixdof);
        if (resp.method == 3)
           disp -= resp.mct(ixdof) .* (w.^2);
        end

     case 3  % Prescribed velocity

        disp -= i * resp.base(ixdof) ./ w;
        if (resp.method == 3)
           disp -= resp.mct(ixdof) .* w;
        end

     case 4 % Prescribed acceleration

        disp -= resp.base(ixdof) ./ (w.^2);
        if (resp.method == 3)
           disp += resp.mct(ixdof);
        end

     end % Load type

  end % Method

end
