# Example: Frequency response analysis of a standard class glider
#          Fourier transform of the gusts
#          Fourier transforms are stored in file gust.bin.
# ----------------------------------------------------------------

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

  set(0, "defaultaxesfontsize", 10);

# Data of the gusts (m, s)

  L  = 30;        % Gust length
  w0 = 1;         % Gust velocity amplitude
  v  = [30, 60];  % Flight velocities
  dt = 0.01;      % Time step

# Initialization

  T    = 5 * L / v(1);   % Duration
  N    = T / dt;         % Length of time series
  fmax = 1 / (2 * dt);   % Nyquist frequency
  df   = 1 / T;          % Frequency increment
  f    = 0 : df : fmax;  % Frequencies
  nf   = length(f);      % Number of frequencies

  W1 = zeros(2, N); W2 = zeros(2, N);

# Loop over velocities and compute Fourier transforms

  for n = 1 : 2

      text{n} = sprintf("v = %2.0f m/s", v(n));
      T0      = L / v(n); % Time to pass the gust

#   Time series

      t0 = 0 : (dt / T0) : 1;
      w1 = w0 * sin(pi * t0).^2;
      w2 = w1 - w0 * sin(2 * pi * t0).^2;

#   Fourier transforms

      W1(n, :) = dt * fft(w1, N); 
      W2(n, :) = dt * fft(w2, N);

  endfor

# Plot data and results

  figure(1, "position", [100, 100, 800, 400],
            "paperposition", [0, 0, 17, 8]);
  plot(t0, w1, t0, w2);
  legend("Gust 1", "Gust 2", "location", "south");
  grid;
  xlabel("t/T_0"); ylabel("w/w_0");
  print(["gust", EXT], FORMAT);

  figure(2, "position", [100, 100, 800, 400],
            "paperposition", [0, 0, 17, 10]);
  subplot(2, 1, 1);
     plot(f, abs(W1(1, 1 : nf)), f, abs(W2(1, 1 : nf)));
     legend("Gust 1", "Gust 2");
     title(text{1});
     grid; 
     xlim([0, 10]);
     ylabel("|W| [m/s/Hz]");
  subplot(2, 1, 2);
     plot(f, abs(W1(2, 1 : nf)), f, abs(W2(2, 1 : nf)));
     title(text{2});
     grid; 
     xlim([0, 10]);
     ylabel("|W| [m/s/Hz]");
     xlabel("f [Hz]");
  print(["gust_ft", EXT], FORMAT);

# Save Fourier transforms for later use

  W{1} = W1; W{2} = W2;

  save -binary gust.bin W df v
