%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Script to calculate the 2d seismic response of a rigid block                                                                                                                                                                                    %                                                                                                                                                                                              
%  created by Michalis Vassiliou in 2008 and evolved by him since then                                                                                                                                                                                   %     
%  mfvassiliou@gmail.com                                                                                                                                                                                                                           %     
%                                                                                                                                                                                                                                                  %                                         
%  Feel free to use it at will, as long as you don't make money out of it. In the latter case, the ETH will not be happy about it.                                                                                                                                % 
%                                                                                                                                                                                                                                                  %                                         
%  Please report any bugs/suggestions to mfvassiliou@gmail.com                                                                                                                                                                                                             %     
%                                                                                                                                                                                                                                                  %                         
%  We have used this script and variations of it in the following papers:                                                                                                                                                                           %
%                                                                                                                                                                                                                                                  %
%  Vassiliou MF, Mackie KR, and Stojadinovi? B. (2016) A finite element model for seismic response analysis of deformable rocking frames. Earthquake Engineering & Structural Dynamics                                                             % 
%  Vassiliou MF, Truniger R, and Stojadinovi?, B. (2015). An analytical model of a deformable cantilever structure rocking on a rigid surface: development and verification. Earthquake Engineering & Structural Dynamics, 44(15), 2775-2794.      % 
%  Truniger R, Vassiliou MF, and Stojadinovi?, B. (2015). An analytical model of a deformable cantilever structure rocking on a rigid surface: experimental validation. Earthquake Engineering & Structural Dynamics, 44(15), 2795-2815.           %     
%  Vassiliou MF and Makris N. (2015). Dynamics of the Vertically Restrained Rocking Column. Journal of Engineering Mechanics, 141(12)                                                                                                              % 
%  Makris, N., & Vassiliou, M. F. (2014). Dynamics of the rocking frame with vertical restrainers. Journal of Structural Engineering, 141(10)                                                                                                      % 
%  Makris N and Vassiliou MF (2014), Are some top-heavy structures more stable? Journal of Structural Engineering, 140(5)                                                                                                                          % 
%  Vassiliou MF, Mackie KR and Stojadinovic B (2014), Dynamic Response Analysis of Solitary Flexible Rocking Bodies: Modeling and Behavior under Pulse-Like Ground Excitation, Earthquake Engineering and Structural Dynamics , 43: 1463–1481      % 
%  Makris N and Vassiliou MF (2012), Planar rocking response and stability analysis of an array of free-standing columns capped with a freely supported rigid beam, Earthquake Engineering and Structural Dynamics, 41(2), 177-196                 % 
%  Makris N and Vassiliou MF (2012), Sizing the slenderness of free-standing rocking columns to withstand earthquake shaking, Archive of Applied Mechanics, 1-15.                                                                                  % 
%  Vassiliou MF and Makris N (2011), Analysis of the rocking response of rigid blocks standing free on a seismically isolated base. Earthquake Engineering and Structural Dynamics, 41(2), 177-196                                                 %                                                                     % 
%                                                                                                                                                                                                                                                  %             
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%alpha = block slenderness
%R = block semidiagonal
%y(:,1) = theta; y(:,2) = thetadot 


close all
clear all

record_dt=0.02;% depends on the ground motion


%Block Characteristics
alpha=atan(1/10); %block slenderness alpha
R=3; %Block semidiagonal R
p=sqrt(3*9.81/4/R); %frequency parameter p=sqrt(3g/4R)
v_lose=1-3/2*(sin(alpha))^2; %coefficient of restitution. I used the Housner one, even though it is not validted experimentally.

%ground motion
ug=textread('elcentro.txt'); %read ground motion
ug=ug*9.81; %make ground motion it in m/s^2 
time_r=record_dt*[1:length(ug)]; 
tspan=time_r(1:end-1);

tic

%initialize
tstart=0;
tout=[0]; %in tout I store the t of the solution
yout=[0 0]; % in yout I store the solution. First column is theta, second column is thetadot

ti=1;
tfinal = time_r(end-1); %we will solve the equation up until tfinal. 
t_land_index=1;


while 1>0 %start of outside loop. I use this loop because the block might uplift, rock, eventually stop rocking (check the next loop for that), and then another spike might uplift it again.

 %find when rocking initiates
 t_up_index=find(abs(ug(t_land_index:end))/tan(alpha)/9.81>1,1,'first')+t_land_index-1;
 tstart=t_up_index*record_dt;

 if size(t_up_index)~=0 %acceleration is enough to uplift
    %padd with zeros the time period that we had no rocking 
    yzeros=zeros(t_up_index-t_land_index,2); 
    tzeros=[record_dt:record_dt:size(yzeros,1)*record_dt]';
    tzeros=tzeros+tout(end);
    yout=[yout;yzeros];
    tout=[tout;tzeros];
    %end of padding with zeros
    
    tstart=t_up_index*record_dt; %this is the time for which we had uplift
    y0 = [0 0]; %set initial condition for the rocking motion. It uplifts with zero theta and thetad
    refine = 4; 
    teout=[]; %in this vector I store the time instants at which an event happens
    yeout=[]; %in this vevtor I store the y values when event happened
    ieout=[]; %in this vector I sotre the indexes that tell me what event happened

    Ec=3/8/p^2*y0(2)+3/4*(cos(alpha-y0(1))-cos(alpha)); %this is the normalized total energy of the system. 
    Emax=Ec;
    
    clear options
    options = odeset('Events',@events_rock_stop_ug,'Refine',refine,'RelTol',1e-12,'AbsTol',1e-14*[1 1],'maxstep',1e-3); % CAUTION: Rocking is extremely non-linear. Always play with tolerances an check convergence!!!

    while (Ec>=1e-3*Emax)&&(tout(end)<0.95*tfinal)%i.e. while it rocks. This while repeats itself for each rocking cycle. 
        % I stop the loop if the total energy of the system becomes smaller
        % than 1e-3 of the maximum energy it had during the "outside while
        % loop" in which we currnetly are. Otherwise it would run forever,
        % because it would stop after infinite number of impacts.
   

        [t,y,te,ye,ie] = ode45(@rock,[tstart:record_dt:tfinal],y0,options,alpha,p,time_r,record_dt,ug); %solve the rocking ode up until we have impact or overturning

        t_land_index=floor(t(end)/record_dt)+1; %find the index (in the time_r vector) of the time instant at which we had an event
        
        % Accumulate output.  
        nt = length(t);
        tout = [tout; t(2:nt)];
        yout = [yout; y(2:nt,:)];
        teout = [teout; te];    % Events at tstart are never reported.
        yeout = [yeout; ye];
        ieout = [ieout; ie];


        if ie==1 %ode solver stopped because theta became 0: we have impact.
            % Set the new initial conditions.
            y0(1) = 0;
            y0(2) = v_lose*y(nt,2);
            tstart = t(nt);
            Ec=3/8/p^2*(y0(2))^2+3/4*(cos(alpha-y0(1))-cos(alpha)); %set the energy at uplift. At that instant we only have kinetic energy.
            if Ec>Emax
                 Emax=Ec; %if this energy is larger that the up-to-now maximum set it as maximum.
            end

        else %it has fallen so I need to break
        break
        end
    end %inside while

    if ie(end)==2 %it has fallen so I need to break the outside loop too
    break
    end

 else %acceleration is not enough to uplift so I need to break the outside loop
    tout(end+1)=time_r(end);
    yout(size(yout,1)+1,:)=[0 0];
    break
 end %end of if statement that checks if the acceleration is enough to uplift

 tstart=tout(end);
     if tstart>time_r(end-1) %it means that it would be pointless to try to find any uplift points from now on
         break
     end
end

toc

figure(1)
hold on
subplot(2,1,1)
plot(tout,yout(:,1)/alpha)
hold on
grid on
ylabel('$\frac{\theta(t)}{\alpha}$','interpreter','latex','fontsize',16)

subplot(2,1,2)
plot(time_r,ug)

