% Script based on the scripts used for 
% Vassiliou, M. F., Burger, S., Egger, M., Bachmann, J. A., Broccardo, M., & Stojadinovic, B. (2017). The three?dimensional behavior of inverted pendulum cylindrical structures during earthquakes. Earthquake Engineering & Structural Dynamics, 46(14), 2261-2280.
% Vassiliou, M. F. (2018). Seismic response of a wobbling 3D frame. Earthquake Engineering & Structural Dynamics, 47(5), 1212-1228.
% Rajah R (2016), Dynamic response of a rectangular plate supported on four cylindrical columns allowed to uplift and wobble, Project Thesis, ETH Zurich 
% Final coding by Michalis F. Vassiliou 
% Use freely and cite if you think the papers are interesting
% Please report bugs (that exist for sure) to mfvassiliou@gmail.com
% 22.08.2018

function [par]=rocking_podium_structure_TH_alphaR
clear
close all
%total height of the columns
h=9.6;
%diameter of the columns
b=1.6;
%column slenderness
alpha=b/h;
%mass ratio: slab over total mass of the columns
gamma=10;

load ag_elcentro.mat

R=sqrt(h^2+b^2)/2;
% Semidiagonal length R: In order to construct the overturning spectra,
% define the range of R. In order to investigate a specific podium
% structure, just insert the semidiagonal length. 
% R = [0.5:0.5:60]';

% Range of scale factor for the overturning spectra. 
scale=1;
% sscale = [1.2:0.1:12]';
%% INPUT
% Geometrical properties
par.g = 9.807;
par.R = R;
par.alpha = alpha;
par.gamma = gamma;
par.phat = sqrt(12*(2*gamma+1)/(48*gamma+15+cos(alpha)^2)*par.g/par.R); 

% Start time, end time, timestep 
time.start=ag(1,1);
time.end=ag(end,1);
time.dt=ag(2,1)-ag(1,1);

% Excitation seperated in its compenent in u and v direction.
% Time vector and timestep given by the ag-file (ElCentro or Kobe)
excitation.ugdd=ag(:,2)*par.g;
excitation.vgdd=ag(:,3)*par.g;
excitation.t=ag(:,1);
excitation.dt=time.dt;

%% INITIALIZATION 
% Define the time vector with the time range in which the calculations 
% will be done. 
N = length(ag);
time_r=time.dt:time.dt:time.dt*(N+500);

% Define the excitation vector with the same length as the time vector
% Build up the square-mean vector
ugdd=zeros(N+500,1);
vgdd=zeros(N+500,1);
ugdd(1:N,1)=excitation.ugdd;
vgdd(1:N,1)=excitation.vgdd;
ugdd_phi=sqrt(ugdd.^2+vgdd.^2);

% For the overturning spectra: Define the zero matrix for the contour plot
C = zeros(length(scale),length(R));

% Define the initial states
fall=0;
contact=0;
finish=0;

% Standstill time 
tContact=time.dt;
%% UPLIFT CONDITION
% Define the uplift condition
a_need=par.g*tan(alpha);

% Needed excitation for uplift
a_max=max(ugdd_phi);
ugdd_need = ugdd./a_max*a_need;
vgdd_need = vgdd./a_max*a_need;


%% SOLVING PROCEDURE & OUTPUT
% Steps is the number of the calculations. 
% Waitbar showing you how long to wait
% waitbar(0,'Please wait...');
steps = length(R)*length(scale);
step = 0;
l = 1;
% disp('0 % done')

% Set j=1
j = 1;

while step < steps

    % while-loop: R increases step for step from 1.2 to 12
    while j < length(scale)+1

        % Scale the ground motion
%         ugdd = ugdd_need .* scale(j);
%         vgdd = vgdd_need .* scale(j);

        % Set i=1
        i = 1;

        % while-loop: for a given scale factor run the solver
        % for the range of R
        while i < length(R)+1

        % Build up the square-mean vector
        ugdd_phi = sqrt(ugdd.^2+vgdd.^2);


        % Locate the first element where uplift occurs and define the
        % corresponding time
        t_ind=find(abs(ugdd_phi)>=1.01*par.g*tan(alpha),1,'first');
        tRock=time.dt*find(abs(ugdd_phi)>=1.01*par.g*tan(alpha));

        % Fill up the output vector yOut with zeros until the first
        % first uplift occurs, define the corresponding time output
        % vector
        yOut=zeros(t_ind,4);
        tOut=time_r(1:t_ind)';

        % Check if there is an uplift
        if isempty(tRock) == 1
            
            % if not, then the tilt angle is zero
            % set i = i+1 to run the calculation for the next R. 
            C(j,i) = 0;
            i = i+1;
            prog = l/steps*100;
%             clc; disp([num2str(prog) '% done'])
            l = l+1;
            step = step+1;
%             waitbar(step/steps)
        else
            
            % if yes, then...
            while 2>0 
                
                % it might be that after standstill there's anothter
                % uplift event, that is, check this case
                tiUplift=find(tRock>tContact,1,'first');
                
                if isempty(tiUplift) == 1
               
                    %  if there is no other uplift event, stop the
                    %  calculation
                    break
                else
                    
                    % Define the time when uplift occurs.
                    % This corresponts to t_start, duration is limited by
                    % t_final
                    tUplift = tRock(tiUplift);
                    tiUpliftOrig = round(tUplift/time.dt) ;
                    t_start = tUplift;
                    t_final = length(time_r)*time.dt+5;
                    
                    % Define the direction where the first uplift occurs
                    if ugdd(tiUpliftOrig)>0 && vgdd(tiUpliftOrig)>0
                    phi_0=atan(vgdd(tiUpliftOrig)/ugdd(tiUpliftOrig))+pi;
                    elseif ugdd(tiUpliftOrig)<0 && vgdd(tiUpliftOrig)>0
                    phi_0=atan(abs(ugdd(tiUpliftOrig))/vgdd(tiUpliftOrig))+3*pi/2;
                    elseif ugdd(tiUpliftOrig)<0 && vgdd(tiUpliftOrig)<0
                    phi_0=atan(abs(vgdd(tiUpliftOrig))/abs(ugdd(tiUpliftOrig)));
                    else
                    phi_0=atan(ugdd(tiUpliftOrig)/abs(vgdd(tiUpliftOrig)))+1*pi/2;
                    end

                    % The firt motion is assumed to be a planar rocking,
                    % this happens during this time interval
                    t_up = [t_start  t_start+time.dt/10];

                    % Initial condition for the planar rocking
                    IC_0=[0,0];

                    % ODE-solver with the equations of motion: uplift_DGL 
                    [t_1, y_1]=ode45(@(t_1, y_1) uplift_DGL(t_1,y_1,time_r,time,...
                        ugdd_phi,par),t_up,IC_0);
                    theta_0=y_1(end,1);
                    thetad_0=y_1(end,2);

                    % Accumulation Output vector
                    tOut = [tOut; t_1(1)];
                    yOut = [yOut; [0 0 0 0]];

                    % Set the options vor ode45
                    refine = 6;
                    options = odeset('Refine',refine,'RelTol',1e-8,...
                        'absTol',1e-8*[1 1 1 1]);
                    IC = [theta_0 thetad_0 phi_0 0];

                    % For the rest of the time.dt: 3D rocking
                    t_3D = [t_start+time.dt/10 t_start+time.dt];
                    [t,y] = ode45(@(t,y) podium_DGL(y,t,time_r,time,ugdd,vgdd,...
                        par),t_3D,IC,options);

                    % Accumulation Output vector
                    nt= length(t);
                    tOut = [tOut; t(1:nt)];
                    yOut = [yOut; y(1:nt,:)];

                    % Update the time
                    t1=t(nt);

                    % Update the initial conditions
                    IC=y(nt,:);

                    % The iteration is continued until a stop
                    % criterion is fulfilled 
                    while 1>0

                    % Start now the iteration (t=0.1s) until t_final
                    t2 = min(t1+0.1,t_final);
                    
                    % Time intervall
                    t_3D=[t1 t2];
                    
                    % Solver
                    [t,y]=ode45(@(t,y) podium_DGL(y,t,time_r,time,ugdd,vgdd,par)...
                        ,t_3D,IC,options);
                    
                    % Update the time
                    nt = length(t);
                    
                    % Accumulation Output vector
                    tOut = [tOut; t(1:nt)];
                    yOut = [yOut; y(1:nt,:)];

                    % Check if there is an overturn event
                    if max(yOut(:,1)) >= pi/2
                        
                        % Stop the while-loop.  
                        fall = 1;
                        break
                    end
                        % if there is no overturn event, check if the tilt
                        % angle is smaller than R*1e-5
                        if t(end)>t_start+0.3 && max(yOut((end-3*nt):end,1))<...
                                par.R*1e-5
                            
                        % if yes, we have (almost) a standstill, the
                        % calculations will be interrupted and the tContact
                        % is updated. 
                        % Accumulatoion Output vector
                        % Stop this while-loop and check if there is
                        % another uplift event after this standstill.
                        tContact=t(nt);
                        contact = 1;
                        tOut=[tOut; tOut(end)];
                        yOut=[yOut; [0 0 0 0]];

                        break
                        else
                            
                        contact = 0;
                        
                        % Check if we reach the final time or not
                        if t(end)<time_r(end)-2*time.dt
                            
                            % if not update the time and the initial
                            % condition for the next iteration
                            t1=t(nt);
                            IC=y(nt,:);
                            finish = 0;
                            continue
                        else
                            
                            % if yes stop the calculations 
                            finish = 1;
                            break
                        end

                        end



                    end

                end

                % for overtrun event, tilt angle is equal to pi/2
                if fall == 1
                    C(j,i) = pi/2;
                    i = i+1;
                    prog = l/steps*100;
                    clc; disp([num2str(prog) '% done'])
                    l = l+1;
                    step = step+1;
%                     waitbar(step/steps)
                    break
                end

                % if the calculations have been stopped due to the timeout,
                % then one chose the max tilt angle occured during this
                % time range
                if finish == 1
                    C(j,i) = max(yOut(:,1));
                    i = i+1;
                    prog = l/steps*100;
%                     clc; disp([num2str(prog) '% done'])
                    l = l+1;
                    step = step+1;
%                     waitbar(step/steps)
                    break
                end

            end

%         plot tilt angle, velocity of the tilt angle, rolling angle,
%         velocity of the rolling angle against time
        
    


        end

        end

        j = j+1;

    end

end


    figure   

        subplot(6,1,1)
        hold on
        plot(tOut,yOut(:,1))
        xlim([0 45])
        ylabel('$\theta $ [rad]','interpreter','latex')
        subplot(6,1,2)
        hold on
        plot(tOut,yOut(:,2))
        xlim([0 45])
        ylabel('$\dot{\theta} $ [rad/s]','interpreter','latex')
        subplot(6,1,3)
        hold on
        plot(tOut,yOut(:,3))
        xlim([0 45])
        ylabel('$\varphi $ [rad]','interpreter','latex')   
        subplot(6,1,4)
        hold on
        plot(tOut,yOut(:,4))
        xlim([0 45])
        ylabel('$\dot{\varphi} $ [rad/s]','interpreter','latex')
        subplot(6,1,5)
        hold on
        plot(ag(:,1),ag(:,2))
        ylabel('$\ddot{u}_g (g)$','interpreter','latex')
        text(25,0.2,'1940 El Centro Array #9, 180^o')
        xlim([0 45])
        subplot(6,1,6)
        hold on
        plot(ag(:,1),ag(:,3))
        ylabel('$\ddot{u}_g (g)$','interpreter','latex')
        xlabel('Time (s)')
        text(25,0.2,'1940 El Centro Array #9, 270^o')
        ylim([-0.5 0.5])
        xlim([0 45])
        
set(gcf,'PaperType','A4','paperunits','normalized','paperposition',[0 0 0.8 0.8])
% print -depsc TH_ElCentro.eps
% Normalize tilt angle by slenderness
C = C./alpha;


end





function dy_1 = uplift_DGL(t_1,y_1,time_r,time,ugdd_phi,par)

% Parameters required
dt=time.dt;
g=par.g;
alpha=par.alpha;
R = par.R;
phat = par.phat;
gamma = par.gamma;

% Interpolation
i=floor(t_1/dt);
if i<length(time_r)-1
    a=ugdd_phi(i)+(t_1-time_r(i))/(time_r(i+1)-time_r(i))*...
    (ugdd_phi(i+1)-ugdd_phi(i));
else
    a=0;
end

% Equation of motion; State space formulation
dy_1(1,1) = y_1(2);


% dy_1(2,1)=-12*g/R*(2*gamma+1)*(sin(alpha-y_1(1))-a/g*cos(alpha-y_1(1)))/(15+(cos(alpha))^2+48*gamma);
dy_1(2,1)=-phat^2*(sin(alpha-y_1(1))-a/g*cos(alpha-y_1(1)));
end

function dy = podium_DGL(y,t,time_r,time,ugdd,vgdd,par)

% Parameters required
R=par.R;
alpha=par.alpha;
dt=time.dt;
g=par.g;

gamma=par.gamma;
phat = par.phat;

% Interpolation 
i=floor(t/dt);
if i < length(time_r)
    u = ugdd(i)+(t-time_r(i))/(time_r(i+1)-time_r(i))*(ugdd(i+1)-ugdd(i));
    v = vgdd(i)+(t-time_r(i))/(time_r(i+1)-time_r(i))*(vgdd(i+1)-vgdd(i));
else
    u = 0;
    v = 0;
end

% Equation of motion; State space formulation
dy(1,1)=y(2);



dy(2,1) = -phat^2*(sin(alpha-y(1))+cos(alpha-y(1))*(u/g*cos(y(3))+v/g*sin(y(3))))+...
    y(4)^2*(sin(y(1))*sin(alpha)^2*(4*gamma+3/2)+cos(alpha)*sin(alpha)*(4*gamma+1)*(1+cos(y(1))-2*cos(y(1))^2)+sin(y(1))*cos(y(1))*(4/3*cos(alpha)^2-5/4*sin(alpha)^2+4*gamma*(cos(alpha)^2-sin(alpha)^2)))...
    *12/(48*gamma+15+cos(alpha)^2);


dy(3,1)=y(4);



dy(4,1) = ((2*gamma+1)*(sin(alpha)+sin(y(1)-alpha))*(u*sin(y(3))-v*cos(y(3)))/R...
-((3+8*gamma)*sin(alpha)^2*sin(y(1))+(8/3*cos(alpha)^2-5/2*sin(alpha)^2+8*gamma*(cos(alpha)^2-sin(alpha)^2))*sin(y(1))*cos(y(1))+2*cos(alpha)*sin(alpha)*(1+4*gamma)*(2*sin(y(1))^2+cos(y(1))-1))*y(2)*y(4))/...
((4/3*cos(alpha)^2-5/4*sin(alpha)^2+4*gamma*(cos(alpha)^2-sin(alpha)^2))*sin(y(1))^2+(3+8*gamma)*sin(alpha)^2*(1-cos(y(1)))+2*sin(alpha)*cos(alpha)*sin(y(1))*(1+4*gamma)*(1-cos(y(1))));
end