# Example: Box Beam
#          Analysis of Excitation Functions
#
# --------------------------------------------------------

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

  pkg load signal

# Output device depending on operating system

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

# Pulse Function
# --------------

# Data

  T =    0.02;   % Duration of motion
  n =      50;   % Length of time series (multiples of T)
  dt = T / 40;   % Time step
  fc =    400;   % Cut-off frequency

# Function values

  t = 0 : dt : T;
  y = motion(t, T);

  y(2, :) *= T;
  y(3, :) *= T^2;

  figure(1, "position", [100, 100, 750, 750],
            "paperposition", [0, 0, 15, 15]);
  subplot(3, 1, 1)
     plot(t, y(1, :));
     grid;
     axis("labely");
     ylabel('\phi');
  subplot(3, 1, 2)
     plot(t, y(2, :));
     grid;
     axis("labely");
     ylabel('T d\phi/dt');
  subplot(3, 1, 3)
     plot(t, y(3, :));
     grid;
     ylabel('T^2 d^2\phi/dt^2');
     xlabel('t [s]');
  print(["motion", EXT], FORMAT);

# Fourier transform of acceleration

  fs = 1 / dt; fa = 0.5 * fs;
  N  = n * (length(t) - 1); 
  Ts = n * T;
  A  = fft(y(3, :), N) * dt;
  f  = (0 : N - 1) / Ts;

  figure(2, "position", [400, 100, 750, 500],
            "paperposition", [0, 0, 15, 10]);
  plot(f, abs(A));
  grid;
  xlim([0, fa]);
  xlabel('f [Hz]'); ylabel('T^2 * |\Phi_a(f)| [1/Hz]');
  print(["FT", EXT], FORMAT);

# Test of cut-off frequency

  yc = resample(y(3, :), fc, fa);
  tc = (0 : length(yc)-1) * dt * fa / fc;

  figure(3, "position", [600, 100, 750, 500],
            "paperposition", [0, 0, 15, 10]);
  plot(t, y(3, :), tc, yc);
  legend("orig.", "filtered");
  grid;
  xlabel('t [s]'); ylabel('T^2 * d^2\phi/dt^2');
  print(["cutoff", EXT], FORMAT);
