【MATLAB】求算凸轮摆杆机构的最大压力角和最小半径并画出凸轮轮廓

function  curve      %曲线方程
a=145;              %机架距离
l=125;               %摆杆长度
r0=60;               %基圆半径
rr=8;                 %滚子半径
step=0.01;          %计算步长
flag=1;             %1为凸轮顺时针转，其他为逆转
fai0=acos((a^2+l^2-r0^2)/(2*a*l))*180/pi;
delta=[0.0001:step:360];
fai=zeros(1,length(delta));
for ii=1:length(delta)
    if (delta(ii)>0 && delta(ii)<=75)
       fai(ii)=15*(delta(ii)/75-sin(2*pi*delta(ii)/75)/(2*pi));
      else if (delta(ii)>75 && delta(ii)<=85)
          fai(ii)=15;
      else if (delta(ii)>85 && delta(ii)<=150)
             fai(ii)= 15*(1-(delta(ii)-85)/65+sin(2*pi*(delta(ii)-85)/65)/(2*pi));
      else if  (delta(ii)<=360 && delta(ii)>150)
              fai(ii)=0;
          end
          end
          end
    end
end
if flag==1
fai=fliplr(fai);
end
  %理论曲线  
  x=a*sin(delta*pi/180)-130*sin((delta+fai+fai0)*pi/180);
  y=a*cos(delta*pi/180)-130*cos((delta+fai+fai0)*pi/180);
  %实际工作曲线
  dx=derivative(x);
  dy=derivative(y);
   ddelta=step*ones(1,length(delta));
   sintheta=(dx./ddelta)/sqrt(((dx./ddelta).^2+(dy./ddelta).^2));
   costheta=-(dy./ddelta)/sqrt(((dx./ddelta).^2+(dy./ddelta).^2));
   xp=x-rr*costheta;
   yp=y-rr*sintheta;
   %理论曲线最小曲率半径
   dx2=derivative(dx);
  
   dy2=derivative(dy);
   
   rou=abs((sqrt(dx.^2+dy.^2).^3)./(dy2.*dx-dy.*dx2));
   min_radius=min(rou)
   %实际工作曲线最大压力角
   dfai=derivative(fai);
   o1p=a./(1+dfai./ddelta);
   o2p=sqrt(o1p.^2+l^2-2*l*o1p.*cos((fai0+fai)*pi/180));
   beta=acos((o2p.^2+l^2-o1p.^2)./(2*l*o2p))*180/pi;
   angle_of_pressure=90-beta;
   max_angle_of_pressure=max(angle_of_pressure)
  % plot(delta,fai);
                             % A(jj,1)=min_radius;
                             % A(jj,2)=max_angle_of_pressure;
                             % jj=jj+1;
                             %    end
%                              A(:,1)
%                              A(:,2)
   plot(x,y,'.');
end

function d=derivative(b)
d=diff(b);
d(length(d)+1)=b(1)-b(length(b));
end