核心代码，记得看注释，哪些地方要改注释已经标明
clc;clear;close all;
%%
pic_num=1;
%滚子直动从动件盘形凸轮
phi=120;
%回程运动角（单位：角度，下同，以下参数需要根据题目修改）
phis=60;
%远休止角
phi_=120;
%推程运动角
phis_=60;
%近休止角
h=14;
%升程
r0=70;
%基圆半径
omega=0.008;
%角速度
r1=12;
%滚子半径
e=-12;
%偏置，正值为右偏置
%从角度转化为弧度
phi=phi*pi/180;
phis=phis*pi/180;
phi_=phi_*pi/180;
phis_=phis_*pi/180;
%一些要给定的初始条件，不用修改
s0=(r0^2-e^2)^(1/2);
beta0=atan(e/s0);
beta=0;
i=1;thetaR=0;
%以下内容均不用修改
while(thetaR<2*pi)
%推程运动
if thetaR-beta+beta0<=phi
s=spatial(thetaR-beta+beta0,h,phi);
v=velocity(omega,h,phi,thetaR-beta+beta0);
beta=atan(e/(s0+s));
p=sqrt((s+s0)^2+e^2);
d(i)=p*cos(beta);
alpha(i)=atan(abs(v/omega+e)/(s+s0));
theta(i)=thetaR;
randa=alpha(i)+beta;
pT=sqrt(p^2+r1^2-2*p*r1*cos(randa));
delta=atan(r1*sin(randa)/(p-r1*cos(randa)));
thetaT=thetaR+delta;
x(i)=pT*cos(thetaT+beta0);
y(i)=pT*sin(thetaT+beta0);
i=i+1;
end
%远休止
if thetaR-beta+beta0>=phi&&thetaR-beta+beta0<=phis+phi
s=14; %这边的s必须要改！！！改成升程！！！
v=0;
beta=atan(e/(s0+s));
p=sqrt((s+s0)^2+e^2);
d(i)=p*cos(beta);
alpha(i)=atan(abs(v/omega+e)/(s+s0));
theta(i)=thetaR;
randa=alpha(i)+beta;
pT=sqrt(p^2+r1^2-2*p*r1*cos(randa));
delta=atan(r1*sin(randa)/(p-r1*cos(randa)));
thetaT=thetaR+delta;
x(i)=pT*cos(thetaT+beta0);
y(i)=pT*sin(thetaT+beta0);
i=i+1;
end
%回程运动
if thetaR-beta+beta0>=phis+phi&&thetaR-beta+beta0<=phis+phi+phi_
s=spatial_(h,thetaR-beta+beta0,phis,phi_,phi);
v=velocity_(h,omega,phi_,thetaR-beta+beta0,phi,phis);
beta=atan(e/(s0+s));
p=sqrt((s+s0)^2+e^2);
d(i)=p*cos(beta);
alpha(i)=atan(abs(v/omega+e)/(s+s0));
theta(i)=thetaR;
randa=alpha(i)+beta;
pT=sqrt(p^2+r1^2-2*p*r1*cos(randa));
delta=atan(r1*sin(randa)/(p-r1*cos(randa)));
thetaT=thetaR+delta;
x(i)=pT*cos(thetaT+beta0);
y(i)=pT*sin(thetaT+beta0);
i=i+1;
end
%近休止
if thetaR-beta+beta0>=phis+phi+phi_&&thetaR-beta+beta0<=phis+phi+phi_+phis_
s=0;
v=0;
beta=atan(e/(s0+s));
p=sqrt((s+s0)^2+e^2);
d(i)=p*cos(beta);
alpha(i)=atan(abs(v/omega+e)/(s+s0));
theta(i)=thetaR;
randa=alpha(i)+beta;
pT=sqrt(p^2+r1^2-2*p*r1*cos(randa));
delta=atan(r1*sin(randa)/(p-r1*cos(randa)));
thetaT=thetaR+delta;
x(i)=pT*cos(thetaT+beta0);
y(i)=pT*sin(thetaT+beta0);
i=i+1;
end
thetaR=thetaR+omega;
omega=omega+0.002;%加速度（这个需要根据题目要求修改，+0.002代表匀加速的加速度为0.002）
end
x(i)=x(1);
y(i)=y(1);
[max,position]=max(alpha); %输出最大压力角及其所在位置
disp(max*180/pi);
disp(theta(position)); %输出最大压力角对应的转角
min=min(d);
disp(min); %输出最小曲率半径
plot(x,y)%绘图
辅助的四个函数，以下四个函数为简谐运动，其他题目要求的等加速或者正弦加速公式请根据高等教育出版社的《机械设计基础》第七版第40,41页来修改这四个函数
一、推程位移函数，distance为化为弧度制的推程转角，height为升程，phine为phine转角，distance_为化为弧度制的回程转角，omega为角速度，distance_s为化为弧度制的远休止角转角
function output = spatial(phine,height,distance)
%UNTITLED 此处显示有关此函数的摘要
% 此处显示详细说明
output=height*(1-cos(phine*pi/distance))/2;
end
二、推程速度函数，
function output = velocity(omega,height,distance,phine)
output=height*pi*omega*sin(pi*phine/distance)/(2*distance);
end
三、回程位移函数
function output = spatial_(height,phine,distance_s,distance_,distance)
%UNTITLED 此处显示有关此函数的摘要
% 此处显示详细说明
output=height*(1+cos(pi/distance_*(phine-distance-distance_s)))/2;
end
四、回程速度函数
function output = velocity_(height,omega,distance_,phine,distance,distance_s)
%VELOCITY_ 此处显示有关此函数的摘要
% 此处显示详细说明
output=-height*pi*omega/(2*distance_)*sin(pi/distance_*(phine-distance-distance_s));
end