# Example: Frequency response analysis of a standard class glider
#          Gust response
#          Files needed:  trf.bin, gust.bin
# ---------------------------------------------------------------

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

  set(0, "defaultaxesfontsize", 10);

# Define response degrees of freedom and wing root element

  rid   = [5, 3; 15, 3; 115, 3];   % Node, Dof

# Define parameters of transient response

  dt = 0.005;   % Time step
  te =     3;   % Duration (for plotting)

# Get solid component with frequency response

  load trf.bin
  gliders = mfs_extract(glider, "solid");
  clear glider

# Get identifier of element at wing root

  eltids = mfs_getset(gliders, "eset", "left_spar");
  eltid  = eltids(1);

# Get frequencies, include zero frequency

  f  = [0, mfs_getresp(gliders, "freqresp", "freq")];
  nf = length(f);

# Get Fourier transforms of the gusts

  load ../loads/gust.bin

# Compute responses

  for lc = 1 : 2        % Velocity loop

#   Get transfer functions (values at f = 0 are zero)

      TU(:, 2 : nf)  = mfs_getresp(gliders, "freqresp", 
                                   "disp", rid, lc);
      TUR(:, 2 : nf) = mfs_getresp(gliders, "freqresp", 
                                   "reldisp", rid, lc);
      TA(:, 2 : nf)  = mfs_getresp(gliders, "freqresp", 
                                   "acce", rid, lc);
      R              = mfs_getresp(gliders, "freqresp",
                                   "resultant", eltid, lc);
      TM(2 : nf)     = R{1}.My;

      for k = 1 : 2    % Gust loop

          titles{lc, k} = sprintf("Gust %1d: v = %2d m/s",
                                  k, v(lc));

#      Compute Fourier transforms of responses

          WG = W{k}(lc, 1 : nf);
          U{lc, k} = TU .* WG; UR{lc, k} = TUR .* WG;
	  A{lc, k} = TA .* WG; M{lc, k}  = TM .* WG;

#      Compute transient responses

          [u{lc, k}, t] = mfs_freq2time(U{lc, k}, df, dt);
          uR{lc, k}     = mfs_freq2time(UR{lc, k}, df, dt);
          a{lc, k}      = mfs_freq2time(A{lc, k}, df, dt);
          m{lc, k}      = mfs_freq2time(M{lc, k}, df, dt);

          if (k == 1)
             uc{lc} = mfs_freq2time(U{lc, k}, df, dt, 2);
          endif

      endfor           % Gust loop

  endfor               % Velocity loop

# Save transient displacements for later comparison

  save -binary response.bin t u

# Plot Fourier transforms of responses

  % Displacements

  figure(1, "position", [100, 100, 800, 800],
            "paperposition", [0, 0, 17, 15]);
  subplot(2, 2, 1)
     plot(f, abs(U{1, 1}));
     legend("WR", "Tail", "WT");
     title(titles{1, 1});
     grid;
     axis("labely");
     ylabel("|U_z| [m/Hz]");
  subplot(2, 2, 2)
     plot(f, abs(U{2, 1}));
     title(titles{2, 1});
     grid;
     axis("labely");
  subplot(2, 2, 3)
     plot(f, abs(U{1, 2}));
     title(titles{1, 2});
     grid;
     ylabel("|U_z| [m/Hz]"); xlabel("f [Hz]");
  subplot(2, 2, 4)
     plot(f, abs(U{2, 2}));
     title(titles{2, 2});
     grid;
     xlabel("f [Hz]");
  print(["FT_U", EXT], FORMAT);

  % Relative displacements

  figure(2, "position", [300, 100, 800, 800],
            "paperposition", [0, 0, 17, 15]);
  subplot(2, 2, 1)
     plot(f, abs(UR{1, 1}));
     legend("WR", "Tail", "WT");
     title(titles{1, 1});
     grid;
     axis("labely");
     ylabel("|U_{Rz}| [m/Hz]");
  subplot(2, 2, 2)
     plot(f, abs(UR{2, 1}));
     title(titles{2, 1});
     grid;
     axis("labely");
  subplot(2, 2, 3)
     plot(f, abs(UR{1, 2}));
     title(titles{1, 2});
     grid;
     ylabel("|U_{Rz}| [m/Hz]"); xlabel("f [Hz]");
  subplot(2, 2, 4)
     plot(f, abs(UR{2, 2}));
     title(titles{2, 2});
     grid;
     xlabel("f [Hz]");
  print(["FT_UR", EXT], FORMAT);

  % Accelerations

  figure(3, "position", [500, 100, 800, 800],
            "paperposition", [0, 0, 17, 15]);
  subplot(2, 2, 1)
     plot(f, abs(A{1, 1}));
     legend("WR", "Tail", "WT");
     title(titles{1, 1});
     grid;
     axis("labely");
     ylabel("|A_z| [m/s^2/Hz]");
  subplot(2, 2, 2)
     plot(f, abs(A{2, 1}));
     title(titles{2, 1});
     grid;
     axis("labely");
  subplot(2, 2, 3)
     plot(f, abs(A{1, 2}));
     title(titles{1, 2});
     grid;
     ylabel("|A_z| [m/s^2/Hz]"); xlabel("f [Hz]");
  subplot(2, 2, 4)
     plot(f, abs(A{2, 2}));
     title(titles{2, 2});
     grid;
     xlabel("f [Hz]");
  print(["FT_A", EXT], FORMAT);

  % Bending moment

  figure(4, "position", [700, 100, 800, 500],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1)
     plot(f, abs(M{1, 1}), f, abs(M{1, 2}));
     legend("Gust 1", "Gust 2");
     title(sprintf("v = %2d m/s", v(1)));
     grid;
     ylabel("|M_y| [Nm/Hz]"); xlabel("f [Hz]");
  subplot(1, 2, 2)
     plot(f, abs(M{2, 1}), f, abs(M{2, 2}));
     title(sprintf("v = %2d m/s", v(2)));
     grid;
     xlabel("f [Hz]");
  print(["FT_M", EXT], FORMAT);

# Plot transient responses

  % Displacements without baseline correction

  figure(5, "position", [100, 200, 800, 800],
            "paperposition", [0, 0, 17, 15]);
  subplot(2, 2, 1)
     plot(t, u{1, 1});
     legend("WR", "Tail", "WT");
     title(titles{1, 1});
     grid; xlim([0, te]);
     axis("labely");
     ylabel("u_z [m]");
  subplot(2, 2, 2)
     plot(t, u{2, 1});
     title(titles{2, 1});
     grid; xlim([0, te]);
     axis("labely");
  subplot(2, 2, 3)
     plot(t, u{1, 2});
     title(titles{1, 2});
     grid; xlim([0, te]);
     ylabel("u_z [m]"); xlabel("t [s]");
  subplot(2, 2, 4)
     plot(t, u{2, 2});
     title(titles{2, 2});
     grid; xlim([0, te]);
     xlabel("t [s]");
  print(["u", EXT], FORMAT);

  % Displacements with baseline correction

  figure(6, "position", [300, 200, 800, 500],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1)
     plot(t, uc{1});
     legend("WR", "Tail", "WT", "location", "east");
     title(sprintf("v = %2d m/s", v(1)));
     grid; xlim([0, te]);
     ylabel("u_z [m]"); xlabel("t [s]");
  subplot(1, 2, 2)
     plot(t, uc{2});
     title(sprintf("v = %2d m/s", v(2)));
     grid; xlim([0, te]);
     xlabel("t [s]");
  print(["uc", EXT], FORMAT);

  % Relative displacements

  figure(7, "position", [500, 200, 800, 800],
            "paperposition", [0, 0, 17, 15]);
  subplot(2, 2, 1)
     plot(t, uR{1, 1});
     legend("WR", "Tail", "WT");
     title(titles{1, 1});
     grid; xlim([0, te]);
     axis("labely");
     ylabel("u_{Rz} [m]");
  subplot(2, 2, 2)
     plot(t, uR{2, 1});
     title(titles{2, 1});
     grid; xlim([0, te]);
     axis("labely");
  subplot(2, 2, 3)
     plot(t, uR{1, 2});
     title(titles{1, 2});
     grid; xlim([0, te]);
     ylabel("u_{Rz} [m]"); xlabel("t [s]");
  subplot(2, 2, 4)
     plot(t, uR{2, 2});
     title(titles{2, 2});
     grid; xlim([0, te]);
     xlabel("t [s]");
  print(["uR", EXT], FORMAT);

  % Accelerations

  figure(8, "position", [700, 200, 800, 800],
            "paperposition", [0, 0, 17, 15]);
  subplot(2, 2, 1)
     plot(t, a{1, 1});
     legend("WR", "Tail", "WT");
     title(titles{1, 1});
     grid; xlim([0, te]);
     axis("labely");
     ylabel("a_z [m/s^2]");
  subplot(2, 2, 2)
     plot(t, a{2, 1});
     title(titles{2, 1});
     grid; xlim([0, te]);
     axis("labely");
  subplot(2, 2, 3)
     plot(t, a{1, 2});
     title(titles{1, 2});
     grid; xlim([0, te]);
     ylabel("a_z [m/s^2]"); xlabel("t [s]");
  subplot(2, 2, 4)
     plot(t, a{2, 2});
     title(titles{2, 2});
     grid; xlim([0, te]);
     xlabel("t [s]");
  print(["a", EXT], FORMAT);

  % Bending moment

  figure(9, "position", [900, 200, 800, 500],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1)
     plot(t, m{1, 1}/1000, t, m{1, 2}/1000);
     legend("Gust 1", "Gust 2");
     title(sprintf("v = %2d m/s", v(1)));
     grid; xlim([0, te]);
     ylabel("M_y [kNm]"); xlabel("t [s]");
  subplot(1, 2, 2)
     plot(t, m{2, 1}/1000, t, m{2, 2}/1000);
     title(sprintf("v = %2d m/s", v(2)));
     grid; xlim([0, te]);
     xlabel("t [s]");
  print(["My", EXT], FORMAT);
