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

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

# Geometry

  geodata

# Discretization

  nsa =   5;      % No. of airfoil spline segments
  nxi =  50;      % No. of panels in x-direction of inner wing
  nxo =  40;      % No. of panels in x-direction of outer wing
  nxa =  10;      % No. of panels in x-direction of aileron
  nyi =  20;      % No. of panels in y-direction of inner wing
  nyo =  24;      % No. of panels in y-direction of outer wing

# Environment and Manoeuvres

  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

# Airfoil: FX 66-S-196 V1

  zs   = dlmread("fx66s196v1.csv", ",", "A100:B143");
  cmbi = mfs_airfoil("fit", "camber", zs, nsa);
  cmbo = mfs_airfoil("fit", "camber", zs, nsa, 0, 1 - fr);
  cmba = mfs_airfoil("fit", "camber", zs, nsa, 1 - fr, 1);

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

  model = struct("type", "aero", "subtype", "vlm");

  yt  = 0.5 * b;           % y-coordinate of wing tip
  ya  = (1 - lr) * yt;     % y-coordinate of aileron

  dihedral = tand(delta);
  zt = yt * dihedral;      % Wing tip
  za = ya * dihedral;      % Begin of aileron

  ca1 = fr * cr;           % inner aileron chord length
  ca2 = fr * ct;           % outer aileron chord length
  co1 = cr - ca1;          % inner outer wing chord length
  co2 = ct - ca2;          % outer outer wing chord length

  aa  = (1 - lr) * at;     % Angle of incidence at start of aileron

# Right wing
# ----------

  points = struct("id",   {1, 2, 3, 4, 5},
                  "coor", {[xle, 0, 0], ...
                           [xle, ya, za], ...
                           [xle, yt, zt], ...
                           [xle + co1, ya, za], ...
                           [xle + co2, yt, zt]});

# Inner wing, outer wing, aileron

  lsf = struct("id",     {10, 11, 12},                   
               "points", {[1, 2], [2, 3], [4, 5]},
               "chord",  {cr, [co1, co2], [ca1, ca2]},
               "camber", {cmbi, cmbo, cmba},
               "nx",     {nxi, nxo, nxa},
               "ny",     {nyi, nyo, nyo},
               "typex",  "linear",
               "typey",  {"linear", "cos>", "cos>"},
               "alpha",  {[0, aa], [aa, at], [aa, at]});

# -------------------------------------------------------------------

# Left wing
# ---------

  points(6 : 10) = struct("id",   {6, 7, 8, 9, 10},
                          "coor", {[xle, 0, 0], ...
                                   [xle, -ya, za], ...
                                   [xle, -yt, zt], ...
                                   [xle + co1, -ya, za], ...
                                   [xle + co2, -yt, zt]});

# Inner wing, outer wing, aileron

  lsf(4 : 6) = struct("id",     {20, 21, 22},
                      "points", {[7, 6], [8, 7], [10, 9]},
                      "chord",  {cr, [co2, co1], [ca2, ca1]},
                      "camber", {cmbi, cmbo, cmba},
                      "nx",     {nxi, nxo, nxa},
                      "ny",     {nyi, nyo, nyo},
                      "typex",  "linear",
                      "typey",  {"linear", "<cos", "<cos"},
                      "alpha",  {[aa, 0], [at, aa], [at, aa]});

# -------------------------------------------------------------------

# Vertical stabilizer
# -------------------

  points(11 : 14) = struct("id",   {31, 32, 33, 34},
                           "coor", {[l - 0.6, 0, 0], ...
                                    [l - 0.5, 0, h], ...
                                    [l, 0, 0], ...
                                    [l, 0, h]});
# Fin and rudder

  lsf(7 : 8) = struct("id",     {31, 32},
                      "points", {[31, 32], [33, 34]},
                      "chord",  {[0.6, 0.5], 0.2},
                      "camber", [],
                      "nx",     {30, 10},
                      "ny",     { 8,  8},
                      "typex",  "linear",
                      "typey",  "cos",
                      "alpha",  0);

# -------------------------------------------------------------------

# Horizontal stabilizer
# ---------------------

  points(15 : 18) = struct("id",   {41, 42, 43, 44},
                           "coor", {[l - 0.4, 0.5 * be, h], ...
                                    [l, 0.5 * be, h], ...
                                    [l - 0.4, -0.5 * be, h], ...
                                    [l, -0.5 * be, h]});

# Fin

  lsf(9 : 10) = struct("id",     {41, 42},
                       "points", {[32, 41], [43, 32]},
                       "chord",  {[0.5, 0.4], [0.4, 0.5]},
                       "camber", [],
                       "nx",     30,
                       "ny",      8,
                       "typex",  "linear",
                       "typey",  {"cos>", "<cos"},
                       "alpha",  -3);

# Elevator

  lsf(11 : 12) = struct("id",     {43, 44},
                        "points", {[34, 42], [44, 34]},
                        "chord",  0.1, 
                        "camber", [],
                        "nx",     6,
                        "ny",     8,
                        "typex",  "linear",
                        "typey",  {"cos>", "<cos"},                       
                        "alpha",  -3); 

# -------------------------------------------------------------------

  model.points = points;
  model.ls     = lsf;

# Control surfaces
# ----------------

# Aileron  (positive: stick to the right)
# Rudder   (positive: right pedal)
# Elevator (positive: pull)

  controls = struct("name",    {"aileron", "rudder", "elevator"},
                    "ls",      { [12, 22],    32,     [43, 44]},
                    "factors", { [-1, 1],      1,     [-1, -1]});

  model.controls = controls;

# Mass properties
# ---------------

  load solid_mass.bin
  model.mass = mp;

# 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

  glidera = mfs_new(fid, model);
  mfs_export("aero.msh", "msh", glidera, "mesh");

  save -binary aero.bin gamma v g glidera

# Perform the trim analysis

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

# Compute and export results

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

# Get load resultants

  [F, M] = mfs_getresp(glidera, "trim", "aeload", mp.cm);
  nconf = columns(F);

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

# Get pressure in some wing sections

  ycmi  = floor(0.5 * nyi); ycmo = floor(0.5 * nyo);
  ycoli = [ycmi, nyi + 1 - ycmi];
  ycolo = [ycmo, nyo + 1 - ycmo];

  [x1, p1r, y(1)] = mfs_xydata(glidera, "trim", "pressure",
                              [21, 22], ycolo(1));
  [x2, p2r, y(2)] = mfs_xydata(glidera, "trim", "pressure",
                              20, ycoli(1));
  [x3, p3r, y(3)] = mfs_xydata(glidera, "trim", "pressure",
                              10, ycoli(2));
  [x4, p4r, y(4)] = mfs_xydata(glidera, "trim", "pressure",
                              [11, 12], ycolo(2));

  save -binary rigid.bin ycoli ycolo p1r p2r p3r p4r

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

  figure(1, "position", [200, 200, 1000, 1000],
            "paperposition", [0, 0, 15, 15]);
  subplot(2, 2, 1);
     plot(x1, p1r * 1e-3);
     title(tittxt{1});
     grid;
     ylim([0, 0.4]);
     ylabel('\Delta p [kPa]');
  subplot(2, 2, 2);
     plot(x4, p4r * 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, p2r * 1e-3);
     title(tittxt{2});
     grid;
     ylim([0.1, 0.5]);
     xlabel('x [m]'); ylabel('\Delta p [kPa]');
  subplot(2, 2, 4);
     plot(x3, p3r * 1e-3);
     title(tittxt{3});
     grid;
     ylim([0.1, 0.5]);
     xlabel('x [m]');
  print(["pr", EXT], FORMAT);

# Compute radius of turn and yaw string angle

  tp = mfs_getresp(glidera, "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);
