clear?all
close?all
clc
tic????????????%計時開始，檢驗程序性能

%?定義x方向參數
kx?=?2.1425e6;???????????????%?N/mx方向剛度
wnx?=?4035*2*pi;?????????%?rad/sw=2*pi*f??轉速大?頻率大?阻尼比（模態參數）
zetax?=?0.01;????????????%x方向阻尼比

%?定義y方向參數?
ky?=?2.1425e6;???????????????%?N/m??y方向剛度
wny?=4035*2*pi;?????????%?rad/s???y方向阻尼比
zetay?=?0.01;??????????%y方向阻尼比

phis?=?0*pi/180;???????????????%?rad??切入角?（逆銑）
phie?=?180*pi/180;??????%?rad??切出角??（逆銑）

%?切削力系數?
%?Ks?=?750e6;?????????????%?N/m^2
%?beta?=?68;??????????????%?deg??螺旋角（以圓柱螺旋立銑刀為背景的）
Kn?=?1550/2550;?
Kt?=?1.482e9;?%?N/m^2??切向切削力系數?
kr=0.889e9;??????%徑向切削力系數?
Nt?=?2;??%齒數
alphaxx?=?0.5*（（cos（2*phie）-2*Kn*phie+Kn*sin（2*phie））-（cos（2*phis）-2*Kn*phis+Kn*sin（2*phis）））;
alphaxy?=?0.5*（（-sin（2*phie）-2*phie+Kn*cos（2*phie））-（-sin（2*phis）-2*phis+Kn*cos（2*p