%對數據進行后處理的函數，此函數中輸出經過線形插值后的主凸輪和回凸輪的向徑和壓力角的值。
function func
[sol,kk]=fun;
sol=sol;                                       %輸出沒有對回凸輪的數據按凸進行處理以前的原始數據。
kk=kk;
%**************************************************************************
%**************************************************************************
%判定因子sequence用來決定是否對回凸輪的數據按凸輪轉角fi進行排序或者按凸輪的向徑角thita進行處理，
%sequence=1:按凸輪的向徑角進行計算,sequence=2:按凸輪轉角fi進行排序。
sequence=2;          
if sequence==1
    sol(360/kk+1,7)=sol(360/kk+1,7)+360;       %按凸輪轉角fi對回凸輪的數據進行重新排序。
    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                                          %按凸輪的向徑角進行計算。
    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;                                     %輸出處理后的結果
%對數據進行線形插值，得到所需的主凸輪和回凸輪向徑和壓力角的值。
a=0:kk:360;                                  %對主凸輪線性插值。
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                                      %輸出主凸輪插值后的結果。
    jj1=jj1;
    re(m,1)=dita;
    re(m,2)=jj1;
    re(m,3)=kk1;
end
re(1,4)=sol(1,10);                           %由于回凸輪數據的第一項不為零，無法迭代計算，所以首先對第一項進行計算。
re(1,5)=sol(1,11); 
for dita=kk:kk:360                           %對回凸輪線性插值。
    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;                            %輸出主凸輪插值后的結果。
    re(m,5)=kk2;
end
%**************************************************************************
%**************************************************************************
%輸出主凸輪和回凸輪插值后的最后結果。
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是用來求出主凸輪和回凸輪的坐標值和壓力角值的。子程序的接口為sol和kk，sol為%(360/ kk+1)*11的矩陣。kk為凸輪轉角的間隔值。
function [sol,kk]=fun
%求出凸輪轉角為0時D點的初始參數DY0和擺桿的初始角度a0。
format short;
kk=0.4;
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));              %擺桿的初始轉角a0。
cx0=l1*cos(a0);
cy0=l1*sin(a0);
a1=acos((cx0-dx)/l2);
dy0=cy0+l2*sin(a1);                           %D點的初始坐標DY0。
result=zeros(360/kk+1,7);
xx=zeros(360/kk+1,2);
yy=xx;
% *********************************************************************************
% *********************************************************************************
%循環開始求出凸輪的坐標值。
for fi=0:kk:360
    if and(fi>=0,fi<120)               %求出D點的位移和速度。
        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);       %調用link子程序求出擺桿角度bita和角加速度bita1。
    res=cam(fi,bita,bita1);            %調用cam子程序求出凸輪在某一凸輪轉角fi坐標值。
    result(fi/kk+1,:)=res;              %通過循環得到0到360度主凸輪和回凸輪的坐標值。
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);                              %輸出主凸輪和回凸輪的轉角φ。
sol(:,2)=tht(:,1);                                 %輸出主凸輪的向徑角θ。
sol(:,3)=result(:,2);                              %輸出主凸輪的X坐標。
sol(:,4)=result(:,3);                              %輸出主凸輪的Y坐標。
sol(:,5)=sqrt(result(:,2).^2+result(:,3).^2);      %輸出主凸輪的向徑值ρ。
sol(:,6)=result(:,4);                              %輸出主凸輪的壓力角α。
sol(:,7)=tht(:,2);                                 %輸出回凸輪的向徑角θ。
sol(:,8)=result(:,5);                              %輸出回凸輪的X坐標。
sol(:,9)=result(:,6);                              %輸出回凸輪的Y坐標。
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));
% **********************************************************************************
% ********************************************************************************
% 計算凸輪的坐標值的子程序
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