# Example: Trim analysis of a standard class glider
#          Determine the position of the neutral point
#
# ---------------------------------------------------------------

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

# Data (kg, m, s)

  v   = [25, 50];    % Flight velocities
  rho = 1.21;        % Mass density of the air
  g   = 9.81;        % Gravity acceleration
  xs  = [1.7, 1.9];  % x-coordinate of center of mass

# Model definition

  model = glider_aero();

# Configurations

  qdyn     = 0.5 * rho * v.^2;
  cname{1} = sprintf("v = %2.0f m/s", v(1));
  cname{2} = sprintf("v = %2.0f m/s", v(2));

  model.config = struct("name", cname, "qdyn", num2cell(qdyn),
                        "ay", 0, "az", g,
                        "pacce", 0, "yacce", 0, "racce", 0,
                        "pitch", 0, "yaw", 0, "roll", 0);

# Trim analysis for two different positions of center of mass

  for n = 1 : 2

      fprintf(fid, "xS = %4.1f m\n\n", xs(n));

      model.mass.cm(1) = xs(n);

      glider = mfs_new(fid, model);
      glider = mfs_trim(glider);
      mfs_print(fid, glider, "trim", "params");

      tp(n) = mfs_getresp(glider, "trim", "params");

      fprintf(fid, "\n");
      for k = 1 : 8
          fprintf(fid, "==========");
      endfor
      fprintf(fid, "\n");

  endfor

# Determination of neutral point

  elev = reshape([tp.elevator], 2, 2);
  qinv = 1 ./ qdyn';

  a  = (elev(2, :) - elev(1, :)) ./ (qinv(2) - qinv(1));
  b  = (a(2) - a(1)) / (xs(2) - xs(1));
  xN = xs - a / b;

  fprintf(fid, "\na = [%7.3f, %7.3f] Pa\n", a);
  fprintf(fid, "b = %7.2f Pa/m\n", b);
  fprintf(fid,
     "\nPosition of neutral point: xN = %7.3f %7.3f m\n", xN);

  fclose(fid);
