摆线针轮减速机摆线方程(Matlab)

最近在研究摆线减速机，手算方程太麻烦，于是抄到了matlab里面

%以下参数手动输入
Rp=30;  %针轮中心圆半径
t=35   ;  %传动比
K1=0.6;  %变幅系数（初取，依照短幅系数表）
delta=0.1 ； %偏置量，保证实体配合     
rrp=2+delta  ;  %针齿半径  需要小于最小曲率且针径系数满足表格要求（针径系数表）
xr=0.15;  %摆线轮齿宽系数 取0.1~0.2

%以下参数自动生成
Zc=t  ;  %摆线轮齿数
Zp=Zc+1    %针轮齿数
a=round(K1*Rp/Zp,1) %中心距
K1=a*Zp/Rp %变幅系数（反求）

rc=a*Zc                                %摆线轮节圆半径
rp=a*Zp                                %针轮节圆半径
rg=a/K1                                %滚圆半径
rbc=rc/K1                              %基圆半径
Pn=2*pi*a/K1                           %啮合齿距
Pmin=Rp*sqrt((27*(1-K1^2)*Zc)/((Zc+2)^3))     %最小曲率半径
Pc=2*pi*(rp-rc)                        %节圆齿距
K2=Rp/rrp*sin(180/Zp/57.3)                 %针径系数

rac=Rp+a-rrp                           %摆线轮齿顶圆半径
rfc=Rp-a-rrp                           %摆线轮齿根圆半径
bc =xr*Rp                              %摆线轮宽度

x=[];
y=[];
for t=0:0.1:135
     x(end+1)=(Rp-rrp*(1+K1^2-2*K1*cos(t))^(-1/2))*cos((1-Zp/Zc)*t) - (a-K1*rrp*(1+K1^2-2*K1*cos(t))^(-1/2))*cos(Zp/Zc*t);
     y(end+1)=(Rp-rrp*(1+K1^2-2*K1*cos(t))^(-1/2))*sin((1-Zp/Zc)*t) + (a-K1*rrp*(1+K1^2-2*K1*cos(t))^(-1/2))*sin(Zp/Zc*t);
    

    

end
plot(x,y)