# Example: Trim analysis of a standard class glider
#
# c) Flexible Trim
#
# Manoeuvres:
#
#  1. Straight level flight
#  2. Truly banked turn (bank angle = 30°)
#  3. Turn with zero side slip (bank angle = 30°)
#  4. Sudden aileron deflection (deflection angle = 1°)
#
# --------------------------------------------------------------------

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

  set(0, "defaultaxesfontsize", 10);

  fid  = fopen("flexible.res", "wt"); 

# ====================================================================
# Model Definition
# ====================================================================

  load solid.bin
  load aero.bin

  model = struct("type", "aeroelastic",
                 "solid", gliders, "aero", glidera);

  clear gliders; clear glidera

# Splines
# -------

  nbw = 15;  % Number of wing spline breaks
  nbv =  3;  % Number of vertical stabilizer spline breaks
  nbh =  3;  % Number of horizontal stabilizer spline breaks

  data = struct("nbreaks", {nbw, nbv, nbh});

# Right and left wing splines

  splines(1 : 2) = struct("id",    {10, 20},
                          "type",  "tb",
                          "lsid",  {10 : 12, 20 : 22},
                          "data",  data(1),
                          "nodes", {"right_wing", "left_wing"});

# Vertical stabilizer

  splines(3) = struct("id",    30,
                      "type",  "tb",
                      "lsid",  [31, 32],
                      "data",  data(2),
                      "nodes", "vertical");

# Horizontal stabilizer

  splines(4) = struct("id",   40,
                      "type", "tb",
                      "lsid", [41, 42, 43, 44],
                      "data", data(3),
                      "nodes", "horizontal");

  model.splines = splines;

# ====================================================================
# Analysis
# ====================================================================

# Create aeroelastic component

  glider = mfs_new(fid, model);

# Compute splines

  glider = mfs_splines(glider);

# Restrained trim analysis
# ------------------------

# Perform the trim analysis

  [rtrims, rtrima] = mfs_trim(glider, "restrained");
  mfs_print(fid, rtrima, "trim", "params");

# Compute and export results

  rtrima = mfs_results(rtrima, "trim", "panel");
  mfs_export("flexible_ra.pos", "msh", rtrima,
             "trim", "pressure", "disp");
  mfs_export("flexible_rs.dsp", "msh", rtrims, "trim", "disp");

# Get load resultants with respect to centre of mass

  load solid_mass.bin
  [Fr, Mr] = mfs_getresp(rtrima, "trim", "aeload", mp.cm);

# Unrestrained trim analysis
# --------------------------

# Perform the trim analysis

  [utrims, utrima] = mfs_trim(glider, "unrestrained");
  mfs_print(fid, utrima, "trim", "params");

# Compute and export results

  utrima = mfs_results(utrima, "trim", "panel");
  mfs_export("flexible_ua.pos", "msh", utrima,
             "trim", "pressure", "disp");
  mfs_export("flexible_us.dsp", "msh", utrims, "trim", "disp");

# Get load resultants with respect to centre of mass

  [Fu, Mu] = mfs_getresp(utrima, "trim", "aeload", mp.cm);

# Print results
# -------------

  nconf = columns(Fr);

  fprintf(fid, "\nLoad resultants with respect to center of mass:\n");

  fprintf(fid, "\n Restrained analysis:\n");
  for n = 1 : nconf
      fprintf(fid, "\n  Configuration %2d:\n", n)
      fprintf(fid, "     F = [%10.3e, %10.3e, %10.3e] kN\n",
              Fr(:, n) / 1000);
      fprintf(fid, "     M = [%10.3e, %10.3e, %10.3e] kNm\n",
              Mr(:, n) * 1e-3);
      mname{n} = ["Man. ", num2str(n)];
  endfor

  fprintf(fid, "\n Unrestrained analysis:\n");
  for n = 1 : nconf
      fprintf(fid, "\n  Configuration %2d:\n", n)
      fprintf(fid, "     F = [%10.3e, %10.3e, %10.3e] kN\n",
              Fu(:, n) / 1000);
      fprintf(fid, "     M = [%10.3e, %10.3e, %10.3e] kNm\n",
              Mu(:, n) * 1e-3);
  endfor

# Compare pressure in wing sections
# ---------------------------------

  load rigid.bin

  [x1, p1fr, y(1)] = mfs_xydata(rtrima, "trim", "pressure",
                              [21, 22], ycolo(1));
  [x2, p2fr, y(2)] = mfs_xydata(rtrima, "trim", "pressure",
                              20, ycoli(1));
  [x3, p3fr, y(3)] = mfs_xydata(rtrima, "trim", "pressure",
                              10, ycoli(2));
  [x4, p4fr, y(4)] = mfs_xydata(rtrima, "trim", "pressure",
                              [11, 12], ycolo(2));

  [x1, p1fu, y(1)] = mfs_xydata(utrima, "trim", "pressure",
                              [21, 22], ycolo(1));
  [x2, p2fu, y(2)] = mfs_xydata(utrima, "trim", "pressure",
                              20, ycoli(1));
  [x3, p3fu, y(3)] = mfs_xydata(utrima, "trim", "pressure",
                              10, ycoli(2));
  [x4, p4fu, y(4)] = mfs_xydata(utrima, "trim", "pressure",
                              [11, 12], ycolo(2));

  for n = 1 : 4
      ytext{n} = sprintf("y = %6.3f m", y(n));
  endfor

  nfig = 1;

  for n = 1 : nconf

      figure(nfig++, "position", [100 + 100 * n, 200, 800, 500],
                     "paperposition", [0, 0, 17, 12]);
      subplot(2, 2, 1);
         plot(x1, [p1r(:, n), p1fr(:, n), p1fu(:, n)] * 1e-3);
         legend("Rigid", "Restrained", "Unrestrained",
                "location", "south");
         title(ytext{1});
         grid;
         ylim([0, 0.4]); 
         ylabel('\Delta p [kPa]');
      subplot(2, 2, 2);
         plot(x4, [p4r(:, n), p4fr(:, n), p4fu(:, n)] * 1e-3);
         title(ytext{4});
         grid;
         ylim([-0.2, 0.4]);
      subplot(2, 2, 3);
         plot(x2, [p2r(:, n), p2fr(:, n), p2fu(:, n)] * 1e-3);
         title(ytext{2});
         grid;
         ylim([0.1, 0.5]);
         xlabel('x [m]'); ylabel('\Delta p [kPa]');
      subplot(2, 2, 4);
         plot(x3, [p3r(:, n), p3fr(:, n), p3fu(:, n)] * 1e-3);
         title(ytext{3});
         grid;
         ylim([0.1, 0.5]);
         xlabel('x [m]');
      print([sprintf("p%d", n), EXT], FORMAT);

  endfor

  dp1 = p1r - p1fr; dp2 = p2r - p2fr;
  dp3 = p3r - p3fr; dp4 = p4r - p4fr;

  figure(nfig++, "position", [200, 400, 700, 500],
                 "paperposition", [0, 0, 17, 12]);
  subplot(2, 2, 1);
     plot(x1, dp1);
     legend(mname, "location", "northeast");
     title(ytext{1});
     grid;
     ylabel('\Delta p_r - \Delta p_f [Pa]');
  subplot(2, 2, 2);
     plot(x4, dp4);
     title(ytext{4});
     grid;
  subplot(2, 2, 3);
     plot(x2, dp2);
     title(ytext{2});
     grid;
     xlabel('x [m]');
     ylabel('\Delta p_r - \Delta p_f [Pa]');
  subplot(2, 2, 4);
     plot(x3, dp3);
     title(ytext{3});
     grid;
     xlabel('x [m]');
  print(["dp", EXT], FORMAT);

# Shear force and bending moment in main spar
# -------------------------------------------

# Get stress resultants and convert results to structure array

  utrims = mfs_results(utrims, "trim", "element");

  Rl = mfs_getresp(utrims, "trim", "resultant",
                   "left_spar", 1 : nconf);
  rl = cell2mat(Rl);

  Rr = mfs_getresp(utrims, "trim", "resultant",
                   "right_spar", 1 : nconf);
  rr = cell2mat(Rr);

# Create array with y-coordinates along spar

  nell = length(Rl); nelr = length(Rr);

  coorl = reshape([rl.coor], 3, nell);
  coorr = reshape([rr.coor], 3, nelr);
  y     = [flip(coorl(2, :)), coorr(2, :)];

# Create arrays with stress resultants Qz and My

  Myl = reshape([rl.My], nconf, nell);
  Myr = reshape([rr.My], nconf, nelr);
  My  = [flip(Myl, 2), Myr];

  Qzl = reshape([rl.Qz], nconf, nell);
  Qzr = reshape([rr.Qz], nconf, nelr);
  Qz  = [flip(Qzl, 2), Qzr];

# Plot the results

  figure(nfig++, "position", [300, 400, 700, 500],
                 "paperposition", [0, 0, 17, 12]);
  subplot(2, 1, 1)
     plot(y, Qz * 1e-3);
     legend(mname);
     grid;
     xlim([-7.5, 7.5]);
     axis("labely");
     ylabel('Q_z [kN]');
  subplot(2, 1, 2)
     plot(y, My * 1e-3);
     grid;
     xlim([-7.5, 7.5]);
     xlabel('y [m]');
     ylabel('M_y [kNm]');
  print(["QM", EXT], FORMAT);

# Compare radius of turn
# ----------------------

# Truly banked turn

  R2 = v^2 / (g * tand(gamma));

# Turn with zero side slip

  tp(1) = mfs_getresp(rtrima, "trim", "params", 3);
  tp(2) = mfs_getresp(utrima, "trim", "params", 3);

  R3 = v^2 ./ (sind(gamma) * [tp.az] - cosd(gamma) * [tp.ay]);

  fprintf(fid, "\nRadius of turn:\n");
  fprintf(fid, "  Truly banked                 : %6.2f m\n",
               R2);
  fprintf(fid, "  Zero side slip (restrained)  : %6.2f m\n",
               R3(1));
  fprintf(fid, "  Zero side slip (unrestrained): %6.2f m\n",
               R3(2));

# Combine results
# ---------------

  mfs_merge("solid.msh", "aero.msh", "glider.msh", "msh");
  mfs_merge("flexible_rs.dsp", "flexible_ra.pos", 
            "glider_r.dsp", "msh");
  mfs_merge("flexible_us.dsp", "flexible_ua.pos", 
            "glider_u.dsp", "msh");

  fclose(fid);
