# Example: Box Beam
#          Analytical frequencies
#
# ------------------------------------------------------------------

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

# Data (N, mm, t)

  E   =  210000;  % Young's modulue
  ny  =     0.3;  % Poisson's ratio
  rho = 7.85e-9;  % Mass density
  L   =    1000;  % Length
  b   =     100;  % Width
  h   =      50;  % Height
  t   =       3;  % Thickness

# Cross section properties

  geom = mfs_beamsection("box", "thin", b, h, t);

# Bending modes

  k1 = 1.8751 / L;  k2 = 4.6941 / L;
  % (cf. Gross/Hauger/Schnell/Wriggers, Technische Mechanik 4, p. 230)
  B  = sqrt(E * geom.Iy / (rho * geom.A));
  w1 = k1^2 * B;      w2 = k2^2 * B;
  f1 = w1 / (2 * pi); f2 = w2 / (2 * pi);

  fprintf(fid, "Bending Modes\n");
  fprintf(fid, "  f1 = %6.2f Hz, f2 = %6.2f Hz\n", f1, f2);

# Longitudinal mode

  fl = 0.25 * sqrt(E / rho) / L;

  fprintf(fid, "\nLongitudinal Modes\n");
  fprintf(fid, "  f1 = %6.2f Hz\n", fl);

# Torsional mode

  G  = E / (2 * (1 + ny));
  Ip = geom.Iy + geom.Iz;

  ft1 = 0.25 * sqrt(G * geom.IT / (rho * Ip)) / L;
  ft2 = 3 * ft1;

  fprintf(fid, "\nTorsional Modes\n");
  fprintf(fid, "  f1 = %6.2f Hz, f2 = %6.2f Hz\n", ft1, ft2);

  fclose(fid);
