function model = glider_aero();

# usage: model = glider_aero();
#
# Input  none
#
# Output model   Structure containing the model definition
#
# This function defined the aerodynamic model of the glider.
#
# Units: kg, m, s
#
# ---------------------------------------------------------------

# Data
# ----

# Geometry

  b     =   15;   % Wing span
  cr    = 0.75;   % Chord length of wing at wing root
  ct    = 0.45;   % Chord length of wing at wing tip
  delta =    2;   % Dihedral angle in degrees
  lr    =  0.4;   % Aileron length ratio
  fr    =  0.2;   % Aileron flap ratio
  at    =   -2;   % Angle of incidence at wing tip (in degrees)
  l     =  6.8;   % Length of fuselage
  be    =  2.0;   % Span of elevator
  h     =  1.3;   % Height of fin
  xle   = 1.55;   % x-position of leading edge

# Discretization

  nsa =   5;      % No. of spline segments to fit airfoil
  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

# 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 type and subtype

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

# Some coordinates and lengths

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

  dihedral = tand(delta);
  zt = yt * dihedral;     % z-coordinate of wing tip
  za = ya * dihedral;     % z-coordinate at start 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: Inner wing, outer wing, aileron

  points(1 : 5) = 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]});

  lsf(1 : 3) = 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: Inner wing, outer wing, aileron

  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]});

  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

  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]});

# Vertical stabilizer: 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,
                      "typex",  "linear",
                      "typey",  "cos",
                      "alpha",  0);

# Horizontal stabilizer: Points

  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]});

# Horizontal stabilizer: 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);

# Horizontal stabilizer: 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
#    Positive aileron : stick to the right
#    Positive rudder  : right pedal
#    Positive elevator: pull

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

# Mass properties

  model.mass = struct("m",  309,
                      "cm", [1.796, 0., 0.103],
                      "JS", [  2467,     0, -72.43;
                                  0, 662.8,      0;
                              -72.43,    0,   3083]);

endfunction
