function g = example4_sim(Z);
% Example 4 - Dynamic Problem of Six Degrees of Freedom
% (The function 'example4_sim' in this file is called by 'g.m')
%
% Department of Civil and Environmental Engineering
% Duke University
% Siu Chung Yau, Henri P. Gavin, January 2006


global Mi K C b h data

% statistical properties of the lognormal random variables
% (mu = mean     V = coefficient of variance)      

mu(1) = 6000;       V(1) = 0.10;   % m_fl (mass of a floor, kg)
mu(2) = 30000000;   V(2) = 0.10;   % k_fl (stiffness of a floor, N/m)
mu(3) = 60000;      V(3) = 0.20;   % c_fl (damping coeff of a floor, N/m/s)

mu(4) = 0.05;       V(4) = 0.20;   % Dy (isolation yield displacement, m)
mu(5) = 20000;      V(5) = 0.20;   % Fy (isolation yield force, N)
mu(6) = 30000000;   V(6) = 0.30;   % Kc (isolation contact force, N/m)

mu(7) = 1.0;        V(7) = 0.20;   % T (pulse excitation, sec)
mu(8) = 1.0;        V(8) = 0.50;   % A (pulse amplitude, m/s/s)

% sample values of normal random variabels are converted to those of
% 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);

% 8 lognormal random variables
m_fl = X(1);               % mass of a floor, N/m/s/s
k_fl = X(2);               % stiffness of a floor, N/m
c_fl = X(3);               % damping coefficient

Dy = X(4);                 % isolation yield displacement, m
Fy = X(5);                 % isolation yield force, N
Dc = 0.5;                  % isolation contact displacement, m
Kc = X(6);                 % isolation contact stiffness, N/m

T = X(7);                  % pulse excitation, sec
A = X(8);                  % pulse amplitude, m/s/s

m1 = 500;                  % mass of block 1, N/m/s/s
m2 = 100;                  % mass of block 2, N/m/s/s

k1 = 2500;                 % stiffness of spring 1, N/m
k2 = 100000;               % stiffness of spring 2, N/m

c1 = 350;                  % damping coefficient of damper 1, N/m/s
c2 = 200;                  % damping coefficient of damper 2, N/m/s

M = [ m_fl    0    0    0    0    0 ;                 % mass matrix
         0 m_fl    0    0    0    0 ;
         0    0 m_fl    0    0    0 ;
         0    0    0 m_fl    0    0 ;
         0    0    0    0   m1    0 ;
         0    0    0    0    0   m2 ];
  
K = [ k_fl     -k_fl      0     0      0    0 ;       % stiffness matrix
     -k_fl 2*k_fl+k1  -k_fl     0    -k1    0 ;
         0     -k_fl 2*k_fl -k_fl      0    0 ;
         0         0  -k_fl  k_fl      0    0 ;
         0       -k1      0     0  k1+k2  -k2 ;
         0         0      0     0    -k2   k2 ];

C = [ c_fl     -c_fl      0     0      0    0 ;       % damping matrix
     -c_fl 2*c_fl+c1  -c_fl     0    -c1    0 ;
         0     -c_fl 2*c_fl -c_fl      0    0 ;
         0         0  -c_fl  c_fl      0    0 ;
         0       -c1      0     0  c1+c2  -c2 ;
         0         0      0     0    -c2   c2 ];

Mi = inv(M);

b = [ 1 0 0 0 0 0 ]';	% location of nonlinear isolator force
h = [ 1 1 1 1 1 1 ]';   % location of inertial earthquake loads
 
param(1) = Dy;        % yield displacement for the isolation bearing
param(2) = Fy;        % yield force for the isolation bearing
param(3) = Dc;        % isolation contact displacement
param(4) = Kc;        % isolation contact stiffness

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

tp = floor(T/delta_t);      % time points in the pulse

if points < 3*tp 
  points = 3*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;    

t = [0:points-1]*delta_t;   % time axis
   
u = zeros(1,points);        % initialize acceleration vector

u(1:tp) = A*sin(2*pi*[1:tp]*delta_t/T); % accleration is a full sine-wave

x0 = zeros(13,1);           % initial system condition

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

displ = x(1:6,:);           % displacement relative to ground
veloc = x(7:12,:);          % velocity relative to ground
Fi    = Fy*x(13,:);         % dimensionless isolation force

Fc = zeros(1,points);       % initialize contact force vector
for p=1:points
    Fc(p) = 0.0;            
    if (displ(1,p) > Dc)    % force applies if the displacement exceeds the 
                            % isolation contact displacement
        Fc(p) = Kc*(displ(1,p) - Dc);
    end
    if (displ(1,p) < -Dc)   % force applies if the displacement exceeds the 
                            % isolation contact displacement in the
                            % opposite direction
        Fc(p) = Kc*(displ(1,p) + Dc);
    end
end

Fnl = [ 1 0 0 0 0 0 ]'*(Fc+Fi);     % all nonlinear forces

Co1 = [ -1  1  0  0  0  0 ;
         0 -1  1  0  0  0 ;
         0  0 -1  1  0  0];

Co2 = [ -1  0  0  0  1  0];
        
a = Mi*(-C*veloc-K*displ-Fnl);  % acceleration

g = 12.5*(0.04-max(max(abs(Co1*displ)))) + (0.5-max(abs(a(6,:)))) + ...
    2*(0.25-max(abs(Co2*displ)));   % limit state value

% ------------------------------------- example4_sim.m  January 23 2006
