clear?all;
clc;
%控制參數，用國際單位制
f=20000;?????????????????????????????%頻率控制，單位為rpm
E=2.01e+011;?????????????????????????%軸段材料的彈性模量單位為N/m^2
P=7850;??????????????????????????????%材料密度Kg/m^3
%安裝軸承部位將軸同樣分成兩段，盤分成兩段2個軸承，2個盤
l=[0.018;0.01;0.008;0.008;0.1;0.012;0.019;0.019;0.432;0.019;0.019;0.012;0.03;0.006;0.006;0.01;0.068;0.009];%分成18段，各軸段長度
d=[0.024;0.025;0.03;0.03;0.033;0.037;0.0419;0.0419;0.033;0.0419;0.0419;0.037;0.033;0.030;0.030;0.025;0.0235;0.018];%各軸段直徑
num=[18;10;8;8;100;12;19;19;432;19;19;12;30;6;6;10;68;9];%各軸段分段數

m_shaft_total=0;
s=0;
m_shaft_L1=0;
for?i=1:18
????m_shaft（i）=pi*d（i）^2/4*l（i）*P;????????%各軸段質量
????m_shaft_total=m_shaft_total+m_shaft（i）;%軸總質量
????if?s>3?&&?s<15
????m_shaft_L1=m_shaft_L1+m_shaft（s）;??????%軸承中間的軸段質量
????end
????J（i）=pi*（d（i）^4）/64;%軸段截面慣性矩
????s=s+1;
end

Dl=0.26;????????????????????????????????????????????????%盤L徑
Bl=0.04;???????????????????????????????????????????????%盤L厚
m_disk_L=pi*Dl^2/4*Bl*P-pi*d（7）^2/4*Bl*P;???????????????????%盤L質量
Ip_L=1/2*m_disk_L*（Dl^2/4+d（7）^2/4）;??????????????????????????%盤L極轉動慣量
Id_L=1/12*m_disk_L*（3*（Dl^2/4+d（7）^2/4）+Bl^2）;?????????????????%盤L直徑轉動慣量

Dr=0.28;????????????????????????????????????????????????%盤R徑
Br=0.04;???????????????????????????????????????????????%盤R厚
m_disk_R=pi*Dr^2/4*Br*P-pi*d（14）^2/4*Br*P;???????????????????%盤R質量
Ip_R=1/2*m_disk_R*（Dr^2/4+d（14）^2/4）;??????????????????????????%盤R極轉動慣量
Id_R=1/12*m_disk_R*（3*（Dr^2/4+d（14）^2/4）+Br^2）;?????????????????%盤R直徑轉動慣量

M_total=m_shaft_total+m_disk_L+m_disk_R;??????????????????%實驗器總重量

bearing_k1=3e+006;
bearing_k2=1e+007;??????????????????%支承剛度，根據實際情況給定bearing_k1為滾棒軸承支承剛度，bearing_k2為3#軸承支承剛度

P_matrix=[1000;
??????????0100;
??????????0010;
??????????0001];?????????????????????%初始矩陣

B_matrix1=[1000;
??????????0100;
??????????0010;
?????????-bearing_k1001];?
?????
B_matrix2=[1000;
??????????0100;
??????????0010;
?????????-bearing_k2001];

for?i=1:1:f
????Omega=i/60*2*pi;???%角速度
????%第一大段，細分成18小段
????n=num（1）;
????N=1;
????L=l（N）/n;
?????M_matrix=[1??????000?;
???????????????0??????100?;
???????????????0??????010?;
???????????m_shaft（N）/n*Omega^2?001];?
????L_matrix=?[1?????????L????????（L^2）/（2*E*J（N））（L^3）/（6*E*J（N））;
???????????????0?????????1???????????L/（E*J（N））??????（L^2）/（2*E*J（N））;
???????????????0?????????0???????????1??????????????????L???????????;
???????????????0?????????0???????????0??????????????????1??????????];
????H_matrix=（M_matrix*L_matrix）^n*P_matrix;

????%第二大段，細分成10小段
????n=num（2）;
????N=2;
?????L=l（N）/n;
?????M_matrix=[1??????000?;
???????????????0??????100?;
???????????????0??????010?;
???????????m_shaft（N）/n*Omega^2?001];?
????L_matrix=?[1?????????L????????（L^2）/（2*E*J（N））（L^3）/（6*E*J（N））;
???????????????0?????????1???????????L/（E*J（N））??????（L^2）/（2*E*J（N））;
???????????????0?????????0???????????1??????????????????L???????????;
????????????