# Example: Swept Wing
#
# Aerodynamic analysis of rigid wing
#  - Definition of the aerodynamic model
#  - Computation of pressure coefficient, lift coefficient and
#    increase of rolling moment due to aileron deflection
#  - Storage of aerodynamic model and component in binary file 
#    aero.bin
#  - Storage of xy-plot data in binary file rigid.plt
#
# ---------------------------------------------------------------------

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

  set(0, "defaultaxesfontsize", 10);

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

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

# Parameters

  params
  frc = 1 - fr;

# Airfoil data

  cmbi = mfs_airfoil("NACA", 5, 30);
  cmbo = mfs_airfoil("NACA", 5, 30, 0, frc);
  cmba = mfs_airfoil("NACA", 5, 30, frc, 1);

# Model type, subtype and symmetry

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

# Points on leading edge of wing

  sa = ya / yt;
  ca = (1 - sa) * cr + sa * ct;
  sw = tand(sweep);

  xr = -0.25 * cr;
  xa = ya * sw - 0.25 * ca;
  xt = yt * sw - 0.25 * ct;

  points(1 : 3) = struct("id",   {1, 2, 3},
                         "coor", {[xr,  0, 0], [xa, ya, 0], ...
                                  [xt, yt, 0]});

# Points on leading edge of aileron

  xa += frc * ca; xt += frc * ct;

  points(4 : 5) = struct("id",   {4, 5},
                         "coor", {[xa, ya, 0], [xt, yt, 0]});

# Lifting surfaces of inner wing, outer wing and aileron

  lsf = struct("id",     {1, 2, 3},
               "points", {[1, 2], [2, 3], [4, 5]},
               "chord",  {[cr, ca], frc * [ca, ct], ...
                          fr * [ca, ct]},
               "camber", {cmbi, cmbo, cmba},
               "nx",     {nxi, nxo, nxa},
               "ny",     {nyi, nyo, nyo},
               "typex",  "linear",
               "typey",  {"linear", "cos>", "cos>"},
               "alpha",  {0, [0, atip], [0, atip]});
 
# Aileron is a control surface

  controls = struct("name", "aileron", "ls", 3);

# Configurations

  for n = 1 : 2
      cfgname = sprintf("Conf. %d: alpha = %5.2f, eta = %5.2f",
                        n, alpha, eta(n));
      config(n) = struct("name", cfgname,
                         "alpha", alpha, "aileron", eta(n));
  endfor

# Add definitions to structure aero

  aero.points   = points;
  aero.ls       = lsf;
  aero.controls = controls;
  aero.config   = config;

# Analysis
# --------

  winga = mfs_new(fid, aero);
  mfs_export("aero.msh", "msh", winga, "mesh", "camber");

  winga = mfs_statresp(winga);
  winga = mfs_results(winga, "statresp", "panel");

  mfs_export("rigid.pos", "msh", winga, "statresp", "pressure");

# Lift coefficient

  [F, M] = mfs_getresp(winga, "statresp", "aeload");
  A      = mfs_getresp(winga, "mesh", "area");
  cLR    = F(3, 1) / A;
  dMxR   = M(1, 2) - M(1, 1);

  fprintf(fid, "Area = %11.4e\n", A);
  fprintf(fid, "cLR  = %11.4e, dMxR = %11.4e\n", cLR, dMxR);

  save -binary aero.bin aero winga cLR dMxR

# Pressure coefficient in selected wing sections

  [xi(:, 1), pri(:, :, 1), y(1)] = ...
     mfs_xydata(winga, "statresp", "pressure", 1, ycoli(1));
  [xi(:, 2), pri(:, :, 2), y(2)] = ...
     mfs_xydata(winga, "statresp", "pressure", 1, ycoli(2));
  [xo(:, 1), pro(:, :, 1), y(3)] = ...
     mfs_xydata(winga, "statresp", "pressure", [2, 3], ycolo(1));
  [xo(:, 2), pro(:, :, 2), y(4)] = ...
     mfs_xydata(winga, "statresp", "pressure", [2, 3], ycolo(2));

  xi = xi / cr; xo = xo / cr;
  y  =  y * 1e-3;  % convert to m

  save -binary rigid.plt pri pro

  for n = 1 : length(y)
      leg{n} = sprintf("y = %4.2f m", y(n));
  endfor
  for n = 1 : length(eta)
      header{n} = sprintf("\\eta = %2.0f\\deg", eta(n));
  endfor

  figure(1, "position", [100, 500, 1000, 500],
            "paperposition", [0, 0, 16, 10]);
  subplot(1, 2, 1)
     plot(xi(:, 1), squeeze(pri(:, 1, 1)),
          xi(:, 2), squeeze(pri(:, 1, 2)),
          xo(:, 1), squeeze(pro(:, 1, 1)),
          xo(:, 2), squeeze(pro(:, 1, 2)));
     legend(leg, "location", "northeast");
     title(header{1});
     ylim([-0.1, 1]);
     grid;
     xlabel('x / c'); ylabel('\Delta c_P');
  subplot(1, 2, 2)
     plot(xi(:, 1), squeeze(pri(:, 2, 1)),
          xi(:, 2), squeeze(pri(:, 2, 2)),
          xo(:, 1), squeeze(pro(:, 2, 1)),
          xo(:, 2), squeeze(pro(:, 2, 2)));
     title(header{2});
     ylim([-0.1, 1]);
     grid;
     xlabel('x / c');
  print(["rigid", EXT], FORMAT);

  fclose(fid);
