%對(duì)數(shù)據(jù)進(jìn)行后處理的函數(shù)，此函數(shù)中輸出經(jīng)過(guò)線形插值后的主凸輪和回凸輪的向徑和壓力角的值。
function func
[sol,kk]=fun;
sol=sol;                                       %輸出沒有對(duì)回凸輪的數(shù)據(jù)按凸進(jìn)行處理以前的原始數(shù)據(jù)。
kk=kk;
%**************************************************************************
%**************************************************************************
%判定因子sequence用來(lái)決定是否對(duì)回凸輪的數(shù)據(jù)按凸輪轉(zhuǎn)角fi進(jìn)行排序或者按凸輪的向徑角thita進(jìn)行處理，
%sequence=1:按凸輪的向徑角進(jìn)行計(jì)算,sequence=2:按凸輪轉(zhuǎn)角fi進(jìn)行排序。
sequence=2;          
if sequence==1
    sol(360/kk+1,7)=sol(360/kk+1,7)+360;       %按凸輪轉(zhuǎn)角fi對(duì)回凸輪的數(shù)據(jù)進(jìn)行重新排序。
    for t=1:1:360/kk+1
        for s=t+1:1:360/kk+1
            if sol(t,7)>sol(s,7)
                m=sol(t,7:11);
                sol(t,7:11)=sol(s,7:11);
                sol(s,7:11)=m;
            end
        end
    end
else                                          %按凸輪的向徑角進(jìn)行計(jì)算。
    for w=2:1:360/kk+1
        if sol(w,7)<=sol(1,7)
            sol(w,7)=sol(w,7)+360;
        end
    end
    media=sol(1,7);
    for z=1:1:360/kk+1
        sol(z,7)=sol(z,7)-media;
    end
end        
sol=sol;                                     %輸出處理后的結(jié)果
%對(duì)數(shù)據(jù)進(jìn)行線形插值，得到所需的主凸輪和回凸輪向徑和壓力角的值。
a=0:kk:360;                                  %對(duì)主凸輪線性插值。
for dita=0:kk:360                            
    m=round(dita/kk+1);
    for v=1:1:360/kk
        while and(dita>=sol(v,2),dita<sol(v+1,2))
            jj1=sol(v,5)+((dita-sol(v,2))/(sol(v+1,2)-sol(v,2)))*(sol(v+1,5)-sol(v,5));
            kk1=sol(v,6)+((dita-sol(v,2))/(sol(v+1,2)-sol(v,2)))*(sol(v+1,6)-sol(v,6));
            break
        end
    end                                      %輸出主凸輪插值后的結(jié)果。
    jj1=jj1;
    re(m,1)=dita;
    re(m,2)=jj1;
    re(m,3)=kk1;
end
re(1,4)=sol(1,10);                           %由于回凸輪數(shù)據(jù)的第一項(xiàng)不為零，無(wú)法迭代計(jì)算，所以首先對(duì)第一項(xiàng)進(jìn)行計(jì)算。
re(1,5)=sol(1,11); 
for dita=kk:kk:360                           %對(duì)回凸輪線性插值。
    m=round(dita/kk+1);
    for v=1:1:360/kk
        while and(dita>=sol(v,7),dita<sol(v+1,7))
            jj2=sol(v,10)+((dita-sol(v,7))/(sol(v+1,7)-sol(v,7)))*(sol(v+1,10)-sol(v,10));
            kk2=sol(v,11)+((dita-sol(v,7))/(sol(v+1,7)-sol(v,7)))*(sol(v+1,11)-sol(v,11));
            break
        end
    end
    re(m,4)=jj2;                            %輸出主凸輪插值后的結(jié)果。
    re(m,5)=kk2;
end
%**************************************************************************
%**************************************************************************
%輸出主凸輪和回凸輪插值后的最后結(jié)果。
re=re
figure;                                    %輸出主凸輪和回凸輪的輪廓線。
plot(sol(:,3),sol(:,4),sol(:,8),sol(:,9))
xlabel('X軸,mm');
ylabel('Y軸,mm');
hold on
figure                                     %輸出主凸輪和回凸輪的向徑值。
plot(re(:,1),re(:,2),re(:,1),re(:,4));
xlabel('向徑角θ,Degree');
ylabel('向徑ρ,mm')
hold on;
figure;                                    %輸出主凸輪和回凸輪的壓力角值。
plot(re(:,1),re(:,3),re(:,1),re(:,5));
xlabel('向徑角θ,Degree');
ylabel('壓力角α,mm');




%子程序fun是用來(lái)求出主凸輪和回凸輪的坐標(biāo)值和壓力角值的。子程序的接口為sol和kk，sol為%(360/ kk+1)*11的矩陣。kk為凸輪轉(zhuǎn)角的間隔值。
function [sol,kk]=fun
%求出凸輪轉(zhuǎn)角為0時(shí)D點(diǎn)的初始參數(shù)DY0和擺桿的初始角度a0。
format short;
kk=1;
l1=216;
l2=280;
l3=185;
l4=56;
bx=0;
by=0;
bx1=0;
by1=0;
dx=175;
a0=acos((l1^2+l3^2-l4^2)/(2*l1*l3));              %擺桿的初始轉(zhuǎn)角a0。
cx0=l1*cos(a0);
cy0=l1*sin(a0);
a1=acos((cx0-dx)/l2);
dy0=cy0+l2*sin(a1);                           %D點(diǎn)的初始坐標(biāo)DY0。
result=zeros(360/kk+1,7);
xx=zeros(360/kk+1,2);
yy=xx;
% *********************************************************************************
% *********************************************************************************
%循環(huán)開始求出凸輪的坐標(biāo)值。
for fi=0:kk:360
    if and(fi>=0,fi<120)               %求出D點(diǎn)的位移和速度。
        s=25*(10*(fi/120)^3-15*(fi/120)^4+6*(fi/120)^5);
        s1=(25*3/(2*pi))*(30*(fi/120)^2-60*(fi/120)^3+30*(fi/120)^4);
    elseif and(fi>=120,fi<250)
        s=25;
        s1=0;
    elseif and(fi>=250,fi<330)
        s=25-25*(10*((fi-250)/80)^3-15*((fi-250)/80)^4+6*((fi-250)/80)^5);
        s1=-1*((25*9)/(4*pi))*(30*((fi-250)/80)^2-60*((fi-250)/80)^3+30*((fi-250)/80)^4);
    else
        s=0;
        s1=0;
    end
    [bita,bita1]=link(s,s1,dy0,a0);       %調(diào)用link子程序求出擺桿角度bita和角加速度bita1。
    res=cam(fi,bita,bita1);            %調(diào)用cam子程序求出凸輪在某一凸輪轉(zhuǎn)角fi坐標(biāo)值。
    result(fi/kk+1,:)=res;              %通過(guò)循環(huán)得到0到360度主凸輪和回凸輪的坐標(biāo)值。
end
xx(:,1)=result(:,2);
xx(:,2)=result(:,5);
yy(:,1)=result(:,3);
yy(:,2)=result(:,6);
tht=zeros(360/kk+1,2);
for u=1:1:2                          %求出主凸輪和回凸輪的向徑角θ。
    for t=1:1:360/kk+1
        if xx(t,u)>=0 
            if yy(t,u)>=0
                tht(t,u)=(atan(yy(t,u)/xx(t,u)))*180/pi;
            else
                tht(t,u)=(2*pi+atan(yy(t,u)/xx(t,u)))*180/pi;
            end
        else
            tht(t,u)=(pi+atan(yy(t,u)/xx(t,u)))*180/pi;
        end
    end
end
sol=zeros(360/kk+1,11); 
sol(:,1)=result(:,1);                              %輸出主凸輪和回凸輪的轉(zhuǎn)角φ。
sol(:,2)=tht(:,1);                                %輸出主凸輪的向徑角θ。
sol(:,3)=result(:,2);                              %輸出主凸輪的X坐標(biāo)。
sol(:,4)=result(:,3);                              %輸出主凸輪的Y坐標(biāo)。
sol(:,5)=sqrt(result(:,2).^2+result(:,3).^2);            %輸出主凸輪的向徑值ρ。
sol(:,6)=result(:,4);                              %輸出主凸輪的壓力角α。
sol(:,7)=tht(:,2);                                %輸出回凸輪的向徑角θ。
sol(:,8)=result(:,5);                              %輸出回凸輪的X坐標(biāo)。
sol(:,9)=result(:,6);                              %輸出回凸輪的Y坐標(biāo)。
sol(:,10)=sqrt(result(:,5).^2+result(:,6).^2);           %輸出回凸輪的向徑值ρ。
sol(:,11)=result(:,7);                             %輸出回凸輪的壓力角α。
sol=sol;
% **********************************************************************************
% **********************************************************************************
%LINK子程序求出擺桿的角位移dita和角加速度dita1。
function [bita,bita1]=link(s,s1,dy0,a0)
    bx=0;
    by=0;
    bx1=0;
    by1=0;
    dx=175;
    dx1=0;
    dy=dy0+s;
    dy1=s1;
    l1=216;
    l2=280;
    d=sqrt((dy-by)^2+(dx-bx)^2);
    dita=atan((dy-by)/(dx-bx));
    gama=acos((d^2+l1^2-l2^2)/(2*d*l1));
    bit=dita-gama;
    cx=l1*cos(bit);
    cy=l1*sin(bit);
    bita=bit-a0;
    bita1=((dx1-bx1)*(cx-dx)+(dy1-by1)*(cy-dy))/((cy-dy)*(cx-bx)-(cy-by)*(cx-dx));
% **********************************************************************************
% ********************************************************************************
% 計(jì)算凸輪的坐標(biāo)值的子程序
function res=cam(fi,bita,bita1)
    fi=(fi/180)*pi;                   
    b=bita;                          
    b1=bita1;
    r01=130;
    r02=75.5;
    r1=21;
    r2=17.5;
    ll1=185;
    ll2=193;
    xy=zeros(1,5);
    res=zeros(1,5);
    for ii=1:1:2
        if ii==1
            l=ll1;
            r=r1;
            r0=r01;
            m=1;
            n=-1;
            k=1;
        else
            l=ll2;
            r=r2;
            r0=r02;
            m=-1;
            n=-1;
            k=1;
            ii=ii+2;
        end
        a=sqrt(130^2+185^2);
        bita=acos((a^2+l^2-r0^2)/(2*a*l));
        labu=atan((l-a*cos(bita))/(2*a*sin(bita)));
        af=atan((l*(1-n*k*b1)-a*cos(bita+k*b))/(a*sin(bita+k*b)));
        x=l*sin(labu+k*b-n*fi)-m*r*cos(labu+k*b-n*fi-af)-a*sin(labu-bita-n*fi);
        y=(l*cos(labu+k*b-n*fi)+m*r*sin(labu+k*b-n*fi-af)-a*cos(labu-bita-n*fi))*n;
        aa=abs(af*180/pi); 
        res(1,1)=fi/pi*180;
        res(1,ii+1)=x;
        res(1,ii+2)=y;
        res(1,ii+3)=aa;
    end