EOM para edificios de N-Pisos / Idealizado como edificio de cortante
clear variables
clc

% Especificamos un array con el conjunto de masas del edificio
m_array=[30000; % Piso 1
         30000; % Piso 2
         30000; % Piso 3
         30000; % Piso 4
         25000; % Piso 5
         25000; % Piso 6
         ];

% Especificamos un array con el conjunto de rigideces laterales del
% edicicio
k_array=[2500000
         2500000
         2500000
         2100000
         2100000
         2100000
         ];

% Especificamos un array con el conjunto de amortiguamientos modales
xi_array=[0.05
          0.05
          0.04
          0.01
          0.01
          0.01
          ];

% Definimos la altura de entrepiso
H_entrepiso=6; %[m]

Definimos un registro sismico

filename='NGA181_ImperialValley06_ElCentroArray6_140.txt';
EQ=Leer_Registro_NGA(filename,4);

Calculos de parametros analiticos

nDof=length(m_array);
% Definimos un vector que contenga todas las alturas de los pisos
Altura_piso=0:H_entrepiso:H_entrepiso*length(m_array);
% Definimos el tiempo que queremos correr la simulacion
Tiempo_max=100; %[s]

% Realizamos el calculo para extender el regitro al tiempo maximo
% establecido

if EQ.t(end,1)<Tiempo_max
    EQ_calc=EQ;
    t=(EQ.t(end,1):EQ.dt:Tiempo_max)';
    z=zeros(length(t),1);
    EQ_calc.A=[EQ.A;z];
    EQ_calc.V=[EQ.V;z];
    EQ_calc.D=[EQ.D;z];
    EQ_calc.t=[EQ.t;t];
    clear z t
else
    EQ_calc=EQ
end

Contruimos las matrices de masas y rigidez
% Ensamblamos la matriz de masas
M=zeros(length(m_array));
for i=1:length(m_array)
    M(i,i)=m_array(i,1);
end

% Ensamblamos la matriz de rigidez (Como edificio de cortante)
K=zeros(length(k_array));
for i=1:length(k_array)
    if i>1
        k=k_array(i,1);
        k_temp=[k,-k;-k,k];
        K(i-1:i,i-1:i)=K(i-1:i,i-1:i)+k_temp;
    else
        K(i,i)=K(i,i)+k_array(i,1);
    end
end

Analisis modal
[fi_hat,omega2]=eig(K,M);
wn=sqrt(omega2);
Tn=2*pi./wn

% Normalizamos la matriz modal para dibujo

fi_draw=zeros(size(fi_hat));
for i=1:length(fi_hat)
    fi_draw(:,i)=fi_hat(:,i)/fi_hat(1,i);
end

% Normalizamos la matriz modal para la masa

fi=zeros(size(fi_hat));

for i=1:length(fi)
    fi(:,i)=fi_hat(:,i)/(fi_hat(:,1)'*M*fi_hat(:,1))^0.5;
end

% Verificamos la normalizacion de la matriz modal
M_hat=fi'*M*fi; % Esta prodia ser definida directamente como una identidad
K_hat=fi'*K*fi; % Esta prodia ser definida directamente como omega2
omega2

fi_draw=[zeros(1,length(fi_draw));fi_draw];

% Dibujos de los modos
len=length(fi_draw)-1;
min_x=ceil(abs(min(min(fi_draw))));
max_x=ceil(max(max(fi_draw)));
max_px=max(min_x,max_x);
hor=ceil((len-1)/2);
ver=len-1-hor;

% Ploteamos los modos de vibracion
for i=1:len
    subplot(ver,hor,i)
    plot(fi_draw(:,i),Altura_piso,'.k','MarkerSize',25)
    hold on
    plot(fi_draw(:,i),Altura_piso,'k')
    grid on
    pbaspect([1 1 1])
    title(['Modo ',num2str(i),' - Tn=',num2str(round(Tn(i,i),2)),'[s]'])
    axis([-max_px-1 max_px+1 0 Altura_piso(1,end)])
    xlabel('Coordenada X')
    ylabel('Coordenada Y')
    set(gcf,'Visible','on','position',[0 0 1920 1080])
    hold on
end

clear ver hor len min_x max_x max_px

Respuesta al registro usando superposicion modal
Peff=zeros(length(m_array),length(EQ_calc.A));
for i=1:length(m_array)
    Peff(i,:)=-M(i,i)*EQ_calc.A*9.81;
end

% Convertimos la fuerza efectiva en coordenadas modales
P_hat=fi'*Peff;

% Definimos el metodo numerico para calcular la respuesta de cada modo
% [xt, vt, at, coef, Tabla, Check] = Newmark(m, k, xi, t, P, delta_t, x0, v0, beta, gamma)
beta=1/4;
gamma=1/2;
x0=0;
v0=0;

% Calculamos la respuesta modal del sistema
vt=zeros(length(fi),length(P_hat));
for i=1:length(fi)
    [vt(i,:),~,~,~,~,~]=Newmark(M_hat(i,i), K_hat(i,i), xi_array(i,1), EQ_calc.t, P_hat(i,:), EQ_calc.dt, x0, v0, beta, gamma);
end

% Ploteamos la respuesta modal
figure
hold on
for i=1:length(fi)
    txt=['Modo #',num2str(i)];
    plot(EQ_calc.t,vt(i,:),'DisplayName',txt)
    grid on
    xlabel('Tiempo [s]')
    ylabel('Desplazamiento [m]')
    title('Respuesta modal')
    set(gcf,'Visible','on','position',[0 0 1920 1080])
end
hold off
legend show

% Calculamos la respuesta en coordenadas de estructura

xt=fi*vt;

% Ploteamos la respuesta estructural
figure
    hold on
for i=1:height(xt)
    txt=['Piso #',num2str(i)];
    plot(EQ_calc.t,xt(i,:),'DisplayName',txt)
    grid on
    xlabel('Tiempo [s]')
    ylabel('Desplazamiento [m]')
    title('Respuesta Estructural')
    set(gcf,'Visible','on','position',[0 0 1920 1080])
end
hold off
legend show

% Ploteamos la respuesta de cada masa de forma independiente

max_draw=max(max(xt));

for i=1:height(xt)
    subplot(height(xt),1,i)

    [Peak,idx]=max(xt(i,:));

    txt=['Piso #',num2str(i)];
    plot(EQ_calc.t,xt(i,:))
    hold on
    plot(EQ_calc.t(idx),xt(i,idx),'or')
    text(EQ_calc.t(idx),xt(i,idx),['( ' num2str(round(EQ_calc.t(idx),2)) ' , ' num2str(round(xt(i,idx),2)) ' )'])
    axis([EQ_calc.t(1,1) EQ_calc.t(end,1) -1.25*max_draw 1.25*max_draw])
    grid on
    xlabel('Tiempo [s]')
    ylabel('Desplazamiento [m]')
    title(txt)
    set(gcf,'Visible','on','position',[0 0 1920 1080])
end