function varargout = mfs_randresp(varargin)

# usage: varargout = mfs_randresp(varargin)
#
#  G      = randresp(Gxx, Gxy, HA, HB)
#  [G, f] = randresp(cmp, Gxx, Gxy, item, select, fpsd, ldcs, cflag)
#
# The function computes the power spectral density of the response
# from the power and the cross spectral densities of the excitation.
#
# In the first form, the transfer functions are supplied to the
# function. The frequencies of the spectral densities and the transfer
# functions must match.
#
# In the second form, the transfer functions are obtained from the
# results stored in the component. 
#
# Input  cmp        Structure with component data
#        Gxx(n, :)  Matrix with power spectral densities
#                   (columns correspond to frequencies,
#                    rows correspond to load cases)
#        Gxy(:, :)  Matrix with cross spectral densities
#                   (columns correspond to frequencies,
#                    for rows see below)
#        HA(n, :)   Matrix with transfer functions of response A
#                   (columns correspond to frequencies,
#                    rows correspond to load cases)
#        HB(n, :)   Matrix with transfer functions of response B
#                   (optional)
#        item       Name of response item
#        select     Item selector (see User Manual)
#        fpsd(:)    List of frequencies corresponding to the spectral
#                   densities (optional)
#        ldcs(n)    List of load cases corresponding to the spectral
#                   densities (optional)
#        cflag      Flag controlling computation of cross spectra
#                    0  do not compute cross spectra (default)
#                    1  compute cross spectra
#
# Output G(:, :)    Power spectral density of the response
#                   or cross spectral density G_AB
#        f          Frequencies (optional)
#
# Matrix Gxy contains the elements above the diagonal of the spectral
# matrix. The first n-1 rows contain the entries of the first row,
# the next n-2 rows those of the second row and so on.
#
# If HB is not defined, HB = HA is assumed and the power spectral
# density G_AA is computed. Otherwise, the cross spectral density
# G_AB is computed.
#
# If fpsd is not defined, it is assumed that the frequencies of the
# spectral densities match those of the transfer functions.
#
# If ldcs is not defined, it is assumed that the n columns of Gxx
# correspond to the first n load cases.
#
# If the selected response item is a nodal point result, then G
# is an array containing the power spectral densities. Columns
# correspond to frequencies and rows to response items.
#
# If the seleced response item is an element result, then G is a
# cell array of structures. The fields contain the power spectral
# densities of the corresponding element results.
#
# Cross spectral densities, if requested, are computed between element
# results at the same result point.
#
# ---------------------------------------------------------------------

# Copyright (c) 2024 by Johannes Wandinger

  t0 = clock();

# Branch according to format

  if (! isstruct(varargin{1}))

     if (nargin < 3 || nargin > 4 || nargout != 1)
        text1 = "mfs_randresp: correct call of function reads\n";
        text2 = "  G = mfs_randresp(Gxx, Gxy, H)\n";
        error([text1, text2]);
     endif

     Gxx = varargin{1}; Gxy = varargin{2}; HA = varargin{3};

     [nofldc, nfreq] = size(Gxx);

     if (nargin == 3)
        HB = [];
     else
        HB = varargin{4};
        if (rows(HB) != nofldc)
           text1 = sprintf("Number of load cases of %s does not match ",
                           inputname(1));
           text2 = sprintf("number of load cases of %s\n", inputname(4));
           error(["mfs_randresp: ", text1, text2]);
        endif
        if (columns(HB) != nfreq)
           text1 = sprintf("Number of frequencies of %s does not match ",
                           inputname(1));
           text2 = sprintf("number of frequencies of %s\n", inputname(4));
           error(["mfs_randresp: ", text1, text2]);
        endif
     endif

     if (rows(HA) != nofldc)
        text1 = sprintf("Number of load cases of %s does not match ",
                        inputname(1));
        text2 = sprintf("number of load cases of %s\n", inputname(3));
        error(["mfs_randresp: ", text1, text2]);
     endif

     if (columns(Gxy) != nfreq)
        text1 = sprintf("Number of frequencies of %s does not match ",
                        inputname(1));
        text2 = sprintf("number of frequencies of %s\n", inputname(2));
        error(["mfs_randresp: ", text1, text2]);
     endif

     if (columns(HA) != nfreq)
        text1 = sprintf("Number of frequencies of %s does not match ",
                        inputname(1));
        text2 = sprintf("number of frequencies of %s\n", inputname(3));
        error(["mfs_randresp: ", text1, text2]);
     endif

     if (rows(Gxy) != 0.5 * (nofldc - 1) * nofldc)
        text1 = sprintf("Number of columns of %s incompatible with ",
                        inputname(2));
        text2 = "number of load cases\n";
        error(["mfs_randresp: ", text1, text2]);
     endif

     time_stamp = ignore_function_time_stamp("all");
     mfs_paths("add", "mfs_randresp.m", "util");
     varargout{1} = mfs_randrespx(Gxx, Gxy, HA, HB);
     mfs_paths("remove");
     ignore_function_time_stamp(time_stamp);

  else

     if (nargout > 2 || nargin < 5 || nargout > 8)
        text1 = "mfs_randresp: correct call of function reads\n";
        text2 = "  [G, <f>] = mfs_randresp(cmp, Gxx, Gxy, item, sel, ";
        text3 = "<fpsd>, <ldcs>, <cflag>\n";
        error([text1, text2, text3]);
     endif

     cmp = varargin{1}; Gxx = varargin{2}; Gxy = varargin{3};

     [nofldc, nfreq] = size(Gxx);

     if (columns(Gxy) != nfreq)
        text1 = sprintf("Number of frequencies of %s does not match ",
                        inputname(2));
        text2 = sprintf("number of frequencies of %s\n", inputname(3));
        error(["mfs_randresp: ", text1, text2]);
     endif

     if (rows(Gxy) != 0.5 * (nofldc - 1) * nofldc)
        text1 = sprintf("Number of columns of %s incompatible with ",
                        inputname(2));
        text2 = "number of load cases\n";
        error(["mfs_randresp: ", text1, text2]);
     endif

     if (! isfield(cmp, "type"))
        error("mfs_randresp: \"%s\" is not a component\n",
              inputname(1));
     endif

     switch cmp.type

     case "solid"

        time_stamp = ignore_function_time_stamp("all");
        mfs_paths("add", "mfs_randresp.m", "solid");
        [G, f, rc] = mfs_randrespx(cmp, Gxx, Gxy, varargin(4 : end));
        mfs_paths("remove");
        ignore_function_time_stamp(time_stamp);

        if (rc)
           error("mfs_randresp: could not process al requests\n");
        else
           varargout{1} = G;
           if (nargout == 2)
              varargout{2} = f;
           endif
        endif

     otherwise

        error("mfs_randresp: model type \"%s\" not supported\n",
              cmp.type);

     endswitch

  endif

  elapsed_time = etime(clock(), t0);
  printf("%10.4f seconds needed to compute random response\n",
         elapsed_time);
  fflush(stdout);

endfunction
