# Example: Box Beam
#          Enhanced modal frequency response analysis
#
# -------------------------------------------------------------------

  colors = [0, 1, 0; 1, 0, 0; 0, 0, 1]; % green; red; blue
  set(0, "defaultaxescolororder", colors);

# 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

# Excitation frequencies

  fmin = 20;  % Minimum excitation frequency
  df   =  5;  % Frequency step
  nb   =  7;  % Additional frequencies in halfpower bandwidth

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

# Frequency response analysis

  [beam, f] = mfs_freqresp(beam, fe, "nband", nb, "loadcase", 1);
  beam      = mfs_freqresp(beam, f, "nband", 0, "loadcase", 2);

# Transfer functions for displacements

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

  Uz1 = mfs_getresp(beam, "freqresp", "disp", rid, 1);
  Uz2 = mfs_getresp(beam, "freqresp", "disp", rid, 2);

  idA = mfs_getset(beam, "nset", "A");
  idB = mfs_getset(beam, "nset", "B");
  lgtext = {sprintf("Point A (%d)", idA), ...
            sprintf("Point B (%d)", idB)};

  figure(1, "position", [100, 200, 1000, 750],
            "paperposition", [0, 0, 15, 15]);
  subplot(2, 1, 1)
     semilogy(f, abs(Uz1));
     title("Loadcase 1");
     legend(lgtext);
     grid;
     ylabel('|U_z| / F [mm/N]'); 
  subplot(2, 1, 2)
     semilogy(f, abs(Uz2));
     title("Loadcase 2");
     grid;
     xlabel('f [Hz]');
     ylabel('|U_z| / a_z [mm/g]'); 
  print(["U", EXT], FORMAT, "-F:10");

# Transfer functions for stresses

  sig1 = mfs_getresp(beam, "freqresp", "stress", {"C", "D"}, 1);
  sig2 = mfs_getresp(beam, "freqresp", "stress", {"C", "D"}, 2);

  sigx1 = [sig1{1}.sigxu; sig1{2}.sigxu];
  sigx2 = [sig2{1}.sigxu; sig2{2}.sigxu];

  lgtext = {sprintf("Element C (%d)", sig1{1}.id), ...
            sprintf("Element D (%d)", sig1{2}.id)}; 

  figure(2, "position", [300, 200, 1000, 750],
            "paperposition", [0, 0, 15, 15]);
  subplot(2, 1, 1)
     semilogy(f, abs(sigx1));
     title("Loadcase 1");
     legend(lgtext);
     grid;
     ylabel('|\sigma_x| / F [MPa/N]'); 
  subplot(2, 1, 2)
     semilogy(f, abs(sigx2));
     title("Loadcase 2");
     grid;
     xlabel('f [Hz]');
     ylabel('|\sigma_x| / a_z [MPa/g]'); 
  print(["sig", EXT], FORMAT, "-F:10");

# Deformation and stress at peak frequencies

  ip1 = mfs_peaks(abs(Uz1(1, :)));
  ip2 = mfs_peaks(abs(Uz2(1, :)));

  beam = mfs_back(beam, "freqresp", "disp", 1, f(ip1));
  beam = mfs_back(beam, "freqresp", "disp", 2, f(ip2));
  beam = mfs_results(beam, "freqresp", "element");

  mfs_export("freq.pos", "msh", beam, "freqresp", 
             "disp", "stress");

# Save results

  save -binary freqresp.bin beam
