# Example: Frequency response analysis of a standard class glider
#          Plot the transfer functions for the response to an
#          elevator input
#          Files needed:  trf.bin
# ---------------------------------------------------------------

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

  set(0, "defaultaxesfontsize", 10);

# Define response degrees of freedom

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

# Load the aeroelastic component

  load trf.bin

# Extract the solid component

  gliders = mfs_extract(glider, "solid");
  clear glider

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

# Extract the transfer functions

  f  = mfs_getresp(gliders, "freqresp", "freq");
  U  = mfs_getresp(gliders, "freqresp", "disp", rid);
  UR = mfs_getresp(gliders, "freqresp", "reldisp", rid);
  A  = mfs_getresp(gliders, "freqresp", "acce", rid);
  R  = mfs_getresp(gliders, "freqresp", "resultant", eltid(1));
  My = R{1}.My;

# Plot the transfer functions

  figure(1, "position", [100, 100, 800, 400],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     semilogy(f, abs(U));
     legend("WR", "Tail", "WT");
     title("Absolute Displacements");
     grid;
     ylabel('|U_z| [m/\deg]'); xlabel("f [Hz]");
  subplot(1, 2, 2);
     semilogy(f, abs(UR));
     title("Relative Displacements");
     grid;
     ylabel('|U_{Rz}| [m/\deg]'); xlabel("f [Hz]");
  print(["TF_U", EXT], FORMAT);

  figure(2, "position", [300, 100, 800, 400],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     semilogy(f, abs(A));
     legend("WR", "Tail", "WT");
     title("Accelerations");
     grid;
     ylabel('|A_z| [(m/s^2)/\deg]'); xlabel("f [Hz]");
  subplot(1, 2, 2);
     semilogy(f, abs(My));
     title("Bending Moment");
     grid;
     ylabel('|M_y| [Nm/\deg]'); xlabel("f [Hz]");
  print(["TF_AM", EXT], FORMAT);
