clc;
clear all
xgn=1024;
s0wn=xgn/2;
global xg 
xg(1,1:xgn)=zeros(1,xgn);
s0w(1,1:s0wn)=zeros(1,s0wn);
fs(1,1:s0wn)=zeros(1,s0wn);
Sa(1,1:s0wn)=zeros(1,s0wn);
dt=0.02;
Td=dt*xgn;
Ts=0.4*Td
t1=0.5*Ts; 
t2=1.2*Ts; 
C=2.5/Ts;
G=9.8
alfa=0.9*G%%8度罕遇地震
kecig=0.05;%阻尼比
Tg=0.25;
gama0=0.9+(0.05-kecig)/(0.5+5*kecig);
yita1=0.02+(0.05-kecig)/8;
yita2=1+(0.05-kecig)/(0.06+1.7*kecig);
D=0.03;
Domiga=2*pi/Td; 

for i=1:s0wn;
 omiga(i)=Domiga*(i); 
end

global Sat
Sat(1,1:s0wn)=zeros(1,s0wn);
for i=1:s0wn
    if 2*pi/omiga(i)<0.1 
        r=(yita2-0.45)*2*pi*alfa/omiga(i)/0.1+0.45*alfa; 
    elseif (2*pi/omiga(i)>=0.1&2*pi/omiga(i)<=Tg) 
        r=alfa*yita2; 
    elseif (2*pi/omiga(i)>=Tg&2*pi/omiga(i)<=5*Tg)  
        r=yita2*alfa*(Tg*omiga(i)/2/pi)^gama0;
    elseif (2*pi/omiga(i)>=5*Tg&2*pi/omiga(i)<=6.0)
        r=(yita2*0.2^gama0-yita1*(2*pi/omiga(i)-5*Tg))*alfa;
     else
        r=0;
    end
    Sat(i)=r;
end

for i=1:s0wn
    ad1=-log(-pi*log(0.999)/omiga(i)/Td);
    s0w(i)=6*kecig.*Sat(i)^2/pi/ad1/omiga(i);
    fzp(i)=2*sqrt(s0w(i)*Domiga);
    xiangweip(i)=rand(1)*2*pi;
    for j=1:xgn
        xg(j)=xg(j)+fzp(i)*cos(omiga(i)*j*dt+xiangweip(i)); 
    end 
end 

v=mean(xg);
for i=1:xgn
    if i*dt<t1 
        fai=(i*dt/t1)^2; 
    elseif (i*dt>=t1&i*dt<t2) 
        fai=1; 
    else i*dt>=t2 
        fai=exp(-C*(i*dt-t2)); 
    end 
    xg(i)=fai*xg(i); 
end 

fushipu=fft(xg);
fzp=abs(fushipu)/s0wn;
xiangweip=angle(fushipu);

for i=1:s0wn
    if (i==1|i==s0wn)
        fs(i)=2*sqrt(s0w(i)*Domiga);
    else 
        fs(i)=sqrt(s0w(i)*Domiga);
    end
   
    if xiangweip(i)<0
        xiangweip(i)=2*pi+xiangweip(i);
    end    
end

for j=1:xgn
        xg(j)=0; 
end 

for i=1:s0wn
    fzp(i)=fs(i);
    for j=1:xgn
        xg(j)=xg(j)+fzp(i)*cos(omiga(i)*j*dt+xiangweip(i)); 
    end 
end 

for cs=1:10
wuca=0.01
Xa=0;
%Xa1=0;
Xa(1,1:s0wn)=zeros(1,s0wn);
Xka(1,1:s0wn)=zeros(1,s0wn);
for j=1:s0wn %   反應譜點
    for i=1:xgn % 時間
        Xa1=0;
        for k=1:s0wn% 第k個分量 
          beita=omiga(k)/omiga(j);
          omigaD=omiga(j)*sqrt(1-kecig^2);
          Cxs=-1/omiga(j)^2/((1-beita^2)^2+(2*kecig*beita)^2);
          G1=2*kecig*beita*Cxs;
          G2=(1-beita^2)*Cxs;
          Bk=-G1*sin(xiangweip(k))-G2*sin(xiangweip(k));
          Ak=(kecig*omiga(j)*Bk-omiga(k)*(G1*cos(xiangweip(k))-G2*sin(xiangweip(k))))/omigaD;
          yy(k)=-omiga(k)^2*(G1*sin(omiga(k)*i*dt+xiangweip(k))+G2*cos(omiga(k)*i*dt+xiangweip(k)))+cos(omiga(k)*i*dt+xiangweip(k));
          zz(k)=(2*kecig^2-1)*omiga(j)^2*(Ak*sin(omigaD*i*dt)+Bk*cos(omigaD*i*dt))*exp(-kecig*omiga(j)*i*dt)-2*kecig*omiga(j)*omigaD*(Ak*cos(omigaD*i*dt)-Bk*sin(omigaD*i*dt))*exp(-kecig*omiga(j)*i*dt)+yy(k);
          Xka1(k)=fzp(k)*zz(k);% 第k個分量在 i*dt 時的結構響應
          Xa1=Xa1+yy(k)*fzp(k)/1.25;%i*dt 時的結構總響應
        end

        if abs(Xa1)>abs(Xa(j))
            Xa(j)=Xa1;
            Xka(j)=Xka1(j);
       end
   end    
   Sa(j)=abs(Xa(j));  
end 
  
   
for j=1:1:s0wn %   反應譜點   
   Ew(j)=abs((Sa(j)-Sat(j))/Sat(j));
   if  Ew(j)>=wuca;
      Rw=Sat(j)/Sa(j); 
        xzdn=1;
              for z=1:xzdn
                     fzp(j)=fzp(j)*Rw^(sign(Xka(j)*Xa(j))); 
              end
   end 
end 
% 幅值譜修正
end

Xa=0;
%Xa1=0;
Xa(1,1:s0wn)=zeros(1,s0wn);
Xka(1,1:s0wn)=zeros(1,s0wn);
for j=1:s0wn %   反應譜點
    for i=1:xgn % 時間
        Xa1=0;
        for k=1:s0wn% 第k個分量 
          beita=omiga(k)/omiga(j);
          omigaD=omiga(j)*sqrt(1-kecig^2);
          Cxs=-1/omiga(j)^2/((1-beita^2)^2+(2*kecig*beita)^2);
          G1=2*kecig*beita*Cxs;
          G2=(1-beita^2)*Cxs;
          Bk=-G1*sin(xiangweip(k))-G2*sin(xiangweip(k));
          Ak=(kecig*omiga(j)*Bk-omiga(k)*(G1*cos(xiangweip(k))-G2*sin(xiangweip(k))))/omigaD;
          yy(k)=-omiga(k)^2*(G1*sin(omiga(k)*i*dt+xiangweip(k))+G2*cos(omiga(k)*i*dt+xiangweip(k)))+cos(omiga(k)*i*dt+xiangweip(k));
          zz(k)=(2*kecig^2-1)*omiga(j)^2*(Ak*sin(omigaD*i*dt)+Bk*cos(omigaD*i*dt))*exp(-kecig*omiga(j)*i*dt)-2*kecig*omiga(j)*omigaD*(Ak*cos(omigaD*i*dt)-Bk*sin(omigaD*i*dt))*exp(-kecig*omiga(j)*i*dt)+yy(k);
          Xka1(k)=fzp(k)*zz(k);% 第k個分量在 i*dt 時的結構響應
          Xa1=Xa1+yy(k)*fzp(k)/1.25;%i*dt 時的結構總響應
        end

        if abs(Xa1)>abs(Xa(j))
            Xa(j)=Xa1;
            Xka(j)=Xka1(j);
       end
   end    
   Sa(j)=abs(Xa(j));  
end 
  

x=omiga;
plot(x,Sa,x,Sat);
y=2.*pi./x
subplot(3,1,1)
loglog(y,Sa,y,Sat),xlabel('周期 秒'),ylabel('加速度 米/秒2'),title('反應譜')

z=0
for i=10:1:s0wn-4
    if z<=abs(Sa(i)-Sat(i))/abs(Sat(i));
        z=abs(Sa(i)-Sat(i))/abs(Sat(i));
        zz=i
    end
end

xg(1,1:xgn)=zeros(1,xgn);
for i=1:s0wn
    for j=1:xgn
        xg(j)=xg(j)+fzp(i)*cos(omiga(i)*j*dt+xiangweip(i)); 
    end 
end 
t=0.02:0.02:Td;
plot(t,xg),xlabel('時間 秒'),ylabel('加速度 米/秒2'),title('加速度時程')

save 4.mat