# Example: Gull wing
#
# Aerodynamic analysis of rigid wing
#  - Definition of the aerodynamic model
#  - Analysis of two configurations
#  - Storage of aerodynamic component including pressure 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");

# Data (N, mm, s)
# ---------------

  ya = 4000;  % y-position of aileron in mm
  yt = 6000;  % y-position of wing tip in mm
  zt = 1000;  % z-position of wing tip in mm

  cr = 1500;  % Chord length at wing root
  ct = 1000;  % Chord length at wing tip
  fr =  0.2;  % Flap chord ratio of aileron

  nxi =  50;  % No. of panels in x-direction, inner wing
  nxo =  40;  % No. of panels in x-direction, outer wing
  nxa =  10;  % No. of panels in x-direction, aileron

  nyi =  25;  % No. of panels in y-direction, inner wing
  nyo =  30;  % No. of panels in y-direction, outer wing

  alpha =        2;   % Angle of attack in degrees
  eta   =   {0, 2};   % Aileron angles in degrees
  v     =     40E3;   % Flight veloctiy in mm/s
  rho   = 1.21E-12;   % Mass density of air in t/mm^3

  frc   = 1 - fr;

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

# Model type, subtype and symmetry

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

# Airfoil data

  ppi = mfs_airfoil("NACA", 5, 30);
  ppo = mfs_airfoil("NACA", 5, 30, 0, frc);
  ppa = mfs_airfoil("NACA", 5, 30, frc, 1);

# Points on leading edge of wing

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

# Points on leading edge of aileron

  points(4 : 5) = struct("id", {12, 13},
                         "coor", {[frc * cr, ya,  0], ...
                                  [frc * ct, yt, zt]});

# Lifting surfaces of inner wing, outer wing, aileron

  lsf = struct("id",     {1, 2, 3},
               "points", {[1, 2], [2, 3], [12, 13]},
               "chord",  {cr, frc * [cr, ct], fr * [cr, ct]},
               "camber", {ppi, ppo, ppa},
               "nx",     {nxi, nxo, nxa},
               "ny",     {nyi, nyo, nyo},
               "typex",  "linear",
               "typey",  {"linear", "cos>", "cos>"});

# Aileron is control surface

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

# Configurations

  qdyn = 0.5 * rho * v^2;

  for l = 1 : 2
      cfgnam{l} = sprintf("Conf. %d: alpha = %5.2f, eta = %5.2f",
                           l, alpha, eta{l});
  endfor

  config = struct("name", cfgnam, "qdyn", qdyn,
                  "alpha", alpha, "aileron", eta);

# Add definitions to model

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

# Analysis
# --------

# Create and export component

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

# Aerodynamic analysis

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

  save aero.bin winga rho

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

# Pressure in selected wing sections

  ycols = [1, floor(0.8 * nyi), 3, floor(0.8 * nyo)];

  [xi(:, 1), pri(:, :, 1), y(1)] = ...
     mfs_xydata(winga, "statresp", "pressure", 1, ycols(1));
  [xi(:, 2), pri(:, :, 2), y(2)] = ...
     mfs_xydata(winga, "statresp", "pressure", 1, ycols(2));

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

  pri = 1000 * pri;   % Pressure in kPa
  pro = 1000 * pro;   

  xi  = xi / 1000;    % Lengths in m
  xo  = xo / 1000;
  y   =  y / 1000;

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

  crm = cr / 1000;

  figure(1, "position", [100, 500, 750, 500],
            "paperposition", [0, 0, 14, 10]);
  subplot(1, 2, 1)   % Configuration 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});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
     xlabel("x [m]"); ylabel("p [kPa]");
  subplot(1, 2, 2)   % Configuration 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});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
     xlabel("x [m]");
  print(["rigid", EXT], FORMAT);

  save rigid.plt ycols crm eta pri pro

  fclose(fid);
