# Example: Swept Wing
#
# Stress analysis at a velocity of 250 km/h
#
# Files needed: aero.bin, aeroelastic.bin
#
# ---------------------------------------------------------------------

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

  set(0, "defaultaxesfontsize", 10);

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

# Data

  params

  v   =      250;    % Flight velocity in km/h
  rho = 1.21e-12;    % Mass density of air in t/mm^3

# Define new configurations

  load aero.bin

  vmms = 1000 * v / 3.6;
  qdyn = 0.5 * rho * vmms^2;

  for n = 1 : 2
      name{n} = sprintf("Conf. %d: qdyn = %10.4e, eta = %5.2f",
                        n, qdyn, eta(n));
  endfor
  aero.config = struct("name", name, "qdyn", qdyn, "alpha", alpha,
                       "aileron", num2cell(eta));

# Create new aerodynamic component

  winga = mfs_new(fid, aero);

# Create new aeroelastic component

  load aeroelastic.bin

  model.aero = winga;

  wing = mfs_new(fid, model);

# Compute and export results

  wing = mfs_splines(wing);

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

  mfs_export("aero.pos", "msh", winga, "statresp",
             "pressure", "disp");
  mfs_export("solid.pos", "msh", wings, "statresp",
             "disp", "stress", "resultant");

  mfs_merge("solid.pos", "aero.pos", "wing.dsp", "msh");

  mfs_print(fid, wings, "statresp", "resultant");

# Compare pressure in wing sections with rigid results

  load rigid.plt

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

  pfi = pfi / qdyn; pfo = pfo / qdyn;
  xi  = xi / cr;    xo  = xo / cr;
  y   = y * 1e-3;

  for m = 1 : 2

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

     figure(m, "position", [100, 500, 750, 500],
               "paperposition", [0, 0, 17, 12]);
     subplot(2, 2, 1)
        plot(xi(:, 1), squeeze(pri(:, m, 1)),
             xi(:, 1), squeeze(pfi(:, m, 1)));
        legend("rigid", "flexible", "location", "southwest");
        title(header{1});
        grid;
        ylim([-1, 1]);
        ylabel('\Delta c_P');
     subplot(2, 2, 2)
        plot(xi(:, 2), squeeze(pri(:, m, 2)),
             xi(:, 2), squeeze(pfi(:, m, 2)));
        title(header{2});
        grid;
        ylim([-1, 1]);
     subplot(2, 2, 3)
        plot(xo(:, 1), squeeze(pro(:, m, 1)),
             xo(:, 1), squeeze(pfo(:, m, 1)));
        title(header{3});
        grid;
        ylim([-1, 1]);
        xlabel('x/c'); ylabel('\Delta c_P');
     subplot(2, 2, 4)
        plot(xo(:, 2), squeeze(pro(:, m, 2)),
             xo(:, 2), squeeze(pfo(:, m, 2)));
        title(header{4});
        grid;
        ylim([-1, 1]);
        xlabel('x/c');
     print(["cp", num2str(m), EXT], FORMAT);

  endfor

# Stresses in beams

  load data.bin

  for n = 2 : 3
      Upper = sprintf("Stringer_%d_Upper", n); 
      Lower = sprintf("Stringer_%d_Lower", n); 
      elemU = mfs_getset(wings, "eset", Upper);
      elemL = mfs_getset(wings, "eset", Lower);
      y1 = 0.5 * stringer(n).b;
      zo = data.(Upper).geom.P(2);
      zu = zo + stringer(n).h - 0.5 * stringer(n).t;
      coor = [y1, zo; -y1, zo; 0, zu];
      mfs_beamstress(fid, wings, "statresp", [elemU(1), elemL(1)], 
                     coor, 2);
  end

  fclose(fid);
