铣削稳定性叶瓣图 MATLAB代码

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