# Airfoil test:
#
#   Airfoil: FX 66-S-196 V1
#   Units:   mm
#
# The airfoil data are approximated by a spline. The approximation is
# tested by computing the pressure coefficient using the discrete
# vortex method.
#
# --------------------------------------------------------------------

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

  set(0, "defaultaxesfontsize", 10);

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

# Data

  ns    =   5;            % No. of spline segments
  ng    = 100;
  nv    =  40;            % No. of discrete vortices
  alpha = [-2, 0, 2, 5];  % Angles of attack

# Read camber from csv file

  data   = dlmread("fx66s196v1.csv", ",", "A100:B143");
  c      = data(end, 1);
  data   = data / c;
  [m, n] = size(data);

# Generate the spline approximation

  zspp = mfs_airfoil("fit", "camber", data, ns);

# Define interval boundaries

  xg = linspace(0, 1, ng + 1);

# Graphical check of interpolation

  zsa = ppval(zspp, xg);
  dzs = ppval(ppder(zspp), xg);

  figure(1, "position", [100, 600, 750, 750],
            "paperposition", [0, 0, 14, 12]);
  subplot(2, 1, 1)
     plot(data(:, 1), 100 * data(:, 2), "color", "green", 
          "marker", "o", "linestyle", "none",
          xg, 100 * zsa, "color", "red");
     legend("Data", "Spline", "location", "south");
     grid;
     ylabel("z/c [%]");
  subplot(2, 1, 2)
     plot(xg, dzs, "color", "red");
     grid;
     ylabel("dz/dx"); xlabel("x/c");
  print(["fx_z", EXT], FORMAT);

# Compute aerodynamic coefficients

  x = linspace(0, 1, nv + 1);
  [cp, xv, cL, cM] = mfs_vortex2d(x, zspp, alpha);

# Output results

  fprintf(fid, "Aerodynamic coefficients:\n\n");
  fprintf(fid, "  alpha      cL        cM\n");
  fprintf(fid, "  ---------------------------\n")
  for n = 1 : length(alpha)
      fprintf(fid, "  %5.2f   %7.4f   %7.4f\n", 
              alpha(n), cL(n), cM(n));
      txt{n} = sprintf("\\alpha = %4.1f\\deg", alpha(n));
  endfor

  figure(2, "position", [200, 500, 750, 500],
            "paperposition", [0, 0, 12, 8]);
  plot(xv, cp);
  legend(txt, "location", "northeast");
  grid;
  ylim([-1, 5]);
  xlabel('x/c'); ylabel('\Delta c_p');
  print(["fx_cp", EXT], FORMAT);

  fclose(fid);
