# Example: Frequency response analysis of a standard class glider
#          Check cut-off frequency of elevator input
# ---------------------------------------------------------------

  addpath("../../../..");
  [EXT, FORMAT] = iniplot();

  set(0, "defaultaxesfontsize", 10);

  pkg load signal

# Cut-off frequency

  fc = 10;

# Data of elevator input (degrees, s)

  eta0 =    5;    % Amplitude of elevator angle
  T0   =  1.0;    % Duration of elevator motion
  dt   = 0.01;    % Time step
  P    =    2;    % Upsampling factor

# Initialization

  fmax = 1 / (2 * dt);  % Nyquist frequency

# Original time series

  t   = 0 : (dt / T0) : 1;
  eta = eta0 * sin(2 * pi * t).^3;

# Filtered time series

  etac = resample(eta, fc, fmax);
  fs   = 2 * fc; fsT0 = fs * T0;
  tc   = (0 : length(etac) - 1) / fs;
  etap = resample(etac, P, 1);
  tp   = (0 : length(etap) - 1) / (P * fs);

# Plot results

  figure(1, "position", [100, 100, 800, 400],
            "paperposition", [0, 0, 17, 7]);
  plot(t,  eta, tc, etac, tp, etap);
  legend("Original", "Filtered", "Upsampled",
         "location", "northeast");
  grid;
  xlim([0, 1]);
  xlabel("t/T_0"); ylabel('\eta_E [\deg]');
  print(["elevator_check", EXT], FORMAT);
