function model = solid(damping_type);

# Input   damping_type    Damping type: "ratios" or "Rayleigh"
# Output  model           Structure with model definition

# -------------------------------------------------------------------
# Example: Frequency response analysis of a standard class glider
#          Solid Model
# -------------------------------------------------------------------

  if (! nargin)
     damping_type = "ratios";      % Default
  endif

# ====================================================================
# Data (kg, m, s)
# ====================================================================

# Geometry

  geodata

# Point masses

  m   =  85;   % Mass of the pilot
  mt  =   5;   % Trim mass

# Material

  mat = struct("type", "iso", "E",   4E10,
               "ny",     0.3, "rho", 2000);

# Analysis

  nofmod = 5;  % Number of normal modes

# x-positions

  xwr = xle + 0.25 * cr;         % Wing root
  xtr = xle + (1 - fr) * cr;     % Trailing edge of inner wing
  xtt = xle + (1 - fr) * ct;     % Trailing edge at wing tip

# y-positions

  yt = 0.5 * b;                  % Wing tip
  ya = yt * (1 - lr);            % Begin of aileron
  y1 = yt / 15;                  % First rib

# z-positions

  dihedral = tand(delta);
  zt = yt * dihedral;            % Wing tip
  za = ya * dihedral;            % Begin of aileron
  z1 = y1 * dihedral;            % First rib

# ====================================================================
# Model Definition
# ====================================================================

  model = struct("type", "solid", "subtype", "3d");

# Fuselage
# --------

# Cross section properties: Radius and thickness

  fuselage.r = [0.2000, 0.3000, 0.4000, 0.4000, 0.3667, ...
                0.2667, 0.1367, 0.1140, 0.1100, 0.1080, ...
                0.1060, 0.1040, 0.1020, 0.1000];
  fuselage.t = [ 0.005,  0.006,  0.006,  0.005,  0.005, ...
                 0.004,  0.003,  0.003,  0.003,  0.003, ...
                 0.003,  0.003,  0.003,  0.003];

# Cabin nodes

  nodes(1).id = 1; nodes(1).coor = [  0, 0, 0]; % Tip
  nodes(2).id = 5; nodes(2).coor = [xwr, 0, 0]; % Wing root

  nodes = mfs_linenodes(nodes, 1, 5, 2 : 4);

# Tail nodes

  nodes(6).id = 15; nodes(6).coor = [l, 0, 0];  % End Node
  nodes = mfs_linenodes(nodes, 5, 15, 6 : 14);

# Elements

  for k = 1 : 14
      geom = mfs_beamsection("ring", "thin",
                             fuselage.r(k), fuselage.t(k));
      geom.v = [0, 0, 1];
      elem(k) = struct("id", k, "type", "b2", 
                       "nodes", [k, k + 1],
                       "geom", geom, "mat", mat);
  endfor

# Wing beam sections
# ------------------

  geom_spar   = mfs_beamsection("box", "thin", 0.4, 0.2, 0.003);
  geom_spar.v = [0, 0, 1];

  geom_te   = mfs_beamsection("bar", 0.001, 0.05);
  geom_te.v = [0, 0, 1];

  geom_le   = mfs_beamsection("bar", 0.001, 0.05);
  geom_le.v = [0, 0, 1];

  geom_ribs   = mfs_beamsection("I", 0.1, 0.1, 0.002, 0.002);
  geom_ribs.v = [0, 0, 1];

# Left Wing
# ---------

# Spar

  nnode = length(nodes) + 1;
  nodes(nnode).id = 109; nodes(nnode++).coor = [xwr, -ya, za];
  nodes(nnode).id = 115; nodes(nnode).coor = [xwr, -yt, zt];
  [nodes, ells1] = mfs_line(nodes, 5, 109, 101 : 108, 101 : 109,
                           "b2", geom_spar, mat);
  [nodes, ells2] = mfs_line(nodes, 109, 115, 110 : 114, 110 : 115,
                           "b2", geom_spar, mat);
  eset.left_spar = 101 : 115;

# Trailing Edge

  nnode = length(nodes) + 1;
  nodes(nnode).id = 121; nodes(nnode++).coor = [xtr, -y1, z1];
  nel = length(elem) + 1;
  elem(nel).id = 121; elem(nel).type = "b2";
  elem(nel).nodes = [6, 121];
  elem(nel).geom = geom_te; elem(nel).mat = mat;

  nodes(nnode).id = 129; nodes(nnode++).coor = [xtr, -ya , za];
  nodes(nnode).id = 135; nodes(nnode).coor = [xtt, -yt, zt];
  [nodes, ellt1] = mfs_line(nodes, 121, 129, 122 : 128, 122 : 129,
                            "b2", geom_te, mat);
  [nodes, ellt2] = mfs_line(nodes, 129, 135, 130 : 134, 130 : 135,
                            "b2", geom_te, mat);

# Leading edge

  nnode = length(nodes) + 1;
  nodes(nnode).id = 141; nodes(nnode++).coor = [xle, -y1, z1];
  nodes(nnode).id = 149; nodes(nnode++).coor = [xle, -ya, za];
  nodes(nnode).id = 155; nodes(nnode).coor = [xle, -yt, zt];
  [nodes, elll1] = mfs_line(nodes, 141, 149, 142 : 148, 142 : 149,
                           "b2", geom_le, mat);
  [nodes, elll2] = mfs_line(nodes, 149, 155, 150 : 154, 150 : 155,
                            "b2", geom_le, mat);

  elem = [elem, ells1, ells2, ellt1, ellt2, elll1, elll2];

# Ribs

  nel1 = length(elem) + 1;
  nel2 = nel1 + 14;

  n1 = 141; n2 = 101; n3 = 121; ne = nel1; ide = 161;
  for k = nel1 : nel2
      elem(ne++) = struct("id", ide++, "type", "b2",
                          "nodes", [n1++, n2],
                          "geom", geom_ribs, "mat", mat);
      elem(ne++) = struct("id", ide++, "type", "b2",
                          "nodes", [n2++, n3++],
                          "geom", geom_ribs, "mat", mat);
  endfor

  nset.left_wing = [5 : 6, 101 : 115, 121 : 135, 141 : 155];

# Right wing
# ----------

# Spar

  nnode = length(nodes) + 1;
  nodes(nnode).id = 209; nodes(nnode++).coor = [xwr, ya, za];
  nodes(nnode).id = 215; nodes(nnode).coor = [xwr, yt, zt];
  [nodes, elrs1] = mfs_line(nodes, 5, 209, 201 : 208, 201 : 209,
                            "b2", geom_spar, mat);
  [nodes, elrs2] = mfs_line(nodes, 209, 215, 210 : 214, 210 : 215,
                           "b2", geom_spar, mat);
  eset.right_spar = 201 : 215;

# Trailing Edge

  nnode = length(nodes) + 1;
  nodes(nnode).id = 221; nodes(nnode++).coor = [xtr,  y1, z1];
  nel = length(elem) + 1;
  elem(nel).id = 221; elem(nel).type = "b2";
  elem(nel).nodes = [6, 221];
  elem(nel).geom = geom_te; elem(nel).mat = mat;

  nodes(nnode).id = 229; nodes(nnode++).coor = [xtr, ya, za];
  nodes(nnode).id = 235; nodes(nnode).coor = [xtt, yt, zt];
  [nodes, elrt1] = mfs_line(nodes, 221, 229, 222 : 228, 222 : 229,
                            "b2", geom_te, mat);
  [nodes, elrt2] = mfs_line(nodes, 229, 235, 230 : 234, 230 : 235,
                            "b2", geom_te, mat);

# Leading edge

  nnode = length(nodes) + 1;
  nodes(nnode).id = 241; nodes(nnode++).coor = [xle,  y1, z1];
  nodes(nnode).id = 249; nodes(nnode++).coor = [xle,  ya, za];
  nodes(nnode).id = 255; nodes(nnode).coor = [xle,  yt, zt];
  [nodes, elrl1] = mfs_line(nodes, 241, 249, 242 : 248, 242 : 249,
                            "b2", geom_le, mat);
  [nodes, elrl2] = mfs_line(nodes, 249, 255, 250 : 254, 250 : 255,
                            "b2", geom_le, mat);

  elem = [elem, elrs1, elrs2, elrt1, elrt2, elrl1, elrl2];

# Ribs

  nel1 = length(elem) + 1;
  nel2 = nel1 + 14;

  n1 = 241; n2 = 201; n3 = 221; ne = nel1; ide = 261;
  for k = nel1 : nel2
      elem(ne++) = struct("id", ide++, "type", "b2",
                          "nodes", [n1++, n2],
                          "geom", geom_ribs, "mat", mat);
      elem(ne++) = struct("id", ide++, "type", "b2",
                          "nodes", [n2++, n3++],
                          "geom", geom_ribs, "mat", mat);
  endfor

  nset.right_wing = [5 : 6, 201 : 215, 221 : 235, 241 : 255];

# Vertical stabilizer
# -------------------

  geom  = mfs_beamsection("box", "thin", 0.10, 0.4, 0.002); 
  geom.v = [1, 0, 0];

  nnode = length(nodes) + 1;
  nodes(nnode).id = 304; nodes(nnode).coor = [l, 0, h];
  [nodes, elvsr] = mfs_line(nodes, 15, 304, 303, 303 : 304,
                           "b2", geom, mat);
   
  geom  = mfs_beamsection("bar", 0.05, 0.002); 
  geom.v = [1, 0, 0];
 
  nnode = length(nodes) + 1;
  nodes(nnode).id = 302; nodes(nnode).coor = [l - 0.4, 0, h];
  [nodes, elvsf] = mfs_line(nodes, 14, 302, 301, 301 : 302,
                           "b2", geom, mat);
   
  elem = [elem, elvsr, elvsf];

  geom  = mfs_beamsection("bar", 0.05, 0.002); 
  geom.v = [0, 0, 1];

  nel = length(elem) + 1;
  id1 = 301; id2 = 303; ide = 305;
  for k = 1 : 2
      elem(nel++) = struct("id", ide++, "type", "b2",
                           "nodes", [id1++, id2++],
                           "geom", geom, "mat", mat);
  endfor

  nset.vertical = [14, 15, 301 : 304];

# Horizontal stabilizer (root is at node 302)
# -------------------------------------------

  d = 0.5 * be;

  geomr   = mfs_beamsection("box", "thin", 0.3, 0.08, 0.002);
  geomr.v = [0, 0, 1];

  geomf   = mfs_beamsection("box", "thin", 0.1, 0.05, 0.002);
  geomf.v = [0, 0, 1];

  nnode = length(nodes) + 1;
  nodes(nnode).id = 402; nodes(nnode++).coor = [l - 0.3, -d, h];
  nodes(nnode).id = 412; nodes(nnode++).coor = [l - 0.3,  d, h];
  nodes(nnode).id = 404; nodes(nnode++).coor = [l, -d, h];
  nodes(nnode).id = 414; nodes(nnode++).coor = [l,  d, h];

  [nodes, elhsfl] = mfs_line(nodes, 302, 402, 401, 401 : 402,
                            "b2", geomf, mat);
  [nodes, elhsrl] = mfs_line(nodes, 304, 404, 403, 403 : 404,
                            "b2", geomr, mat);

  [nodes, elhsfr] = mfs_line(nodes, 302, 412, 411, 411 : 412,
                            "b2", geomf, mat);
  [nodes, elhsrr] = mfs_line(nodes, 304, 414, 413, 413 : 414,
                            "b2", geomr, mat);

  elem = [elem, elhsfl, elhsrl, elhsfr, elhsrr];

  geom   = mfs_beamsection("bar", 0.002, 0.04);
  geom.v = [0, 0, 1];

  nel = length(elem) + 1;
  elem(nel++) = struct("id", 421, "type", "b2",
                       "nodes", [402, 404],
                       "geom", geom, "mat", mat);
  elem(nel++) = struct("id", 422, "type", "b2",
                       "nodes", [412, 414],
                       "geom", geom, "mat", mat);

  nset.horizontal = [302, 304, 401 : 404, 411 : 414];

# Pilot
# -----

  mass.m = m;
  elem(nel++) = struct("id",  500, "type", "m1",
                       "nodes", 3, "geom", mass, "mat", []);

# Trim mass
# ---------

  mass.m = mt;
  elem(nel) = struct("id",  510, "type", "m1",
                     "nodes", 1, "geom", mass, "mat", []);

  model.nodes = nodes; model.elements = elem;
  model.nset  = nset;  model.eset     = eset;   

# Damping
# -------

# Rayleigh damping is adjusted sucht that it gives D = 0.02
# at 1 Hz and at 10 Hz

  switch damping_type

  case "ratios"
     model.damping = struct("type", "ratios", "data", 0.02);

  case "Rayleigh"
     w1 = 2 * pi; w2 = 20 * pi;
     W  = 0.5 * [w1, 1 / w1; w2, 1 / w2];
     model.damping = struct("type", "Rayleigh",
                            "data", W \ [0.02; 0.02]);
     
  otherwise
     error("Unknown damping type %s\n", damping_type);

  endswitch

endfunction
