gun.m

Matlab数值解算法实现代码

function gun%(rou,thita,inte)
%θρψωαε
inte=0.5;
inte=inte/180*pi;
ther=round(2*pi/inte);
re=func1;
for i=1:1:ther
    rou(i)=re(i,:);
end
rou=rou;
thita=[0:inte:2*pi-inte];
% thita=thita/180*pi;
% inte=inte/180*pi;
lab=185;
loa=(130^2+185^2)^0.5;
rr=21;                                       %滚子半径。
rm=1.5;                                      %测量仪器的刃口半径。
lbta=1;                                      %外凸轮为1，内凸轮为－1。
yita=1;
%**************************************************************************
%**************************************************************************
%求出凸轮的法向压力角kexi。
for i=1:1:ther
    if i==1
        yy(i)=rou(i+1)*sin(thita(i+1))-rou(ther)*sin(thita(ther));
        xx(i)=rou(i+1)*cos(thita(i+1))-rou(ther)*cos(thita(ther));
    elseif i==ther
        yy(i)=rou(1)*sin(thita(1))-rou(ther-1)*sin(thita(ther-1));
        xx(i)=rou(1)*cos(thita(1))-rou(ther-1)*cos(thita(ther-1));
    else
        yy(i)=rou(i+1)*sin(thita(i+1))-rou(i-1)*sin(thita(i-1));
        xx(i)=rou(i+1)*cos(thita(i+1))-rou(i-1)*cos(thita(i-1));
    end
    if xx(i)>0
        if yy(i)>0
            kexi(i)=atan(yy(i)/xx(i))-0.5*pi;
        else
            kexi(i)=2*pi+atan(yy(i)/xx(i))-0.5*pi;
        end
    else
        kexi(i)=pi+atan(yy(i)/xx(i))-0.5*pi;
    end
end
kexi(1)=kexi(1)-2*pi;
kexi=kexi;
%**************************************************************************
%**************************************************************************
% km=(0:inte:2*pi-inte)*180/pi;
% k1=km(:);
% k2=kexi(:)*180/pi;
% resu(:,1)=k1(:,1);
% resu(:,2)=xx(:);
% resu(:,3)=yy(:);
% resu(:,4)=k2(:,1);
% resu(:,5)=rou(:);
% resu=resu
% plot(resu(:,1),resu(:,4));
% figure;
% plot(resu(:,1),resu(:,3));
%**************************************************************************
%**************************************************************************
%求出凸轮的理论轮廓线的向径值roub和向径角thitab。
for i=1:1:ther
    y1(i)=rou(i)*sin(thita(i))+lbta*(rr-rm)*sin(kexi(i));
    x1(i)=rou(i)*cos(thita(i))+lbta*(rr-rm)*cos(kexi(i));
    if x1(i)>0
        if y1(i)>0
            thitab(i)=atan(y1(i)/x1(i));
        else
            thitab(i)=2*pi+atan(y1(i)/x1(i));
        end
    else
        thitab(i)=pi+atan(y1(i)/x1(i));                                             %凸轮理论轮廓线向径角。
    end
    roub(i)=rou(i)*cos(thita(i)-thitab(i))+lbta*(rr-rm)*cos(kexi(i)-thitab(i));     %凸轮理论轮廓线向径值。
end
thitab(1)=thitab(1)-2*pi;
%**************************************************************************
%**************************************************************************
% km=(0:inte:2*pi-inte)*180/pi;
% k1=km(:);
% resu(:,1)=k1;
% resu(:,2)=y1(:);
% resu(:,3)=x1(:);
% resu(:,4)=thitab(:)*180/pi;
% resu(:,5)=roub(:);
% resu=resu
% plot(resu(:,1),resu(:,4),resu(:,1),resu(:,5));
% plot(resu(:,1),resu(:,5));
% hold on
% plot(resu(:,1),resu(:,4));
%**************************************************************************
%**************************************************************************
%求出凸轮在转动时的对应的转角fi和摆杆的转角kx、摆杆的角速度dkx_fi。
for i=1:1:ther
    beitab=acos((roub(1)^2+loa^2-lab^2)/(2*roub(1)*loa));
    beita(i)=acos((roub(i)^2+loa^2-lab^2)/(2*roub(i)*loa));
    fi(i)=thitab(i)+yita*(beitab-beita(i));                            %求凸轮的转角。
    kxb=acos((loa^2+lab^2-roub(1)^2)/(2*loa*lab));
    kxx(i)=acos((loa^2+lab^2-roub(i)^2)/(2*loa*lab));
    kx(i)=kxx(i)-kxb;                                                  %求摆杆的转角。
    lap(i)=lab*((sin(fi(i)-yita*(beitab+kxb+kx(i))-kexi(i))))/sin(fi(i)-yita*beitab-kexi(i));
    dkx_fi(i)=yita*(1-loa/lap(i));                                     %求出摆杆的角速度。
end
%**************************************************************************
%**************************************************************************
% km=(0:inte:2*pi-inte)*180/pi;
% ffi=fi(:)*180/pi;
% kkx=kx(:)*180/pi;
% ddk=dkx_fi(:);
% k1=km(:);
% resu(:,1)=k1;
% resu(:,2)=ffi;
% resu(:,3)=kkx;
% resu(:,4)=ddk;
% resu=resu
% % plot(resu(:,1),resu(:,2),resu(:,1),resu(:,3));
% % hold on
% figure
% plot(resu(:,1),resu(:,3))
% hold on
% figure
% plot(resu(:,1),resu(:,4));
%**************************************************************************
%**************************************************************************
%求摆杆的类角加速度。
for i=1:1:ther
    if i==1
        ddkx_fi(i)=(dkx_fi(i+1)-dkx_fi(ther))/(fi(i+1)-fi(ther));
    elseif i==ther
        ddkx_fi(i)=(dkx_fi(1)-dkx_fi(ther-1))/(fi(1)-fi(ther-1));
    else
        ddkx_fi(i)=(dkx_fi(i+1)-dkx_fi(i-1))/(fi(i+1)-fi(i-1));
    end
end
%**************************************************************************
%**************************************************************************
%输出结果。
km=(0:inte:2*pi-inte)*180/pi;
ffi=fi(:)*180/pi;
kkx=kx(:)*180/pi;
ddk=dkx_fi(:);
ddkk=ddkx_fi(:);
k1=km(:);
resu(:,1)=k1;   
resu(:,2)=ffi;              %凸轮的转角。
resu(:,3)=kkx;              %摆杆的转角。
resu(:,4)=ddk;              %摆杆的角速度。
resu(:,5)=ddkk;             %摆杆的角加速度。
resu=resu                   %输出结果。
plot(resu(:,1),resu(:,3));
xlabel('凸轮转角φ,Degree');
ylabel('摆杆摆角ψ,Degree');
figure
plot(resu(:,1),resu(:,4));
xlabel('凸轮转角φ,Degree');
ylabel('摆杆角速度ψ,mm/s');
figure
plot(resu(:,1),resu(:,5));
xlabel('凸轮转角φ,Degree');
ylabel('摆杆角加速度ε,mm/s*s');