# Example: Box Beam
#          Transient response due to base motion
#
# -------------------------------------------------------------------

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

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

# Output device depending on operating system

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

# Define excitation parameters

  T    = 0.02;   % Duration of prescribed motion

# Define simulation parameters

  ntip =   10;   % Number of time intervals per period
  nt1  =    5;   % Simulation time as multiple of longest period

# Get results from normal modes analysis

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

  load modes.bin

# Degrees of freedom for output

  rid  = {"Point_A", 3; "Point_B", 3};

# Get simulation time and time step from natural frequencies

  f  = mfs_getresp(beam, "modes", "freq");
  Ts = ceil(100 * nt1 / f(1)) / 100;
  fm = 1000 * ceil(f(end) / 1000);
  dt = 1 / (ntip * fm);

  fprintf(fid, "\nSimulation parameters:\n\n");
  fprintf(fid, "  Simulation time = %12.6e s \n", Ts);
  fprintf(fid, "  Time step       = %12.6e s\n\n", dt);

# Define enforced motion

  motdef = struct("lc", 1, "func", "motion", "params", T);

# Perform modal transient response analysis

  beam = mfs_transresp(beam, [dt, Ts], "load", motdef);

# Plot results at points A and B

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

  figure(1, "position", [100, 200, 1000, 800],
            "paperposition", [0, 0, 15, 15]);
  subplot(3, 1, 1)
     plot(t, u);
     legend("Point A", "Point B",
            "location", "southeast");
     grid;
     axis("labely");
     xlim([0, Ts]);
     ylabel('u [mm]');
  subplot(3, 1, 2)
     plot(t, v);
     grid;
     axis("labely");
     xlim([0, Ts]);
     ylabel('v [m/s]');
  subplot(3, 1, 3)
     plot(t, a);
     grid;
     xlim([0, Ts]);
     ylabel('a [g]'); xlabel('t [s]');
  print(["tresp", EXT], FORMAT);

# Save results

  dlmwrite("tresp.csv", [t; u; v; a], "precision", 5);

# Plot stresses in elements C and D

  sig  = mfs_getresp(beam, "transresp", "stress",
                     {"Elem_C", "Elem_D"});
  sigx = [sig{1}.sigxu; sig{2}.sigxu];

  figure(2, "position", [300, 200, 1000, 750],
            "paperposition", [0, 0, 15, 10]);
  plot(t, sigx);
  legend("Elem. C", "Elem. D");
  grid;
  xlim([0, Ts]);
  ylabel('\sigma_x [MPa]'); xlabel('t [s]');

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

# Perform backtransformation and export displacements

  tb = 0 : 5 * dt : Ts;
  beam = mfs_back(beam, "transresp", "disp", 1, tb);
  mfs_export("tresp.dsp", "msh", beam, "transresp", "disp");

# Save component for restart

  save -binary tresp.bin beam

  fclose(fid);
