# Example: Instrument panel
#          Random response analysis using extracted transfer
#          functions
#
# -------------------------------------------------------------------

  colors = [1, 0, 0; 0, 1, 0; 0, 0, 1];
  set(0, "defaultaxescolororder", colors);

# Data

  g  =    9.81;   % Gravity acceleration in m/s^2
  a0 = 0.3 * g;   % Acceleration level

# Output device depending on operating system

  if (isunix)
     DEVICE = "svg";
  else
     DEVICE = "png";
  endif
  FORMAT = ["-d", DEVICE];
  EXT    = [".", DEVICE];

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

# Results from frequency response analysis

  load freqresp.bin

# Frequencies, power and cross spectral densities

  M   = dlmread("Gxx.csv");
  f   = M(1, :);
  Gxx = M(2 : end, :);
  Gxy = dlmread("Gxy.csv");

# Collect transfer functions

  rid = {"CM_Altimeter", 1; "CM_AirSpeed",  1;
         "CM_Vario",     1; "CM_Altimeter", 3;
         "CM_AirSpeed",  3; "CM_Vario",     3};

  for lc = 1 : 4
      T{lc} = mfs_getresp(panel, "freqresp", "acce", rid, lc);    
  endfor

  for n = 1 : 6
      for lc = 1 : 4
          H{n}(lc, :) = T{lc}(n, :);
      endfor
  endfor

# Compute response spectra

  for n = 1 : 6
      G(n, :) = mfs_randresp(Gxx, Gxy, H{n});
  endfor

# Plot input spectra

  figure(1, "position", [100, 100, 1000, 500],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     semilogy(f, Gxx([1, 3], :))
     title("Longitudinal");
     grid;
     xlabel('f [Hz]'); ylabel('G [m^2/(s^4 Hz)]');
  subplot(1, 2, 2);
     semilogy(f, Gxx([2, 4], :))
     title("Vertical");
     legend("Right", "Left");
     grid;
     xlabel('f [Hz]');
  print(["InputG", EXT], FORMAT);

# Plot response spectra

  figure(2, "position", [100, 300, 1000, 500],
            "paperposition", [0, 0, 17, 10]);
  subplot(1, 2, 1);
     semilogy(f, G(1 : 3, :))
     title("Longitudinal");
     legend("Altimeter", "Air-Speed Indic.", "Variometer",
            "location", "southeast");
     grid;
     xlabel('f [Hz]'); ylabel('G [m^2/(s^4 Hz)]');
  subplot(1, 2, 2);
     semilogy(f, G(4 : 6, :))
     title("Vertical");
     grid;
     xlabel('f [Hz]');
  print(["ResponseG_1", EXT], FORMAT);

# Statistical parameters of excitation

  psix = trapz(f, Gxx, 2);
  RMSx = sqrt(psix);
  N0x  = sqrt(trapz(f, f.^2 .* Gxx, 2) ./ psix);
  arel = a0 ./ RMSx;
  expa = exp(-0.5 * arel.^2);
  Nax  = 2 * expa .* N0x;
  Px   = 100 * erfc(arel / sqrt(2));

  fprintf(fid, "Statistical Parameters of Excitation\n\n"); 
  fprintf(fid, "             right    left\n");
  fprintf(fid, "  ---------------------------------\n");
  fprintf(fid, "  RMS a_x   %7.4f  %7.4f  m/s²\n", RMSx([1, 3]));
  fprintf(fid, "      a_z   %7.4f  %7.4f  m/s²\n", RMSx([2, 4]));
  fprintf(fid, "  N_0 a_x   %7.4f  %7.4f  1/s\n",  N0x([1, 3]));
  fprintf(fid, "      a_z   %7.4f  %7.4f  1/s\n",  N0x([2, 4]));
  fprintf(fid, "  N_a a_x   %7.4f  %7.4f  1/s\n",  Nax([1, 3]));
  fprintf(fid, "      a_z   %7.4f  %7.4f  1/s\n",  Nax([2, 4]));
  fprintf(fid, "  P   a_x   %7.4f  %7.4f  %%\n",   Px([1, 3]));
  fprintf(fid, "      a_z   %7.4f  %7.4f  %%\n",   Px([2, 4]));

# Statistical parameters of response

  psiy = trapz(f, G, 2);
  RMSy = sqrt(psiy);
  N0y  = sqrt(trapz(f, f.^2 .* G, 2) ./ psiy);
  arel = a0 ./ RMSy;
  expa = exp(-0.5 * arel.^2);
  Nay  = 2 * expa .* N0y;
  Py   = 100 * erfc(arel / sqrt(2));

  fprintf(fid, "\nStatistical Parameters of Response\n\n"); 
  fprintf(fid, "              Alti.   Speed    Vario\n");
  fprintf(fid, "  ------------------------------------------\n");
  fprintf(fid, "  RMS a_x   %7.4f  %7.4f  %7.4f  m/s²\n",
          RMSy(1 : 3));
  fprintf(fid, "      a_z   %7.4f  %7.4f  %7.4f  m/s²\n",
          RMSy(4 : 6));
  fprintf(fid, "  N_0 a_x   %7.4f  %7.4f  %7.4f  1/s\n",
          N0y(1 : 3));
  fprintf(fid, "      a_z   %7.4f  %7.4f  %7.4f  1/s\n",
          N0y(4 : 6));
  fprintf(fid, "  N_a a_x   %7.4f  %7.4f  %7.4f  1/s\n",
          Nay(1 : 3));
  fprintf(fid, "      a_z   %7.4f  %7.4f  %7.4f  1/s\n",
          Nay(4 : 6));
  fprintf(fid, "  P   a_x   %7.4f  %7.4f  %7.4f  %%\n",
          Py(1 : 3));
  fprintf(fid, "      a_z   %7.4f  %7.4f  %7.4f  %%\n",
          Py(4 : 6));

  fclose(fid);
