# Example: Wing with winglet and controls
#
# --------------------------------------------------------------------

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

  set(0, "defaultaxesfontsize", 10);

  file = mfilename();

# Wing data

  b     = 10000;      % Half span (mm)
  s     =   0.7;      % Aileron starts at 70 % of half span
  xw    =   200;      % x-coordinate of winglet tip (mm)
  h     =  1000;      % Winglet height (mm)

  c1    =  1000;      % Chord length at wing root
  c2    =   800;      % Chord length at wing tip
  cw    =   600;      % Chord length at winglet tip
  fr    =   0.2;      % Flap chord ratio

  alpha_tip = -2;     % Twist at wing tip (degrees)

  alpha   = [0, 2];     % Angle of attack (degrees)
  flap    = [-2, 0, 2]; % Flap angle (degrees)
  aileron = [-4, 0, 4]; % Aileron angle (degrees)

  delta  = 1;         % Dihedral angle (degrees)

# Discretization data

  nxw = 40;   % Number of wing panels in x-direction
  nxf = 10;   % Number of flap panels in x-direction
  nyi = 20;   % Number of inner wing panels in y-direction
  nyo = 15;   % Number of outer wing panels in y-direction
  nyw = 15;   % Number of winglet panels in y-direction

  ns  = 10;   % Number of spline segments

# Additional data

  mass   =  650;      % Mass (kg)
  g      = 9.81;      % Gravity acceleration (m/s^2)
  rho    = 1.21;      % Mass density of air (kg/m^3)

  d      = tand(delta);
  frc    = 1 - fr;

# Airfoil data: Horstmann-Quast HQ17

  data   = dlmread("hq17.csv", ",", "A108:B155");

  airfoil_wing = mfs_airfoil("fit", "camber", data, ns, 0, frc);
  airfoil_flap = mfs_airfoil("fit", "camber", data, ns, frc, 1);

# Model definition
# ----------------

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

# Leading edge points

  ya = s * b; za = d * ya;
  zt = d * b;

  points(1 : 4) = struct("id",   {1, 2, 3, 4},
                         "coor", {[0, 0, 0], [0, ya, za], ...
                                  [0, b, zt], [xw, b, zt + h]});

# Hinge points

  xr = frc * c1; xt = frc * c2;

  points(5 : 7) = struct("id",   {11, 12, 13},
                         "coor", {[xr, 0, 0], [xr, ya, za], ...
                                  [xt, b, zt]});

  model.points = points;

# Lifting surfaces of wing and winglet

  lsf(1 : 3) = struct("id",     {10, 20, 30},
                      "points", {[1, 2], [2, 3], [3, 4]},
                      "chord",  {frc * c1, frc * [c1, c2], ...
                                 [c2, cw]},
                      "camber", {airfoil_wing, airfoil_wing, []},
                      "nx",     nxw,
                      "ny",     {nyi, nyo, nyw},
                      "typex",  "linear",
                      "typey",  {"linear", "cos>", "cos"},
                      "alpha",  {0, [0, alpha_tip], 0});

# Lifting surfaces of flap and aileron

  lsf(4 : 5) = struct("id",     {11, 21},
                      "points", {[11, 12], [12, 13]},
                      "chord",  {fr * c1, fr * [c1, c2]},
                      "camber", airfoil_flap,
                      "nx",     nxf,
                      "ny",     {nyi, nyo},
                      "typex",  "linear",
                      "typey",  {"linear", "cos>"},
                      "alpha",  {0, [0, alpha_tip]});
                    
  model.ls = lsf;

# Control surfaces

  model.controls = struct("name",    {"flap", "aileron"},
                          "ls",      {[11, 21], [11, 21]},
                          "factors", {[1, 1], [0.1, 1]});

# Configurations

  nc = 1;
  for l = 1 : length(alpha)
    for m = 1 : length(flap)
      for n = 1 : length(aileron)
        name = sprintf("Conf. %2d: alpha = %4.1f, ",
                       nc, alpha(l));
        name = [name, sprintf("flap = %4.1f, ", flap(m))];
        name = [name, sprintf("ail. = %4.1f", aileron(n))];
        config(nc++) = struct("name",    name,
                              "alpha",   alpha(l),
                              "flap",    flap(m),
                              "aileron", aileron(n));
      endfor
    endfor
  endfor

  model.config   = config;

# Analysis
# --------

  fid = fopen([file, ".res"], "wt");

# Create and export component

  wing = mfs_new(fid, model);
  mfs_export([file, ".msh"], "msh", wing, "mesh", "camber");

# Wing area

  A     = mfs_getresp(wing, "mesh", "area", [10, 11, 20, 21]);
  Awing = sum(A);

# Compute and export results

  wing = mfs_statresp(wing);
  wing = mfs_results(wing, "statresp", "panel");
  mfs_export([file, ".pos"], "msh", wing, "statresp", "pressure");

# Aerodynamic coefficients

  [F, M] = mfs_getresp(wing, "statresp", "aeload");

  cL = F(3, :) / Awing;
  cM = M(2, :) / (Awing * c1);

# Flight velocity

  v = (1e6 * mass * g / (rho * Awing)) ./ cL;
  v = 3.6 * sqrt(v);

# Output results

  fprintf(fid, "\nWing Area = %10.5e mm^2\n", Awing);

  fprintf(fid, "\n                                        ");
  fprintf(fid, "            cL       cM     v [km/h] ");
  fprintf(fid, "Mx [Nmm/MPa]\n");
  for nc = 1 : length(config)
      fprintf(fid, "%s: ", config(nc).name);
      fprintf(fid, "%7.4f  %7.4f  %7.2f  %10.5e\n",
              cL(nc), cM(nc), v(nc), M(1, nc));
  end

# xy-plots of pressure coefficient

  cfg = [5, 8, 9];

  for k = 1 : 3
      lgtxt{k} = sprintf("Config. %2.0d", cfg(k));
  endfor

  [x1, cp1, y1] = ...
      mfs_xydata(wing, "statresp", "pressure", [10, 11], 1, cfg);
  [x2, cp2, y2] = ...
      mfs_xydata(wing, "statresp", "pressure", [10, 11], 
                 floor(0.7 * nyi), cfg);
  [x3, cp3, y3] =...
      mfs_xydata(wing, "statresp", "pressure", [20, 21], 1, cfg);
  [x4, cp4, y4] = ...
      mfs_xydata(wing, "statresp", "pressure", [20, 21], 
                 floor(0.7 * nyo), cfg);

  figure(1, "position", [100, 500, 1000, 1000],
            "paperposition", [0, 0, 15, 15]);
  subplot(2, 2, 1)
     plot(x1, cp1);
     ylim([-0.2, 1]);
     title(sprintf("y = %4.0f mm", y1));
     legend(lgtxt, "location", "south");
     grid;
     ylabel('\Delta c_p');
  subplot(2, 2, 2)
     plot(x2, cp2);
     ylim([-0.2, 1]);
     title(sprintf("y = %4.0f mm ", y2));
     grid;
  subplot(2, 2, 3)
     plot(x3, cp3);
     ylim([-0.2, 1.2]);
     title(sprintf("y = %4.0f mm ", y3));
     grid;
     ylabel('\Delta c_p');
     xlabel('x [mm]');
  subplot(2, 2, 4)
     plot(x4, cp4);
     ylim([-0.2, 1.2]);
     title(sprintf("y = %4.0f mm ", y4));
     grid;
     xlabel('x [mm]');
  print(["wing_cp", EXT], FORMAT);

  fclose(fid);
