# Example: Gull wing
#
# --------------------------------------------------------------------

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

# Data (mm, t, 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 =   {0, 0, 2, 2};   % Angles of attack in degrees
  eta   =   {0, 5, 0, 5};   % Aileron angles in degrees
  v     =     40E3;         % Flight velocity in mm/s
  rho   = 1.21E-12;         % Mass density of air in t/mm^3

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

# Model type, subtype and symmetry

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

# Airfoil data

  frc = 1 - fr;

  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).id = 1; points(1).coor = [0,  0,  0];
  points(2).id = 2; points(2).coor = [0, ya,  0];
  points(3).id = 3; points(3).coor = [0, yt, zt];

# Points on leading edge of aileron

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

# Lifting surface of inner wing

  lsf(1).id = 1;     lsf(1).points = [1, 2]; 
  lsf(1).chord = cr; lsf(1).camber = ppi;
  lsf(1).nx = nxi;   lsf(1).ny = nyi;
  lsf(1).typex = "linear";

# Lifting surface of outer wing

  lsf(2).id = 2;                 lsf(2).points = [2, 3]; 
  lsf(2).chord = frc * [cr, ct]; lsf(2).camber = ppo;
  lsf(2).nx = nxo;               lsf(2).ny = nyo;
  lsf(2).typex = "linear";       lsf(2).typey = "cos>";

# Lifting surface of aileron

  lsf(3).id = 3;                lsf(3).points = [12, 13]; 
  lsf(3).chord = fr * [cr, ct]; lsf(3).camber = ppa;
  lsf(3).nx = nxa;              lsf(3).ny = nyo;
  lsf(3).typex = "linear";      lsf(3).typey = "cos>";

# Aileron is control surface

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

# Configurations

  qdyn = 0.5 * rho * v^2;

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

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

# Add definitions to the model

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

# Analysis
# --------

# Create and export component

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

# Compute and export results

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

# Aerodynamic coefficients

  [F, M] = mfs_getresp(wing, "statresp", "aeload");
  As     = mfs_getresp(wing, "mesh", "area", 1 : 3);

  tg = zt / (yt - ya);
  cs = 1 / sqrt(1 + tg^2);
  A  = As(1) + (As(2) + As(3)) * cs;

  cL = F(3, :) / (qdyn * A);
  cM = M(2, :) / (qdyn * A * cr);

  dalpha = pi * (alpha{3} - alpha{1}) / 180; 
  deta   = pi * (eta{2} - eta{1}) / 180;

  dcLda = (cL(3) - cL(1)) / dalpha;
  dcLde = (cL(2) - cL(1)) / deta;

  fprintf(fid, "\nWing area = %10.5e, ", A);
  fprintf(fid, "dcLda = %7.4f, dcLde = %7.4f\n\n", 
               dcLda, dcLde);

  for n = 1 : length(config)
      fprintf(fid, "%s: ", config(n).name);
      fprintf(fid, "cL = %7.4f, cM = %7.4f\n", cL(n), cM(n));
  endfor

  fclose(fid);
