function rmax = example3_sim(Z);
% Example 3 - Dynamic Problem of One Degree of Freedom
% (The function 'example3_sim' in this file is called by 'g.m')
%
% Input:
%   rmax - Maximum displacement over time
% Output:
%   Z    - Sample values of the normal random variables
%
% Department of Civil and Environmental Engineering
% Duke University
% Siu Chung Yau, Henri P. Gavin, January 2006


% statistical properties of the lognormal random variables
% (mu = mean     V = coefficient of variance)      
mu(1) = 1000;   V(1) = 0.05;       % m  [mass, kg]
mu(2) = 4000;   V(2) = 0.10;       % k  [stiffness, N/m]
mu(3) = 250;    V(3) = 0.20;       % c  [damping coeff, N/m/s]

mu(4) = 400;    V(4) = 0.10;       % fy [yield force, N]
mu(5) = 0.01;   V(5) = 0.15;       % ry [yield displacement, m]

mu(6) = 400;    V(6) = 0.20;       % Fp [force ampltude, N]
mu(7) = 0.5;    V(7) = 0.50;       % Tp [force period, sec]

% the sample values of normal random variables are transformed into
% lognormal random variables
sigma_logx = sqrt(log(V.^2+1));         
mu_logx = log(mu) - 0.5*sigma_logx.^2;  
X = exp(mu_logx + Z.*sigma_logx);      

% 7 lognormal random variables
m  = X(1);      % mass
k  = X(2);      % stiffness
c  = X(3);      % damping coefficient
fy = X(4);      % yield force
ry = X(5);      % yield displacement
Fp = X(6);      % force amplitude
Tp = X(7);      % force period

delta_t = 0.005;            % time step (sec)
points = 1000;              % number of time steps

tp = floor(Tp/delta_t);     % time points in the pulse
if points < 2*tp 
  points = 2*tp;            % time points in the pulse depends on the
                            % random variable Tp (force period)
                            %
                            % increase the number of time steps if the time
                            % points in the pulse is too many
end;   

fext = zeros(1,points);     % initialize external force vector

t = [0:points-1]*delta_t;   % time axis

fext(1:tp) = Fp*0.5*(1-cos(2*pi*[1:tp]*delta_t/Tp));    % full cosine-wave
                                                        % external force 
                                                        % impulse

param = [m k c fy ry];

x0 = zeros(3,1);    % initial condition of the system

[t,x] = ode4('example3', t, x0, fext, param);    % solve the problem

displ = x(1,:);     % displacement over time
veloc = x(2,:);     % velociety over time
fnl   = fy*x(3,:);  % nonlinear hysteretic force over time
accel = (-k.*displ - c.*veloc - fnl + fext)./m; % acceleration over time

rmax = max(abs(displ)); % maximum dieplacement over time

% end ------------------------------------- example3_sim.m  January 23 2006
