%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Programa de analisis modal de edificios %
%Calcula auto valores y auto vectores    %
%Para la vibracion libre de la estructura%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Funcion aceleracion calcula periodos y aceleracones de diseño
%Ademas calcula las fuerzas de piso por cada modo de vibracion.

%-------------------------------------------------------------
function[F,omega,V,UI]=aceleracion(npis,M,FI,OM,hpis)
for i=1:npis
    T(i)=2.*pi/(sqrt(OM(1,i)));
end
sprintf('Periodos')
T'
pause

%calculo del factor de participacion modal
for i=1:npis
    sum2=0;
    sum22=0;
    num1=0;
    num11=0;
    for j=1:npis
        num1=M(j,j)*FI(j,i);
        num11=M(j,j)*FI(j,i)*FI(j,i);
        sum2=sum2+num1;
        sum22=sum22+num11;
        
    end
    FPM(i)=sum2/sum22;
    
end
sprintf('factor de participacion modal')
FPM
pause
%factor de distribucion sismica
sprintf('factor de distribucion sismica')
FDS=FI*FPM'



%---------------------------------------------------------------
%determinacion de las aceleraciones de diseño
KC=input('introduzca el coeficiente de contribucion')
TS=input('indique el tipo de suelo')
if TS==1.0
    CT=1.0;
end
if TS==2
    CT=1.3;
end
if TS==3
    CT=1.6;
end
if TS==4
    CT=2.0;
end
%calculo de los periodos caracteristicos del espectro
%norma española NCSE-2002
TA=KC*CT/10;
TB=KC*CT/2.5;

%calculo de las aceleraciones de diseño de cada modo
for i=1:npis
    if T(i)<TB
        if T(i)>TA
            AD(i)=2.5;
        elseif T(i)<TA
            AD(i)=1.0+(1.5*T(i)/TA);
        end
    end
    
    
    if T(i)>TB
        AD(i)=KC*CT/T(i);  
    end
    
end
sprintf('Aceleraiciones de diseño normalizadas')
AD


%determinacion de la aceleracion sismica de calculo
AB=input('introduzca la aceleracion basica, conforme al mapa de zonificacion')
RO=input('introduzca el valor del factor de importancia')
g=9.81;
FC=AB*RO*g;

%se determina el coeficiente de amplificacion del terreno S
for i=1:npis
    if FC<0.1*g
        ESE=CT/1.25;
    end
    if FC>0.1*g|FC<0.4*g
        ESE=(CT/1.25)+3.33*((RO*AB)-0.1)*(1-CT/1.25);
    end
    if FC>0.4*g
        ESE=1.0;
    end
end
%aceleracion de calculo ac
ESE
pause
AC=ESE*RO*AB*g;

%----------------------------------------------------------------
%aceleracion de diseño elastico de cada modo
for i=1:npis
A(i)=AC*AD(i);
end
sprintf('Aceleraciones de diseño')
A
%aceleracion de diseño inalastico
omega=input('introduzca la fraccion de amortiguamiento en %')
miu=input('introduzca el factor de ductilidad')
vi=(5/omega)^.4;
coer=vi/miu;
for i=1:npis
    AI(i)=A(i)*coer;
end
sprintf('Aceleraiciones inelasticas de diseño')
AI=AI

%----------------------------------------------------------------

%determinacion de las aceleraciones modales de piso
for i=1:npis
    for j=1:npis
        AMP(j,i)=FI(j,i)*FPM(i)*AI(i);
    end
end
AMP
%determinacion de los desplazamientos modales de piso (elasticos)
for i=1:npis
    for j=1:npis
        UE(j,i)=AMP(j,i)/(OM(i));
    end
end 
UE';
%clc
%determinacion de los desplazamientos modales de piso (inelasticos)
sprintf('Desplazamientos modales inelasticos en mm')
UI=1000*miu*UE
pause

%----------------------------------------------------------------

%----------------------------------------------------------------
%Verificacion de las derivas de piso (segun norma venezolana)
UI=UI*.8;
drift=zeros(npis,npis);
delta=zeros(npis,npis);
for i=1:npis
    delta(1,i)=(1/1000)*UI(1,i);
    for j=2:npis
        delta(j,i)=(1/1000)*(UI(j,i)-UI(j-1,i));
    end
end

for i=1:npis
    for j=1:npis
        drift(j,i)=delta(j,i)/hpis(j);
    end
end

for i=1:npis
    for j=1:npis
        if drift(j,i)>0.015
            sprintf ('la deriva del piso %d supera el valor normativo',j)
            sprintf ('modo %d ',i)
            sprintf ('REDIMENSIONAR')
        end
    end
end
sprintf('Derivas de piso')
drift
sprintf('Derivas de piso en %')
deriva=100*drift

pause


%----------------------------------------------------------------

%----------------------------------------------------------------

%determinacion de las fuerzas modales de piso
for i=1:npis
    for j=1:npis
        F(j,i)=M(j,j)*FI(j,i)*FPM(i)*AI(i);
    end
end
sprintf('Fuerzas de piso, por modos, en Newtons')
F=F
pause

%-----------------------------------------------------------------

%determinacion de las fuerzas cortantes de piso
for i=1:npis
    sum=0;
    for j=npis:-1:1
        sum=sum+F(j,i);
        V(j,i)=sum;
    end
end
sprintf('Fuerzas cortantes de piso, por modos, en Newtons')
V
pause