# Example: Gull wing
#
# Aeroelastic analysis of flexible wing
#
# Files needed: aeroelastic.bin, rigid.plt
#
# ---------------------------------------------------------------------

  addpath("../../..");
  [EXT, FORMAT] = iniplot();

  set(0, "defaultaxesfontsize", 10);

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

  load aeroelastic.bin

# Analysis
# --------

# Aeroelastic static response

  [wings, winga] = mfs_statresp(wing);

  mfs_print(fid, wings, "statresp", {"disp", [1301, 1341]});

  winga = mfs_results(winga, "statresp", "panel");
  mfs_export("solid.dsp", "msh", wings, "statresp", "disp");
  mfs_export("aero.prs", "msh", winga, 
             "statresp", "pressure", "disp");

# Load resultants

  [FA, MA] = mfs_getresp(winga, "statresp", "aeload");
  nc = columns(FA);

  fprintf(fid, "\nLoad resultants of aerodynamic component: \n\n");
  for k = 1 : nc
      fprintf(fid, "  Configuration %2d:\n", k)
      fprintf(fid, "     F = [%10.3e, %10.3e, %10.3e] kN\n",
              FA(:, k) / 1000);
      fprintf(fid, "     M = [%10.3e, %10.3e, %10.3e] kNm\n",
              MA(:, k) * 1e-6);
  endfor

# Comparison of pressure in wing sections

  load rigid.plt

  [xi(:, 1), pfi(:, :, 1), y(1)] = ...
     mfs_xydata(winga, "statresp", "pressure", 1, ycols(1));
  [xi(:, 2), pfi(:, :, 2), y(2)] = ...
     mfs_xydata(winga, "statresp", "pressure", 1, ycols(2));

  [xo(:, 1), pfo(:, :, 1), y(3)] = ...
     mfs_xydata(winga, "statresp", "pressure", [2, 3], ycols(3));
  [xo(:, 2), pfo(:, :, 2), y(4)] = ...
     mfs_xydata(winga, "statresp", "pressure",
                [2, 3], ycols(4));

  pfi = 1000 * pfi; pfo = 1000 * pfo;
  xi  = xi / 1000;  xo  = xo / 1000;  y   =  y / 1000;

  for k = 1 : length(y)
      header{k} = sprintf("\\eta = %1d\\deg, y = %4.2f m", 
                          eta{1}, y(k));
  endfor

  figure(1, "position", [100, 500, 750, 500],
            "paperposition", [0, 0, 14, 10]);
  subplot(2, 2, 1)
     plot(xi(:, 1), squeeze(pri(:, 1, 1)),
          xi(:, 1), squeeze(pfi(:, 1, 1)));
     legend("rigid", "flexible", "location", "northeast");
     title(header{1});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
     ylabel("p [kPa]");
  subplot(2, 2, 2)
     plot(xi(:, 2), squeeze(pri(:, 1, 2)),
          xi(:, 2), squeeze(pfi(:, 1, 2)));
     title(header{2});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
  subplot(2, 2, 3)
     plot(xo(:, 1), squeeze(pro(:, 1, 1)),
          xo(:, 1), squeeze(pfo(:, 1, 1)));
     title(header{3});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
     xlabel("x [m]"); ylabel("p [kPa]");
  subplot(2, 2, 4)
     plot(xo(:, 2), squeeze(pro(:, 1, 2)),
          xo(:, 2), squeeze(pfo(:, 1, 2)));
     title(header{4});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
     xlabel("x [m]");
  print(["flexible1", EXT], FORMAT);

  for k = 1 : length(y)
      header{k} = sprintf("\\eta = %1d\\deg, y = %4.2f m", 
                          eta{2}, y(k));
  endfor

  figure(2, "position", [250, 250, 750, 500],
            "paperposition", [0, 0, 14, 10]);
  subplot(2, 2, 1)
     plot(xi(:, 1), squeeze(pri(:, 2, 1)),
          xi(:, 1), squeeze(pfi(:, 2, 1)));
     legend("rigid", "flexible", "location", "northeast");
     title(header{1});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
     ylabel("p [kPa]");
  subplot(2, 2, 2)
     plot(xi(:, 2), squeeze(pri(:, 2, 2)),
          xi(:, 2), squeeze(pfi(:, 2, 2)));
     title(header{2});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
  subplot(2, 2, 3)
     plot(xo(:, 1), squeeze(pro(:, 2, 1)),
          xo(:, 1), squeeze(pfo(:, 2, 1)));
     title(header{3});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
     xlabel("x [m]"); ylabel("p [kPa]");
  subplot(2, 2, 4)
     plot(xo(:, 2), squeeze(pro(:, 2, 2)),
          xo(:, 2), squeeze(pfo(:, 2, 2)));
     title(header{4});
     xlim([0, crm]); ylim([0, 0.9]);
     grid;
     xlabel("x [m]");
  print(["flexible2", EXT], FORMAT);

  fclose(fid);
