% 鉸鏈四桿機構運動設計(調用qbyg.m)
x0=[50 120 200 0.5];
k=1.25;                    % 行程速比系數
theta=pi*(k-1)/(k+1);      % 極位夾角
yg=250;                    % 搖桿長度
pusai=pi/6;                % 搖桿擺角
gamin=2*pi/9;              % 最小傳動角
x=fsolve(@qbyg,x0);
disp '           ********  已知條件  ********'  
fprintf (1,'         行程速比系數    k = %3.4f  \n',k)
fprintf (1,'         極位夾角    theta = %3.4f 度 \n',theta*180/pi)
fprintf (1,'         搖桿長度       yg = %3.4f mm \n',yg)
fprintf (1,'         搖桿擺角    pusai = %3.4f 度 \n',pusai*180/pi)
fprintf (1,'         最小傳動角  gamin = %3.4f 度 \n',gamin*180/pi)
disp '           ********  計算結果  ********'  
fprintf (1,'         曲柄長度        a = %3.4f mm \n',x(1))
fprintf (1,'         連桿長度        b = %3.4f mm \n',x(2))
fprintf (1,'         機架長度        d = %3.4f mm \n',x(3))
fprintf (1,'         搖桿位置角  gamin = %3.4f 度 \n',x(4)*180/pi)

% 鉸鏈四桿機構非線性參數方程組
function f=qbyg(x)
k=1.25;                    % 行程速比系數
theta=pi*(k-1)/(k+1);      % 極位夾角
yg=250;                    % 搖桿長度
pusai=pi/6;                % 搖桿擺角
gamin=2*pi/9;              % 最小傳動角
% x(1)是曲柄長度;x(2)是連桿長度;x(3)是機架長度;x(4)是搖桿初始位置角
f1=(x(2)+x(1))^2+(x(2)-x(1))^2-2*(x(2)+x(1))*(x(2)-x(1))*cos(theta)-(2*yg*sin(pusai/2))^2;
f2=yg^2+x(3)^2-2*yg*x(3)*cos(x(4))-(x(2)-x(1))^2;
f3=yg^2+x(3)^2-2*yg*x(3)*cos(x(4)+pusai)-(x(2)+x(1))^2;
f4=yg^2+x(2)^2-2*yg*x(2)*cos(gamin)-(x(3)-x(1))^2;
f=[f1;f2;f3;f4];