# Example: Frequency response analysis of a standard class glider
#          Response to elevator input
#          Files needed:  trf.bin, elevator.bin
# ---------------------------------------------------------------

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

  set(0, "defaultaxesfontsize", 10);

# Define response degrees of freedom

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

# Define parameters of transient response

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

# Get solid component with frequency response

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

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

# Get Fourier transform of elevator input

  load ../loads/elevator.bin

# Inverse Fourier transform for second gust and second velocity

  spectrum = struct("df", df, "lc", 1, "spec", ETA);

  gliders = mfs_freq2time(gliders, spectrum, 1, dt, 2);

# Get results

  t  = mfs_getresp(gliders, "transresp", "time");
  u  = mfs_getresp(gliders, "transresp", "disp", rid);
  uR = mfs_getresp(gliders, "transresp", "reldisp", rid);
  a  = mfs_getresp(gliders, "transresp", "acce", rid);
  R  = mfs_getresp(gliders, "transresp", "resultant", eltid(1));
  My = R{1}.My;

# Plot results

  figure(1, "position", [100, 100, 800, 400],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     plot(t, u);
     legend("WR", "Tail", "WT", "location", "east");
     title("Absolute Displacements");
     grid;
     xlim([0, te]);
     ylabel('u_z [m]'); xlabel("t [s]");
  subplot(1, 2, 2);
     plot(t, uR * 100);
     title("Relative Displacements");
     grid;
     xlim([0, te]);
     ylabel('u_{Rz} [cm]'); xlabel("t [s]");
  print(["u", EXT], FORMAT);

  figure(2, "position", [300, 100, 800, 400],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     plot(t, a);
     legend("WR", "Tail", "WT");
     title("Accelerations");
     grid;
     xlim([0, te]);
     ylabel('a_z [m/s^2]'); xlabel("t [s]");
  subplot(1, 2, 2);
     plot(t, My / 1000);
     title("Bending Moment");
     grid;
     xlim([0, te]);
     ylabel('M_y [kNm]'); xlabel("t [s]");
  print(["aM", EXT], FORMAT);

# Perform backtransformation and export displacements

  tb = t(find(t <= te));
  gliders = mfs_back(gliders, "transresp", "disp", 1, tb);
  mfs_export("elevator.dsp", "msh", gliders,
             "transresp", "disp");
