% program about base isolation structure
% 行末注有@表明這行參數在不同的算例中有可能會變化
XZ=0;                                                % 陡線中點橫坐標值
XS=0.02;                                             % 隔震層屈服位移                @
cs1=0.05;                                            % 上部結構第一振型阻尼比        @
cs2=0.05;                                            % 上部結構第二振型阻尼比        @
k1=1.5;                                              % 隔震層彈性剛度                @
k2=0.6;                                              % 隔震層塑性剛度                @
ug=load('d:\data\input\wave.txt','-ascii');          % 輸入地震波
ug=ug';                                              % 轉置
ug=ug(:);                                            % 轉為單下標標識
n1=length(ug);                                       % 輸入地震波記錄點數
n1=2000;                                              % 計算點數
ugm=341.7/max(abs(ug));                              % 調整加速度峰值
ug=0.01*ug*ugm;                                      % 單位轉化（1厘米/秒*秒＝0.01米/秒*秒）
m0=[41 53.2 61.3 54.5 54.5 54.5 32.4]*1e+3;          % 輸入質量陣 %%%@  黃建文
n=length(m0);                                        % 質點數
M=diag(m0);                                          % 化為對角陣
P=zeros(n,1);                                        % 為各層的恢復力列向量(方程右邊項）
p(1)=0;                                              % 為隔震層剛度判別因子
k0=[k1 115.3 115.3 74.9 74.9 74.9 74.6]*1e+6;        % 輸入剛度陣                   @
[K]=matrix(k0,n);                                    % 調用子程序求剛度陣
[x,d]=eig(K,M);                                      % 求結構的自振特性
d=sqrt(d);%
w=sort(diag(d));                                     % 結構的自振頻率
T=2*3.1415926./w;                                    % 結構的自振周期
c1=2*w(1)*w(2)*(cs1*w(2)-cs2*w(1))/(w(2)^2-w(1)^2);  % 上部結構瑞雷阻尼系數
c2=2*(cs2*w(2)-cs1*w(1))/(w(2)^2-w(1)^2);            % 上部結構瑞雷阻尼系數
%C=c1*M+c2*K;                                        % 阻尼陣
c0=[43.5 171.5 171.5 111.4 111.4 111.4 111.4]*1e+3;  % 輸入各層阻尼值
[C]=matrix(c0,n)                                     % 調用子程序求阻尼陣
dt=0.005;                                             % 時間步長                      @
ug0=[ug(1)*ones(n,1)];%
F0=-M*ug0;                                           % 計算初始力
X0=zeros(n,1);                                       % 計算初始位移
X1=zeros(n,1);                                       % 計算初始速度
X2=inv(M)*(F0-C*X1-K*X0);                            % 計算初始加速度
XX0(1:n,1)=X0;%
XX1(1:n,1)=X1;%
XX2(1:n,1)=X2;%

%計算wilson法的幾個系數
sita=1.37;%
a0=6/(sita*dt)^2;%
a1=3/(sita*dt);%
a2=2*a1;                  
a3=sita*dt/2;%
a4=a0/sita;%
a5=-a2/sita;%
a6=1-3/sita;%
a7=dt/2;%
a8=dt*dt/6;%
Ki=K+a0*M+a1*C;                                      % 形成有效剛度矩陣



% 對每一時間步長作如下計算
i=2;
while(i<n1+1);
    ug0=[ug(i)*ones(n,1)];                           % 地震加速度
    F_t=-M*ug0-P;%
    Ft0=F0+sita*(F_t-F0)+M*(a0*X0+a2*X1+2*X2)+C*(a1*X0+2*X1+a3*X2);%
    Xt0=inv(Ki)*Ft0;                                 % (t+sita*dt)的位移
    
    
    % 計算t+dt的加速度，速度，位移
    X_2=X2+(a0*Xt0-a0*X0-2*a1*X1-3*X2)/sita;          %計算t+dt的加速度
    X_1=X1+a7*(X_2+X2);                               %速度
    X_0=X0+dt*X1+a8*(X_2+2*X2);                       %位移
    X_2=-inv(M)*C*X_1-inv(M)*K*X_0-ug(i)-inv(M)*P;    %直接由動力方程計算加速度，消除累積誤差
    
    

       %%%%%%%%%% 求出剛度變化判斷因子p %%%%%%%%%%
    if (X_0(1)-XZ)>=XS&X_1(1)>0;%
        p(i)=1;                                       % 上緩線
        P(1)=(k1-k2)*(1e+6)*XS;                       % 隔震層的恢復力列向量(方程右邊項）    @
        Q(i)=k2*(1e+6)*X_0(1)+P(1);                   % 隔震層的恢復力列向量                 @
    elseif (X_0(1)-XZ)<=-XS&X_1(1)<0%
        p(i)=-1;                                      % 下緩線
        P(1)=(k2-k1)*(1e+6)*XS;                       % 隔震層的恢復力列向量(方程右邊項）    @
        Q(i)=k2*(1e+6)*X_0(1)+P(1) ;                  % 隔震層的恢復力列向量                 @
    else p(i)=0;                                      % 陡線
        P(1)=(k2-k1)*(1e+6)*XZ;%                      % 隔震層的恢復力列向量(方程右邊項）    @
        Q(i)=k1*(1e+6)*X_0(1)+P(1);                   % 隔震層的恢復力列向量                 @
      end
   
      
    
     %%%%%%%%%%   隔震層第二剛度賦值   %%%%%%%%%%
    if (p(i)==1)|(p(i)==-1);                                %
        XZ=X_0(1)-p(i)*XS;                                  %
        k0=[k2 115.3 115.3 74.9 74.9 74.9 74.6]*1e+6;      % 輸入剛度陣            @
        [K]=matrix(k0,n);
        [x,d]=eig(K,M);%
        d=sqrt(d);%
        w=sort(diag(d));
        c1=2*w(1)*w(2)*(cs1*w(2)-cs2*w(1))/(w(2)^2-w(1)^2);  % 上部結構瑞雷阻尼系數
        c2=2*(cs2*w(2)-cs1*w(1))/(w(2)^2-w(1)^2);            % 上部結構瑞雷阻尼系數
        %C=c1*M+c2*K;                                        % 阻尼陣
        Ki=K+a0*M+a1*C;                                      % 形成有效剛度矩陣
    end
    
     if (p(i)==0);                                  % 
        k0=[k1 115.3 115.3 74.9 74.9 74.9 74.6]*1e+6;         % 輸入剛度陣            @
        [K]=matrix(k0,n);
        [x,d]=eig(K,M);                                        %
        d=sqrt(d);                                             %
        c1=2*w(1)*w(2)*(cs1*w(2)-cs2*w(1))/(w(2)^2-w(1)^2);    % 上部結構瑞雷阻尼系數
        c2=2*(cs2*w(2)-cs1*w(1))/(w(2)^2-w(1)^2);              % 上部結構瑞雷阻尼系數
        %C=c1*M+c2*K;                                          % 阻尼陣
        Ki=K+a0*M+a1*C;                                        % 形成有效剛度矩陣
    end;
 %變量傳遞
    X0=X_0;%
    X1=X_1;%
    X2=X_2;%
    F0=-M*ug0-P;%
    XX0(1:n,i)=X0;%
    XX1(1:n,i)=X1;%
    XX2(1:n,i)=X2;%
    XX3(1:n,i)=X2+ug(i);%
    i=i+1;%
end