% spec(per,zai,dt,p,nsteps)
% ******************
% Calculation of elastic response spectra 
% The program calculates Sd and Sv as well as the pseudo
% velocity and acceleration spectra PSV and PSA
% The integration is carried out using the Newmark-Beta algorithm
% at this time CAA is used.
% Zero initial conditions are assumed and for consistency the first
% value on the acceleration record should be zero. 
% ***********
% Input Data
% ***************************************************************
%   per = vector with periods at which the spectral ordinates are
%         computed.
%   zai = pecent of critical damping
%  dt = integration time step
%   p = earthquake time history
% ***************************************************************
% function [psa,psv,sd,sv]=spec(per,zai,dt,p,nsteps)

function [psa,psv,sd,sv,sa]=spec(per,zai,dt,p,nsteps)

beta=0.25;
gama=0.5;
nper=length(per);
for ii=1:nper;
omega=2*pi/per(ii);
c=2*omega*zai;
k=omega^2;

gb=gama/beta;
kd=1/(beta*dt^2)+gb*c/dt+k;
kdi=1/kd;
nc=length(p);
if nc<nsteps+1;
disp('el vector "p" debe tener al menos nsteps+1 puntos')

break;
else;
end;
  d(1)=0;
  v(1)=0;
  a(1)=0;
  t(1)=0;
for i=1:1:nsteps;
     j=i+1;
    dp=p(j)-p(i);
    vec1=(0.5*a(i)+1/dt*v(i))/beta;
    vec2=(gb/2-1)*dt*a(i)+gb*v(i);    
    dld=dp+vec1+c*vec2;
    dd=kdi*dld;
    dv=(1-0.5*gb)*dt*a(i)-gb*v(i)+gb/dt*dd;
    d(j)=d(i)+dd;
    v(j)=v(i)+dv;
    a(j)=(p(j)-c*v(j)-k*d(j));
    at(j)=(a(j)+p(j));
    t(j)=t(i)+dt;
end;

sd(ii)=max(abs(d));
sv(ii)=max(abs(v));
sa(ii)=max(abs(at));
psv(ii)=omega*sd(ii);
psa(ii)=(omega^2)*sd(ii);
end;
