# Example: Box Beam
#          Modal frequency response analysis with subsequent
#          transformation to the time domain
#
# -------------------------------------------------------------------

  colors = [1, 0, 0; 0, 1, 0; 0, 0, 1; 1, 0, 1];
  set(0, "defaultaxescolororder", colors);

# Output device depending on operating system

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

# Define excitation parameters

  T   = 0.02;  % Duration of prescribed motion

# Define simulation parameters

  fc  =  400;  % Cut-off frequency
  df  =    1;  % Frequency resolution
  dtb = 1e-4;  % Time step for backtransformation
  Tp  = 0.07;  % Time to plot

# Compute parameters for Fourier transform

  dt = 1 / (2 * fc);  % Time step
  Ts = 1 / df;        % Duration of time series 

# Get component with results

  if (! isfile("modes.bin"))
     error("File modes.m must be executed before\n");
  endif

  load modes.bin

# Perform frequency response analysis

  f = 0 : df : fc; nfreq = length(f);
  beam = mfs_freqresp(beam, f, "nband", 0);

# Plot transfer functions

  rid  = {"Point_A", 3; "Point_B", 3};

  H = mfs_getresp(beam, "freqresp", "acce", rid) / 9810;

  figure(1, "position", [100, 100, 600, 400],
            "paperposition", [0, 0, 15, 10]);
  semilogy(f(2 : end), abs(H(:, 2 : end)));
  legend("Point A", "Point B", "location", "southeast");
  grid;
  xlabel('f [Hz]'); ylabel('H [g/mm]');
  print(["fresp_tf", EXT], FORMAT);

# Get Fourier transform of prescribed displacements

  t   = 0 : dt : Ts;
  uva = motion(t, T);
  Az  = fft(uva(3, :)) * dt;
  wq  = (2 * pi * i * f).^2;
  Uz(2 : nfreq) = Az(2 : nfreq) ./ wq(2 : nfreq);

  figure(2, "position", [300, 100, 750, 500]);
  plot(f, abs(Uz));
  grid; xlim([0, fc]);
  xlabel('f [Hz]'); ylabel('|Uz| [mm/Hz]');

# Perform transformation to time domain

  spectrum = struct("df", df, "lc", 1, "spec", Uz);
  beam = mfs_freq2time(beam, spectrum, 3, dtb, 2);

# Plot time histories

  t = mfs_getresp(beam, "transresp", "time", 3);
  u = mfs_getresp(beam, "transresp", "disp", rid, 3);
  v = mfs_getresp(beam, "transresp", "velo", rid, 3) / 1000;
  a = mfs_getresp(beam, "transresp", "acce", rid, 3) / 9810;

  figure(4, "position", [500, 200, 1000, 800],
            "paperposition", [0, 0, 15, 15]);
  subplot(3, 1, 1)
     plot(t, u);
     legend("Point A", "Point B", "location", "southeast");
     axis("labely");
     grid; xlim([0, Tp]);
     ylabel('u [mm]');
  subplot(3, 1, 2)
     plot(t, v);
     axis("labely");
     grid; xlim([0, Tp]);
     ylabel('v [m/s]');
  subplot(3, 1, 3)
     plot(t, a);
     grid; xlim([0, Tp]);
     ylabel('a [g]'); xlabel('t [s]');
  print(["fresp", EXT], FORMAT);

# Save results

  dlmwrite("fresp.csv", [t; u; v; a], "precision", 5);
