# Example: Frequency response analysis of a standard class glider
#          Plot the transfer functions for the gust response
#          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

# Get identifier of element at wing root

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

# Extract the transfer functions

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

# Plot the transfer functions

  figure(1, "position", [100, 200, 800, 500],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     semilogy(f, abs(U{1}));
     legend("WR", "Tail", "WT");
     title("v = 30 m/s");
     grid;
     ylabel("|U_z| [m/(m/s)]"); xlabel("f [Hz]");
  subplot(1, 2, 2);
     semilogy(f, abs(U{2}));
     title("v = 60 m/s");
     grid;
     xlabel("f [Hz]");
  print(["TF_U", EXT], FORMAT);

  figure(2, "position", [400, 200, 800, 500],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     semilogy(f, abs(UR{1}));
     legend("WR", "Tail", "WT", "location", "southwest");
     title("v = 30 m/s");
     grid;
     ylabel("|U_{Rz}| [m/(m/s)]"); xlabel("f [Hz]");
  subplot(1, 2, 2);
     semilogy(f, abs(UR{2}));
     title("v = 60 m/s");
     grid;
     xlabel("f [Hz]");
  print(["TF_UR", EXT], FORMAT);

  figure(3, "position", [700, 200, 800, 500],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     semilogy(f, abs(A{1}));
     legend("WR", "Tail", "WT", "location", "south");
     title("v = 30 m/s");
     grid;
     ylabel("|A_z| [(m/s^2)/(m/s)]"); xlabel("f [Hz]");
  subplot(1, 2, 2);
     semilogy(f, abs(A{2}));
     title("v = 60 m/s");
     grid;
     xlabel("f [Hz]");
  print(["TF_A", EXT], FORMAT);

  figure(4, "position", [1000, 200, 800, 500],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     semilogy(f, abs(My{1}));
     title("v = 30 m/s");
     grid;
     ylabel("|M_y| [Nm/(m/s)]"); xlabel("f [Hz]");
  subplot(1, 2, 2);
     semilogy(f, abs(My{2}));
     title("v = 60 m/s");
     grid;
     xlabel("f [Hz]");
  print(["TF_My", EXT], FORMAT);
