tic
%！傳遞矩陣法！
clc;clear;

%參數輸入：
n?=?20;%建筑總層數
h?=?3.6;%H/n;%層高（m）
H?=?n*h;%建筑總高（m）
m?=?23437*ones（n1）;%各層沿高度的分布質量（kg/m）
M?=?84375*ones（n1）;%各層集中質量（kg）
EI?=?2.6447*10^12*ones（n1）;%各層抗彎剛度（N/m2）
GA?=?4.0387*10^10*ones（n1）;%各層抗剪剛度（N/m2）
v?=?0.02;%阻尼
g?=?9.8;%重力加速度（m/s2）
k_?=?1.2;%剪應力不均勻系數


%********************************************地震反應計算*******************************************************
%參數輸入：
EI_?=?EI*（1+v*1i）;
GA_?=?GA*（1+v*1i）;
%傅里葉變換
???load?a3
???T=max（a（:1））;
???a（:2）?=?a（:2）/1000;%統一單位m
???Ag=fft（a（:2））;
???N=length（a）;
??%計算omega
???omega=zeros（N1）;
???for?j=1:N
?????omega（j）=2*pi/T*j;
???end
???%對Ag進行處理
??for?j=1:N????????????????
??????if?j>=2&&j<=N/2
??????????Ag（j）=2*Ag（j）;
??????end
??????if?j>=N/2+2&&j<=N
??????????Ag（j）=0;
??????end
??end
??
%參數計算:
alpha?=?k_*h*h*（omega.^2*（m./GA_）‘）;%alpha（:n）表示第n層的所有omega對應的值，下同
beta?=?h^4*（omega.^2*（m./EI_）‘）;%beta（N*n）
lamda1?=?（（alpha.^2/4+beta）.^0.5?-?alpha/2）.^0.5;%lamda1（N*n）
lamda2?=?（（alpha.^2/4+beta）.^0.5?+?alpha/2）.^0.5;%lamda2（N*n）

d?=?ones（N1）*（h*h*EI_.^-1）‘;%d（N*n）
deta?=?（lamda1.^2?+?lamda2.^2）.^-1;%deta（N*n）
c0?=?deta.*（lamda2.^2.*cosh（lamda1）?+?lamda1.^2.*cos（lamda2））;
c1?=?deta.*（lamda2.^2./lamda1.*sinh（lamda1）?+?lamda1.^2./lamda2.*sin（lamda2））;
c2?=?deta.*（cosh（lamda1）?-?cos（lamda2））;
c3?=?deta.*（sinh（lamda1）./lamda1?-?sin（lamda2）./lamda2）;

%場傳遞矩陣tf（4*4*N*n）
tf（11::）?=?c0?-?alpha.*c2;
tf（12::）?=?-beta.*c3/h;
tf（13::）?=?-beta.*c2./d;
tf（14::）?=?beta.*（alpha.*c3?-?c1）./d/h;
tf（21::）?=?-h*（c1?-?alpha.*c3）;
tf（22::）?=?c0;
tf（23::）?=?h*beta.*c3./d;
tf（24::）?=?beta.*c2./d;
tf（31::）?=?-d.*c2;
tf（32::）?=?d.*c1/h;
tf（33::）?=?c0;
tf（34::）?=?beta.*c3/h;
tf（41::）?=?h*d.*（alpha.*c1?-?（beta+alpha.^2）.*c3）./beta;
tf（42::）?=?d.*c2;
tf（43::）?=?h*（c1?-?alpha.*c3）;
tf（44::）?=?c0?-?alpha.*c2;

%點傳遞矩陣tp（4*4*N*n）
tp（11::）?=?ones（Nn）;
tp（22::）?=?ones（Nn）;
tp（33::）?=?ones（Nn）;
tp（44::）?=?ones（Nn）;
tp（12::）?=?zeros（Nn）;?
tp（12::）?=?zeros（Nn）;???
tp（13::）?=?zeros（Nn）;???
tp（14::）?=?-（omega.^2）*M‘?;
tp（21::）?=?zeros（Nn）;???
tp（23::）?=?zeros（Nn）;
tp（24::）?=?zeros（Nn）;???
tp（31::）?=?zeros（Nn）;???
tp（32::）?=?zeros（Nn）;???
tp（34::）?=?zeros（Nn）;?
tp（41::）?=?zeros（Nn）;???
tp（42::）?=?zeros（Nn）;???
tp（43::）?=?zeros（Nn）;

%求總傳遞矩陣Tn
Tn_?=?1;
for?i=1:N
????for?j=1:n
????Tn_?=?tp（::ij）*tf（::ij）*Tn_;
????%Tn（::ij）?=?Tn_;?
????end
????TN（::i）?=?Tn_;
????Tn_?=?1;
end

%求地震反應kesi
Ug?=?-Ag./（omega.^2）;%傅里葉變換后的地面位移
temp?=?（TN（24:）.*TN（12:）?-?TN（14:）.*TN（22:））./（TN（22:）.*TN（11:）?-?TN（12:）.*TN（21:））;
Q_I（:1）?=?temp（:）.*Ug（:）;
%第一層層底剪力
temp?=?（TN（14:）.*TN（21:）?-?TN（24:）.*TN（11:））./（TN（22:）.*TN（11:）?-?TN（12:）.*TN（21:））;
M_I（:1）?=?temp（:）.*Ug（:）;
%第一層層底彎矩


kesi_I（1:1）?=?Q_I（:1）;?
kesi_I（2:1）?=?M_I（:1）;?
kesi_I（3:1）?=?zeros（N1）;
k