# Example: Plate
#          Modal transient response analysis
#          Initial deformation
#
# -------------------------------------------------------------------

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

# Output device depending on operating system

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

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

# Define excitation parameters

  T      =  2e-3;  % Duration of decay

# Define simulation paramters

  dt     =  1e-4;  % Time step
  Ts     =  0.15;  % Simulation time

# Define excitation

  loads = struct("lc", 1, "func", "decay", "params", T); 

# Get results from normal modes analysis

  if (! isfile("modes.bin"))
     error("File modes.m must be executed before\n");
  end

  load modes.bin

# Compute initial deformation

  plate = mfs_statresp(plate);
  plate = mfs_newset(plate, "nset", "union",
                     {"Resp_A", "Resp_B", "Resp_C"}, "Resp");
  ws    = mfs_getresp(plate, "statresp", "disp", "Resp");

# Perform transient response analysis

  plate = mfs_transresp(plate, [dt, Ts], "load", loads, 
                        "statresp", 1);

# Plot results

  rid = {"Resp_A", 3; "Resp_B", 3; "Resp_C", 3};

  t = mfs_getresp(plate, "transresp", "time");
  w = mfs_getresp(plate, "transresp", "disp", rid);
  v = mfs_getresp(plate, "transresp", "velo", rid) / 1000;
  a = mfs_getresp(plate, "transresp", "acce", rid) / 9810;

  figure(1, "position", [100, 200, 1000, 750],
            "paperposition", [0, 0, 15, 12]);
  plot(t, w);
  legend("Pnt. A", "Pnt. B", "Pnt. C",
         "location", "southeast");
  grid;
  xlim([0, Ts]);
  xlabel('t [s]'); ylabel('w [mm]');

  print(["decay_modal_w", EXT], FORMAT);

  figure(2, "position", [300, 200, 1000, 750],
            "paperposition", [0, 0, 15, 12]);
  plot(t, v);
  legend("Pnt. A", "Pnt. B", "Pnt. C",
         "location", "southeast");
  grid;
  xlim([0, Ts]);
  xlabel('t [s]'); ylabel('v [m/s]');

  print(["decay_modal_v", EXT], FORMAT);

  figure(3, "position", [500, 200, 1000, 750],
            "paperposition", [0, 0, 15, 12]);
  plot(t, a);
  legend("Pnt. A", "Pnt. B", "Pnt. C",
         "location", "southeast");
  grid;
  xlim([0, Ts]);
  xlabel('t [s]'); ylabel('a [g]');

  print(["decay_modal_a", EXT], FORMAT);

# Save results

  dlmwrite("decay_modal.csv", [t; w; v; a], "precision", 5);

# Check accuracy of initial conditions

  ew = 100 * (ws(:, 3) - w(:, 1)) ./ ws(:, 3);
  fprintf(fid, "  Displacement error:\n");
  fprintf(fid, "     Point A: %11.4e %%\n", ew(1));
  fprintf(fid, "     Point B: %11.4e %%\n", ew(2));
  fprintf(fid, "     Point C: %11.4e %%\n", ew(3));

# Compute and export deformed shapes

  tb = t(1 : 5 : length(t));
  plate = mfs_back(plate, "transresp", "disp", 1, tb);
  mfs_export("decay_modal.dsp", "msh", plate,
             "transresp", "disp");

  fclose(fid);
