# Example: Swept Wing
#
# Aeroelastic analysis of flexible wing
#
# Files needed: aero.bin, aeroelastic.bin
#
# ---------------------------------------------------------------------

  set(0, "defaultaxesfontsize", 10);

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

  params

# Static divergence
# -----------------

  load aeroelastic.bin

  [wings, winga, nfound] = mfs_diverg(wing, 5);

  if (nfound)
     nout = mfs_eset2nset(wings, "Stringer_1_Upper");
     mfs_print(fid, wings, "diverg", "qdyn", {"disp", nout});
     mfs_print(fid, winga, "diverg", "qdyn", "disp");
     qd = mfs_getresp(winga, "diverg", "qdyn", 1);
  else
     fclose(fid);
     return;
  end

# Lift coefficient and aileron efficiency
# ---------------------------------------

  load aero.bin

# Define new configurations

  qr   = [0.05 : 0.05 : 0.5];
  nq   = length(qr);
  qdyn = qr * qd;

  nc = 1;

  for k = 1 : 2
      for l = 1 : nq
          name = sprintf("Conf. %d: qdyn = %10.4e, eta = %5.2f",
                          nc, qdyn(l), eta(k));
          config(nc++) = struct("name",    name,
                                "qdyn",    qdyn(l),
                                "alpha",   alpha,
                                "aileron", eta(k));
      end
  end

  aero.config = config;

# Create new aerodynamic component

  winga = mfs_new(fid, aero);

# Create new aeroelastic component

  model.aero = winga;

  wing = mfs_new(fid, model);

# Compute splines and results

  wing = mfs_splines(wing);

  [wings, winga] = mfs_statresp(wing);
  winga = mfs_results(winga, "statresp", "panel");

# Compute lift coefficient

  A      = mfs_getresp(winga, "mesh", "area");
  [F, M] = mfs_getresp(winga, "statresp", "aeload");

  fprintf(fid, "     qdyn         Fx          Fy          Fz      ");
  fprintf(fid, "    Mx          My          Mz\n");
  for k = 1 : nq
      fprintf(fid, "  %10.4e  %10.4e  %10.4e  %10.4e", 
              qdyn(k), F(:, k));
      fprintf(fid, "  %10.4e  %10.4e  %10.4e\n", M(:, k));
  end

  cLF = F(3, 1 : nq) ./ (A * qdyn);
  cLr = cLF / cLR;

# Compute aileron efficiency

  dMxF = M(1, nq + 1 : end) - M(1, 1 : nq);
  dMxF = dMxF ./ qdyn;
  etaR = dMxF / dMxR;

# Plot results

  figure(1, "position", [100, 500, 1000, 500],
            "paperposition", [0, 0, 14, 10]);

  subplot(1, 2, 1);
     plot(qr, cLr, "color", "red");
     grid;
     xlabel('q_\infty / q_D');
     ylabel("c_L/c_{LR}");

  subplot(1, 2, 2);
     plot(qr, etaR, "color", "red");
     grid;
     xlabel('q_\infty / q_D');
     ylabel('\eta_R');

  print("flexible.svg", "-dsvg");

  fclose(fid);
