% 對心直動凸輪機構壓力角的計算(調用TLYLJ.M)
disp ' '
disp '        ******** 對心直動凸輪機構壓力角的計算 ********'
disp ' '
disp '                ======== 已 知 條 件 ========'
disp ' '
rb = input('           基圓半徑(mm)   rb = ');
h = input('           推程升程(mm)    h = ');
k=h/rb;hd=pi/180;
fai = input('         推程運動角(度)  fai = ');
fprintf (1,'         運動結構系數      k = %3.4f \n',k) 
YDGL = input('   運動規律類型：等速-"ZX";等加減速-"PW";余弦加速-"JX";正弦加速-"BX" == ');
disp ' '
if YDGL=='ZX'
    disp '               ======== 等速運動（直線）規律 ========'
    fm=0;
    alfm=atan(k/(fai*hd));
elseif YDGL=='PW'
    disp '            ======== 等加減速運動（拋物線）規律 ========'
    if k<=2
        fm=fai*hd/2;
        alfm=atan(4*k/(fai*hd*(2+k)));
    elseif k>2
        fm=fai*hd/sqrt(2*k);
        alfm=atan(sqrt(2*k)/(fai*hd));
    end
elseif YDGL=='JX'
    disp '          ======== 余弦加速度運動（簡諧曲線）規律 ========'
    fm=fai*hd*acos(k/(2+k))/pi;
    alfm=atan(k*pi/(2*fai*hd*sqrt(1+k)));
elseif YDGL=='BX'
    disp '            ======== 正弦加速度運動（擺線）規律 ========'
    x=fsolve(@TLYLJ,fai*hd/2);           % 使用fsolve求解漸開線函數方程
    fm=x/pi*(fai*hd);
alfm=atan(k*(1-cos(2*pi*fm/(fai*hd)))/(fai*hd+k*fm-k*fai*hd*sin(2*pi*fm/(fai*hd))/(2*pi)));
end
fprintf (1,'              最大壓力角      alfm = %3.4f 度 \n',alfm/hd)
fprintf (1,'      最大壓力角的位置角        fm = %3.4f 度 \n',fm/hd)

% 壓力角漸開線函數
function f=TLYLJ(x)
global k                       % 定義全局變量
f=tan(x)-x-pi/k;