clc
clear all
disp('BABAK.HM,,,,,,,,,,,Please,,,,,,,,,wait')
c=0;
v=0:0.05:40;
s=length(v);
rsd=zeros(10,s);
rsv=zeros(10,s);
rsa=zeros(10,s);
rsd1=[];
rsv1=[];
rsa1=[];
zetaa=[0,0.02,0.05,0.1];
for q=1;
    c=c+1;
    zeta=zetaa(1,c);
    if q==1:4;
        disp('zeta=0')
    end
    if q==2;
        disp('zeta=0.02')
    end
        if q==3;
            disp('zeta=0.05')
        end
            if q==4;
                disp('zeta=0.1')   
            end
delta=0.005;
for j=2:2:20;
a=importdata('d:\gather.txt');
a=a(:,j);
a(isnan(a))=[];
g=981;%cm/s^2%
b=max(abs(a));
a=(a/b)*g*0.35;%Hampaye kardan%
m=length(a);
y(1:m)=0;
yd(1:m)=0;
ydd(1:m)=0;
yddabs(1:m)=0;
b=0;
for v=0.05:0.05:40.05;
    wn=2*pi*v;
    b=b+1;
    for i=1:(m-1);
        wd=wn*sqrt(1-zeta^2);
    %pish mohasebeh%
    F=exp(-zeta*wn*delta);
    G=zeta/(sqrt(1-(zeta^2)));
    H=sin(wd*delta);
    I=cos(wd*delta);
    K=(1-(2*zeta^2))/(wd*delta);
    J=(2*zeta)/(wn*delta);
    L=wn/(sqrt(1-(zeta^2)));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    A=F*((G*H)+I);
    B=(F*H)/wd;
    C=(1/wn^2)*(J+(F*((K-G)*H-(1+J)*I)));
    D=(1/wn^2)*(1-J+F*((-K*H)+(J*I)));
    AP=-F*L*H;
    BP=F*(I-(G*H));
    CP=(1/wn^2)*((-1/delta)+F*(((L+(G/delta))*H)+(I/delta)));
    DP=(1/(delta*wn^2))*(1-F*((G*H)+I));
    y(i+1)=(A*y(i))+(B*yd(i))+(C*a(i))+(D*a(i+1));
    yd(i+1)=(AP*y(i))+(BP*yd(i))+(CP*a(i))+(DP*a(i+1));
    ydd(i+1)=a(i+1)-(2*zeta*wn*yd(i+1))-(wn^2*y(i+1));
    yddabs(i+1)=(2*zeta*wn*yd(i+1))+(wn^2*y(i+1));
    end
    rsd(j/2,b)=max(abs(y(1,1:m)));
    rsv(j/2,b)=max(abs(yd(1,1:m)));
    rsa(j/2,b)=((max(abs(yddabs(1,1:m))))/981);
    
end
end
rsd1(c,:)=sum(rsd)/10;
rsv1(c,:)=sum(rsv)/10;
rsa1(c,:)=sum(rsa)/10;
end
v=0.05:0.05:40.05;
t=zeros(1,length(v));
for p=1:length(v);
    t(p)=1/v(p);
end
figure(1)
plot(t,rsd1)
figure(2)
plot(t,rsv1)
figure(3)
plot(t,rsa1)
v=0.05:0.05:40.05;
w=2*pi*v;
for i=1:4
prsv(i,:)=w.*rsd1(i,:);
prsa(i,:)=w.*prsv(i,:);
prsa(i,:)=prsa(i,:)./981;
end
figure(4)
semilogx(t,prsv)
figure(5)
semilogx(v,prsv)
figure(6)
semilogx(v,prsa)