function [G, f, rc] = mfs_randrespx(cmp, Gxx, Gxy, args)

# usage: [G, f, rc] = mfs_randrespx(cmp, Gxx, Gxy, args)
#
# Input  cmp        Structure with component data
#        Gxx(n, :)  Matrix with power spectral densities
#        Gxy(:, :)  Matrix with cross spectral densities
#        args{1}    Response item
#        args{2}    Nodal point or element selector
#        args{3}    List of spectral density frequencies (optional)
#        args{4}    List of load cases (optional)
#        args{5}    Flag controlling computation of cross spectra
#                   (optional)
#
# Output G(:, :)    Response spectra
#        f(:)       List of frequencies of response spectra
#        rc         Return code: 0 = no errors
#                                1 = errors
#
# For nodal point results, G is a matrix with rows corresponding to
# the selected responses and columns corresponding to frequencies.
#
# For element results, G is a cell array of structures, with each cell
# containing the results of one element. The structures contain the
# same fields as those of a frequency response analysis. The fields
# contain the power spectral densities. Additional fields contain the
# cross spectral densities between the different element results. The
# names of these fields are of the form 'field1'_'field2'.
# 
# --------------------------------------------------------------------

# Initialize

  G = []; f = []; rc = 1;

# Check which arguments are defined

  nargs = length(args);
  if (nargs < 2 || nargs > 5)
     error("mfs_randresp: wrong number of input arguments (%d)\n",
           nargs);
  endif

  item = args{1}; sel = args{2};

  is_fpsd = 0; is_ldcs = 0; cflag = 0;

  if (nargs > 2)
     is_fpsd = ! isempty(args{3});
     if (is_fpsd)
        fpsd = args{3};
     endif
  endif

  if (nargs > 3)
     is_ldcs = ! isempty(args{4});
     if (is_ldcs)
        ldcs = args{4};
     endif
  endif

  if (nargs == 5)
     cflag = ! isempty(args{5});
  endif

# Check load cases

  [nofldc, nfpsd] = size(Gxx);

  if (is_ldcs)
     if (length(ldcs) != nofldc)
        printf("*E* mfs_randresp: incompatible number of load cases\n");
        printf("    Gxx has %d columns, size of ldcs is %d\n",
               nofldc, length(ldcs));
        return;
     endif
  else
     ldcs = 1 : nofldc;
  endif

# Check frequencies of transfer functions

  for l = 1 : nofldc
      [out, rc] = mfs_getrespx(cmp, "freqresp", "freq", {ldcs(l)});
      if (rc) return; endif
      f     = out{1};
      nf(l) = length(f);
      ff{l} = f;
  endfor

  nfreq = nf(1); f = ff{1};

  rc = 1;

  for l = 2 : nofldc
      if (nf(l) != nfreq)
         printf("*E* mfs_randresp: frequencies of load cases %d\n");
         printf("    and %d do not match\n", ldcs(1), ldcs(l));
         return;
      else
         ix = find(ff{l} - f);
         if (ix)
            printf("*E* mfs_randresp: frequencies of load cases %d\n");
            printf("    and %d do not match\n", ldcs(1), ldcs(l));
            return;
         endif
      endif
  endfor

# Interpolate spectra to frequencies of transfer functions

  if (is_fpsd)
     for l = 1 : rows(Gxx)
        Gxxf(l, :) = interp1(fpsd, Gxx(l, :), f);
     endfor
     for l = 1 : rows(Gxy)
        Gxyf(l, :) = interp1(fpsd, Gxy(l, :), f);
     endfor
     Gxx = Gxxf;
     Gxy = Gxyf;
  elseif (nfpsd != nfreq)
     printf("*E* mfs_randresp: frequencies of spectra not compatible\n")
     printf("    with frequencies of transfer functions\n");
     return;
  endif

# Branch according to response type

  switch item

  case {"disp", "velo", "acce"}

     % Nodal point results

     for l = 1 : nofldc
         [out, rc] = mfs_getrespx(cmp, "freqresp", item,
                                  {sel, ldcs(l)});
         if (rc) return;  endif
         H(:, l, :) = out{1};
     endfor

     nresp = size(H, 1);
     G     = zeros(nresp, nfreq);

     for k = 1 : nresp
         l = 1;
         h = squeeze(H(k, :, :));
         for m = 1 : nofldc
             G(k, :) += Gxx(m, :) .* (h(m, :) .* conj(h(m, :)));
             for n = m + 1 : nofldc
                 G(k, :) += 2 * real(h(m, :) .* Gxy(l++, :) .* conj(h(n, :)));
             endfor
         endfor
     endfor

  case {"stress", "resultant"}

     % Element results

     for l = 1 : nofldc
         [out, rc] = mfs_getrespx(cmp, "freqresp", item,
                                  {sel, ldcs(l)});
         if (rc) return;  endif
         H{l} = out{1};
     endfor

     nresp = length(out{1});
     G     = cell(nresp, 1);

     for k = 1 : nresp
         for l = 1 : nofldc
             HS(l) = H{l}{k};
         endfor
         G{k} = eltspectra(Gxx, Gxy, HS, nofldc, cflag);
     endfor

  otherwise

     prinf("*E* mfs_randresp: item \"%s\" not supported\n", item);

  endswitch

endfunction

function G = eltspectra(Gxx, Gxy, HS, nofldc, cflag)

# usage: G = eltspectra(Gxx, Gxy, HS, nofldc, cflag)
#
# Input  Gxx(nofldc, :)  Matrix with power spectral densities
#        Gxy(:, :)       Matrix with cross spectral densities
#        HS(nofldc)      Structure array with transfer functions
#        nofldc          Number of load cases
#        cflag           0 do not compute cross spectral densities
#                        1 compute cross spectral densities
# Output G               Structure with power and cross spectral
#                        densities
#
# The function computes the power and cross spectral densities of the
# results of one element. Results are store in structure G which has
# the same fields a for frequency response results.
#
# Field names of cross spectral densities are a combination of the
# field names of the two results, with "_" between them.
#
# --------------------------------------------------------------------

# Initialize output structure

  G     = HS(1);
  names = fieldnames(G);
  if (G.rtype < 3)
     f1 = 6;     % Skip header
  else
     f1 = 7;     % Skip header and TE
  endif
  rnames = names(f1 : end);

# Compute power spectral densities

  nofres = length(rnames);
  nofpnt = G.nofpnt;
  nfreq  = G.nofcol;

  for k = 1 : nofres
      for l = 1 : nofldc
          H(:, l, :) = getfield(HS(l), rnames{k});
      endfor
      Gyy = zeros(nofpnt, nfreq);
      for p = 1 : nofpnt
          l = 1;
          h = squeeze(H(p, :, :));
          for m = 1 : nofldc
              Gyy(p, :) += Gxx(m, :) .* (h(m, :) .* conj(h(m, :)));
              for n = m + 1 : nofldc
                  Gyy(p, :) += 2 * real(h(m, :) .* Gxy(l++, :) .* conj(h(n, :)));
              endfor
          endfor
          G.(rnames{k}) = Gyy;
      endfor
  endfor

  if (! cflag) return; endif

# Compute cross spectral densities

  for kA = 1 : nofres

      for l = 1 : nofldc
          HA(:, l, :) = getfield(HS(l), rnames{kA});
      endfor

      for kB = kA + 1 : nofres

         for l = 1 : nofldc
             HB(:, l, :) = getfield(HS(l), rnames{kB});
         endfor

         rname = [rnames{kA}, "_", rnames{kB}];

         GAB = zeros(nofpnt, nfreq);

         for p = 1 : nofpnt
             l = 1;
             hA = squeeze(HA(p, :, :));
             hB = squeeze(HB(p, :, :));
             for m = 1 : nofldc
                 GAB(p, :) += Gxx(m, :) .* (hA(m, :) .* conj(hB(m, :)));
                 for n = m + 1 : nofldc
                     GAB(p, :) += hA(m, :) .* Gxy(l, :) .* conj(hB(n, :));
                     GAB(p, :) += hA(n, :) .* conj(Gxy(l++, :)) .* conj(hB(m, :));
                 endfor
             endfor
         endfor

         G.(rname) = GAB;

      endfor

  endfor

endfunction
