# Example: Plate
#          Modal frequency response analysis with subsequent
#          transformation to time domain
#
# -------------------------------------------------------------------

  colors = [1, 0, 0; 0, 1, 0; 0, 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.01;  % Duration of pulse
  delay = [0, 4, 2] * T;  % Delays

# Define simulation paramters

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

# Compute parameters for Fourier transform

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

# Get results from normal modes analysis

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

  load modes.bin

# Perform frequency response analysis

  f = 0 : df : fc; nfreq = length(f);
  for lc = 1 : 3
      plate = mfs_freqresp(plate, f, "loadcase", lc, "nband", 0);
  endfor

# Plot transfer functions

  rid = {"Resp_A", 3; "Resp_B", 3; "Resp_C", 3};

  HA = mfs_getresp(plate, "freqresp", "disp", rid, 1);
  HB = mfs_getresp(plate, "freqresp", "disp", rid, 2);
  HC = mfs_getresp(plate, "freqresp", "disp", rid, 3);

  figure(1, "position", [100, 200, 1000, 800],
            "paperposition", [0, 0, 15, 15]);
  subplot(3, 1, 1);
     semilogy(f, abs(HA));
     legend("Pnt. A", "Pnt. B", "Pnt. C",
            "location", "northeastoutside");
     grid;
     axis("labely");
     xlim([0, fc]);
     ylabel('|H_A| [mm/N]');
  subplot(3, 1, 2)
     semilogy(f, abs(HB));
     grid;
     axis("labely");
     xlim([0, fc]);
     ylabel('|H_B| [mm/N]');
  subplot(3, 1, 3)
     semilogy(f, abs(HC));
     grid;
     xlim([0, fc]);
     ylabel('|H_C| [mm/N]');
     xlabel('f [Hz]');
  print(["impact_tf", EXT], FORMAT);

# Get Fourier transforms of loads

  t  = 0 : dt : Ts;
  l  = pulse(t, [T, 0]);
  L  = fft(l) * dt;

  LA = L(1 : nfreq);
  LB = LA .* exp(-2 * pi * i * f * delay(2));
  LC = LA .* exp(-2 * pi * i * f * delay(3));

  figure(2, "position", [300, 200, 1000, 700],
            "paperposition", [0, 0, 15, 12]);
  plot(f, abs([LA; LB; LC]));
  legend('L_A', 'L_B', 'L_C');
  grid;
  xlim([0, fc]);
  xlabel('f [Hz]'); ylabel('|L| [N/Hz]');

# Compute Fourier transform of response

  W = HA .* LA + HB .* LB + HC .* LC;
  V = 2 * pi * i * f .* W / 1000;
  A = 2 * pi * i * f .* V / 9.81;

  figure(3, "position", [500, 200, 1000, 800],
            "paperposition", [0, 0, 15, 15]);
  subplot(3, 1, 1);
     semilogy(f(2 : end), abs(W(:, 2 : end)));
     legend("Pnt. A", "Pnt. B", "Pnt. C",
            "location", "northeastoutside");
     grid;
     axis("labely");
     xlim([0, fc]);
     ylabel('|W| [mm/Hz]');
  subplot(3, 1, 2)
     semilogy(f(2 : end), abs(V(:, 2 : end)));
     grid;
     axis("labely");
     xlim([0, fc]);
     ylabel('|V| [m/s/Hz]');
  subplot(3, 1, 3)
     semilogy(f(2 : end), abs(A(:, 2 : end)));
     grid;
     xlim([0, fc]);
     ylabel('|A| [g/Hz]');
     xlabel('f [Hz]');
  print(["impact_WVA", EXT], FORMAT);

# Perform transformation to time domain

  [w, t] = mfs_freq2time(W, df, dtb);
  [v, t] = mfs_freq2time(V, df, dtb);
  [a, t] = mfs_freq2time(A, df, dtb);

# Plot time histories

  figure(4, "position", [100, 100, 1000, 750],
            "paperposition", [0, 0, 15, 12]);
  plot(t, w);
  legend("Pnt. A", "Pnt. B", "Pnt. C", "autoupdate", "off");
  grid;
  xlim([0, Tp]);
  line([delay(2), delay(2)], ylim(), "color", "black");
  line([delay(3), delay(3)], ylim(), "color", "black");
  xlabel('t [s]'); ylabel('w [mm]');

  print(["impact_fresp_w", EXT], FORMAT);

  figure(5, "position", [300, 100, 1000, 750],
            "paperposition", [0, 0, 15, 12]);
  plot(t, v);
  legend("Pnt. A", "Pnt. B", "Pnt. C", "autoupdate", "off");
  grid;
  xlim([0, Tp]);
  line([delay(2), delay(2)], ylim(), "color", "black");
  line([delay(3), delay(3)], ylim(), "color", "black");
  xlabel('t [s]'); ylabel('v [m/s]');

  print(["impact_fresp_v", EXT], FORMAT);

  figure(6, "position", [500, 100, 1000, 750],
            "paperposition", [0, 0, 15, 12]);
  plot(t, a);
  legend("Pnt. A", "Pnt. B", "Pnt. C", "autoupdate", "off");
  grid;
  xlim([0, Tp]);
  line([delay(2), delay(2)], ylim(), "color", "black");
  line([delay(3), delay(3)], ylim(), "color", "black");
  xlabel('t [s]'); ylabel('a [g]');

  print(["impact_fresp_a", EXT], FORMAT);

# Save results

  dlmwrite("impact_fresp.csv", [t; w; v; a], "precision", 5);
