# Example: Trim analysis of a standard class glider
#
# 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("trim.res", "wt"); 

# Data (kg, m, s)
# ---------------

  v     =   30;   % Flight velocity
  rho   = 1.21;   % Mass density of air
  g     = 9.81;   % Gravity acceleration
  gamma =   30;   % Bank angle in degrees
  eta   =    1;   % Aileron deflection in degrees
  ds    =  1.0;   % Distance of yaw string from center of mass

# Model Definition
# ----------------

  model = glider_aero();

# Configurations

  qdyn = 0.5 * rho * v^2;

  % Straight level flight

  config(1).name  = "Manoeuvre 1: Straight level flight";
  config(1).qdyn  = qdyn;
  config(1).ay    = 0;
  config(1).az    = g;
  config(1).pacce = 0;
  config(1).yacce = 0;
  config(1).racce = 0;
  config(1).pitch = 0;
  config(1).yaw   = 0;
  config(1).roll  = 0;

  % Truly banked turn

  tg = tand(gamma);
  cs = 1 / sqrt(1 + tg^2); sn = tg * cs;
  w  = tg * g / v;

  config(2).name  = "Manoeuvre 2: Truly banked turn";
  config(2).qdyn  = qdyn;
  config(2).ay    = 0;
  config(2).az    = g / cs;
  config(2).pacce = 0;
  config(2).yacce = 0;
  config(2).racce = 0;
  config(2).pitch = w * sn / v;
  config(2).yaw   = w * cs / v;
  config(2).roll  = 0;

  % Turn with zero side slip

  config(3).name  = "Manoeuvre 3: Turn with zero side slip";
  config(3).qdyn  = qdyn;
  config(3).pacce = 0;
  config(3).yacce = 0;
  config(3).racce = 0;
  config(3).roll  = 0;

  ss = sn * sn; cc = cs * cs; sncs = sn * cs; vv = v * v;

  lincon{1} = struct("ay", sn, "az", cs, "rhs", g);
  lincon{2} = struct("ay", sncs, "az", -ss, "pitch", vv);
  lincon{3} = struct("ay", cc, "az", -sncs, "yaw", vv);
  lincon{4} = struct("beta", 1, "yaw", ds);

  config(3).lincon = lincon;

  % Sudden aileron deflection

  config(4).name  = "Manoeuvre 4: Sudden aileron deflection";
  config(4).qdyn  = qdyn;
  config(4).ay    = 0;
  config(4).az    = g;
  config(4).pacce = 0;
  config(4).yacce = 0;
  config(4).pitch = 0;
  config(4).yaw   = 0;
  config(4).roll  = 0;
  config(4).aileron = eta;

  model.config = config;

# Analysis
# --------

# Create and export component

  glider = mfs_new(fid, model);
  mfs_export("glider.msh", "msh", glider, "mesh",
             "mesh", "normals");

# Perform the trim analysis

  glider = mfs_trim(glider);
  mfs_print(fid, glider, "trim", "params");

# Compute and export results

  glider = mfs_results(glider, "trim", "panel");
  mfs_export("trim.pos", "msh", glider, "trim", "pressure");

# Get load resultants

  [F, M] = mfs_getresp(glider, "trim", "aeload", model.mass.cm);
  [n, nconf] = size(F);

  fprintf(fid, 
          "\nLoad resultants with respect to center of mass:\n");
  for k = 1 : nconf
      fprintf(fid, "\n  Configuration %2d:\n", k)
      fprintf(fid, "     F = [%10.3e, %10.3e, %10.3e] kN\n",
              F(:, k) / 1000);
      fprintf(fid, "     M = [%10.3e, %10.3e, %10.3e] kNm\n",
              M(:, k) * 1e-3);
  endfor

# Plot pressure in some wing sections

  nyi   = model.ls(1).ny;   nyo  = model.ls(2).ny;
  ycmi  = floor(0.5 * nyi); ycmo = floor(0.5 * nyo);
  ycoli = [ycmi, nyi + 1 - ycmi];
  ycolo = [ycmo, nyo + 1 - ycmo];

  [x1, p1, y(1)] = mfs_xydata(glider, "trim", "pressure",
                              [21, 22], ycolo(1));
  [x2, p2, y(2)] = mfs_xydata(glider, "trim", "pressure",
                              20, ycoli(1));
  [x3, p3, y(3)] = mfs_xydata(glider, "trim", "pressure",
                              10, ycoli(2));
  [x4, p4, y(4)] = mfs_xydata(glider, "trim", "pressure",
                              [11, 12], ycolo(2));

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

  figure(1, "position", [200, 200, 1000, 1000],
            "paperposition", [0, 0, 15, 15]);
  subplot(2, 2, 1);
     plot(x1, p1 * 1e-3);
     title(tittxt{1});
     grid;
     ylim([0, 0.4]);
     ylabel('\Delta p [kPa]');
  subplot(2, 2, 2);
     plot(x4, p4 * 1e-3);
     legend("Man. 1", "Man. 2", "Man. 3", "Man. 4",
            "location", "south");
     title(tittxt{4});
     grid;
     ylim([-0.2, 0.4]);
  subplot(2, 2, 3);
     plot(x2, p2 * 1e-3);
     title(tittxt{2});
     grid;
     ylim([0.1, 0.5]);
     xlabel('x [m]'); ylabel('\Delta p [kPa]');
  subplot(2, 2, 4);
     plot(x3, p3 * 1e-3);
     title(tittxt{3});
     grid;
     ylim([0.1, 0.5]);
     xlabel('x [m]');
  print(["p_wing", EXT], FORMAT);

# Plot pressure in vertical stabilizer section

  [xv, pv, ~, zv] = mfs_xydata(glider, "trim", "pressure",
                               [31, 32], 4);

  figure(2, "position", [500, 200, 500, 500],
            "paperposition", [0, 0, 10, 10]);
  plot(xv, pv * 1e-3)
  title(sprintf("z = %6.3f m", zv));
  legend("Man. 1", "Man. 2", "Man. 3", "Man. 4");
  grid;
  xlabel('x [m]'); ylabel('\Delta p [kPa]');
  print(["p_vert", EXT], FORMAT);

# Compute radius of turn and yaw string angle

  tp = mfs_getresp(glider, "trim", "params");

  R2 = vv / (g * tg);
  R3 = vv / (sn * tp.az(3) - cs * tp.ay(3));

  fprintf(fid, "\nRadius of truly banked turn       : %6.2f m\n",
          R2);
  fprintf(fid, "Radius of turn with zero side slip: %6.2f m\n",
          R3);

  ysa = (tp.beta + ds * tp.yaw) * 180 / pi;

  fprintf(fid, "\nYaw string angles:\n");
  fprintf(fid, "  %6.3f, %6.3f, %6.3f, %6.3f deg.\n", ysa);

  fclose(fid);
