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


global Mi K C b h

displ = x(1:6);     % displacement relative to ground
veloc = x(7:12);    % velocity relative to ground
z     = x(13);      % dimenionless isolation force

Dy   = param(1);       % yield displacement for the isolation bearing
Fy   = param(2);       % yield force for the isolation bearing
Dc   = param(3);       % contact displacement
Kc   = param(4);       % contact stiffness

Fc = 0.0;              % force due to contact/impact
if (displ(1) > Dc)     % force applies if the displacement exceeds the 
                       % isolation contact displacement
    Fc = Kc*(displ(1) - Dc);
end
if (displ(1) < -Dc)    % force applies if the displacement exceeds the
                       % isolation contact displacement in the opposite
                       % direction
    Fc = Kc*(displ(1) + Dc);
end

Fi = Fy*z;            % base isolation force

Fnl = b'*(Fc+Fi);     % nonlinear forces :
                      % base isolation force and the isolation contact force

accel  = Mi*(-K*displ - C*veloc - Fnl) - h*quake_accel;
z_dot  = (1.0 - z^2*(0.5*sign(z*veloc(1)) + 0.5))*veloc(1)/Dy;

x_dot = [ veloc ;      % velocity
          accel ;      % acceleration
          z_dot ];     % time derivative of the dimensionless isolation force

% ----------------------------------------- example4.m  January 23 2006
