# Example: Box Beam
#          Comparison of results from direct, modal and enhanced
#          modal frequency response analysis
#          Only loadcase 1 is considered.
#
# --------------------------------------------------------------
  %          green     red      blue   magenta
  colors = [0, 1, 0; 1, 0, 0; 0, 0, 1; 1, 0, 1];
  set(0, "defaultaxescolororder", colors);

  fid = fopen("compare.res", "wt");

# Output device depending on operating system

  if (isunix)
     DEVICE = "svg";
  else
     DEVICE = "png";
  endif
  FORMAT = ["-d", DEVICE];
  EXT    = [".", DEVICE];

# Results from normal modes analysis

  load modes.bin

# Elements for output

  elts = {"C", "D"};

# Degrees of freedom for displacements

  rid = {"A", 3; "B", 3};

# Excitation frequencies

  fmin  =  0;  % Zero frequency = static solution
  df    = 10;

  fe = fmin : df : fmax;  % fmax from modes.bin

# Error estimation

  rederr = mfs_reductionerror(fid, beam, fe, 1);

# Frequency response analyses

  beamd = mfs_freqresp(beam, fe, "loadcase", 1,
                       "method", "direct");
  beamm = mfs_freqresp(beam, fe, "nband", 0, "loadcase", 1,
                       "method", "modal");
  beame = mfs_freqresp(beam, fe, "nband", 0, "loadcase", 1);

# Results

  Ed = mfs_getresp(beamd, "freqresp", "strnegy", 1);
  Em = mfs_getresp(beamm, "freqresp", "strnegy", 1);
  Ee = mfs_getresp(beame, "freqresp", "strnegy", 1);

  Ud = mfs_getresp(beamd, "freqresp", "disp", rid, 1);
  Um = mfs_getresp(beamm, "freqresp", "disp", rid, 1);
  Ue = mfs_getresp(beame, "freqresp", "disp", rid, 1);

  sigd = mfs_getresp(beamd, "freqresp", "stress", elts);
  sigm = mfs_getresp(beamm, "freqresp", "stress", elts);
  sige = mfs_getresp(beame, "freqresp", "stress", elts);

  for n = 1 : length(elts)
      sigxd(n, :) = sigd{n}.sigxu;
      sigxm(n, :) = sigm{n}.sigxu;
      sigxe(n, :) = sige{n}.sigxu;
  endfor

# Static results

  Es = 2 * Ed(1);
  us = abs(Ud(:, 1)); sigxs = abs(sigxd(:, 1));

# Error in strain energy

  em = abs(Em(2 : end) - Ed(2 : end)) / Es;
  ee = abs(Ee(2 : end) - Ed(2 : end)) / Es;
  f  = fe(2 : end);

  figure(1, "position", [100, 400, 750, 500],
            "paperposition", [0, 0, 15, 10]);
  semilogy(f, rederr.rel(2 : end), f, em,
           f, rederr.erel(2 : end), f, ee);
  legend("Estimated modal", "Computed modal",
         "Estimated enhanced", "Computed enhanced",
         "location", "southeast");
  grid;
  xlabel('f [Hz]'); 
  print(["see", EXT], FORMAT, "-F:10");

# Displacement error

  eum = abs(Um(:, 2 : end) - Ud(:, 2 : end)) ./ us;
  eue = abs(Ue(:, 2 : end) - Ud(:, 2 : end)) ./ us;

  figure(2, "position", [300, 400, 750, 500],
            "paperposition", [0, 0, 15, 10]);
  semilogy(f, eum(1, :), f, eum(2, :),
           f, eue(1, :), f, eue(2, :));
  legend("A modal", "B modal",
         "A enhanced", "B enhanced",
         "location", "southeast");
  grid;
  xlabel('f [Hz]'); 
  ylabel('\Delta U_z / u_S');
  print(["eu", EXT], FORMAT, "-F:10");

# Stress error

  esm = abs(sigxm(:, 2 : end) - sigxd(:, 2 : end)) ./ sigxs;
  ese = abs(sigxe(:, 2 : end) - sigxd(:, 2 : end)) ./ sigxs;

  figure(3, "position", [500, 400, 750, 500],
            "paperposition", [0, 0, 15, 10]);
  semilogy(f, esm(1, :), f, esm(2, :),
           f, ese(1, :), f, ese(2, :));
  legend("C modal", "D modal",
         "C enhanced", "D enhanced",
         "location", "southeast");
  grid;
  xlabel('f [Hz]'); 
  ylabel('\Delta \sigma_x / \sigma_{Sx}');

  print(["esig", EXT], FORMAT, "-F:10");

  fclose(fid);
