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

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

  set(0, "defaultaxesfontsize", 10);

  pkg load signal

# Cut-off frequency

  fc = 10;

# Data of the gust (m, s)

  L  = 30;        % Gust length
  w0 = 1;         % Gust velocity amplitude
  v  = 60;        % Flight velocity
  dt = 0.01;      % Time step
  P  = 5;         % Upsampling factor

# Initialization

  fmax = 1 / (2 * dt);  % Nyquist frequency
  T0   = L / v;         % Time to pass the gust

# Original time series

  t = 0 : (dt / T0) : 1;
  w = w0 * (sin(pi * t).^2 - sin(2 * pi * t).^2);

# Filtered time series

  wc = resample(w, fc, fmax);
  fs = 2 * fc; fsT0 = fs * T0;            % Sampling rate
  tc = (0 : length(wc) - 1) / fsT0;       % t/T0
  wp = resample(wc, P, 1);
  tp = (0 : length(wp) - 1) / (P * fsT0); % t/T0

# Plot results

  figure(1, "position", [100, 100, 800, 400],
            "paperposition", [0, 0, 17, 10]);
  plot(t, w, tc, wc, tp, wp);
  legend("Original", "Filtered", "Upsampled", 
         "location", "south");
  grid; xlim([0, 1]);
  xlabel("t/T_0"); ylabel("w [m/s]");
  print(["gust_check", EXT], FORMAT);
