四连杆之MATLAB相关程式.docx

四連桿之 MATLAB 相關程式

-26-

台大生機系 馮丁樹教授


4-25-2002

第三章中之四連桿分析可以參考相關資料，本節則針對四連桿之動作程式加 以說明。目前所設計之程式有 f4bar.m、drawlinks.m、 fb_angle_limits.m、 drawlimits.m 等四個程式，其功能分別說明如下：

圖一、四連桿之關係位置及各桿名稱

一、f4bar 函數：
f4bar 函數之呼叫格式如下：
function [values,form] = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)

輸入參數：

>r(1:4) = 各桿之長度，r(1)為固定桿，其餘分別為曲桿、結合桿及被動 桿。
>theta1 = 第一桿之水平角，或為四連桿之架構角，以角度表示。
>theta2 = 驅動桿之水平夾角，以角度表示。一般為曲桿角，但若為結合

桿驅動，則為結合桿之水平夾角。
>td2 = 驅動桿(第二桿或第三桿)之角速度(rad/sec)。
>tdd2 = 驅動桿(第二桿或第三桿)之角加速度(rad/sec^2)。
>sigma = +1 or -1. 組合模式，負值表示閉合型，正值為分支型，但有 時需視實際情況而定。
>driver = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿)

輸出變數：
>form = 組合狀態， 0 ：表示無法組合； 1：可以正確組合

>values = 輸出矩陣，其大小為 4 X 7，各行之資料分配如下：

1
2(deg)
3(rad/s)
4(rad/s2)
5
6
7
I
桿 1 位置
θ1
ω1
α1
VQ
|VQ|
∠VQ
II
桿 2 位置
θ2
ω2
α2
VP
|VP|
∠VP
III
桿 3 位置
θ3
ω3
α3
AQ
|AQ|
∠AQ
IV
桿 4 位置
θ4
ω4
α4
AP
|AP|
∠AP

其中第一行之連桿位置向量，屬於單桿的位置向量，以格式以複數表示。第二 行為各桿之水平夾角，以度表示；第三及第四行為各桿之角度速度及角加速度， 以單位時間之弧度表示。第五至七行則為 P 點與 Q 點之速度與加速度量，第五 行為向量，第六行為絕對量，第七行為夾角，以度數表示。

值得一提的是第一行之向量表示法屬於複數之型式，故若要得到其絕對值僅 需在 MATLAB 指令檔中，以 abs()這一個函數指令即可求得，而以函數 angle() 則可求得其夾角，雖然第二行與第七行之輸出亦有相對應之夾角。

例一：為第二桿為驅動桿，r=[3 2 4 2]，theta1=0;theta2=60;td2=10; tdd2=0;sigma=-1;driver=0(crank)

>>[val,form]=f4bar([3 2 4 2],0,60,10,0,-1,0)

val = Columns 1 through 3

3
0

0
1 + 1.7321i
60

10
3.8682 - 1.0182i
-14.746

5.4078
1.8682 + 0.71389i
20.913

16.549
Columns 4 through 6

0
1 +
1.7321i
2
0
1.8682 +
0.71389i
2
-127.58
173.21 -
100i
200
-236.27
364.19 -
953.09i
1020.3
Column 7

60
20.913
-30

-69.087
form = 1 (表示可以組合)

本例中，有框線者表示其為輸入值。但第一行則已經轉換為複數型式。未來 複數型式要轉為 x-y 座表時，只要使用函數 real()及 imag()兩指令，即可進行轉 換。

例二：為第三桿(coupler)為驅動桿，r=[3 2 4 2]，theta1=0; theta2=60; td2=10; tdd2=0; sigma=-1; driver=1(crank)

>>[val,form]=f4bar([3 2 4 2],0,60,10,0,-1,1)

val =

Columns 1 through 3

3

0

0
1.3321 -
1.4919i
-48.239

-8.9487

2 +
3.4641i
60

10
0.33205 +
1.9722i
80.443

24.333

Columns 4 through 6

0
1.3321 -
1.4919i
2
-582.55
0.33205 +
1.9722i
2
0
-988.55 -
882.66i
1325.3
496.46
188.64 -
31.759i
191.29

Column 7

-48.239
80.443
-138.24
-9.5568
form = 1

程式內容：f4bar.m

function [values,form] = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)

%function [values,form] = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)
% program designed by Din-sue Fon, NTU, revised from Waldron's
% This function analyzes a four-bar linkage when the crank is the
% driving link. The input values are:
% theta1,theta2 are angles in degrees
%r(1) = length of vector 1 (frame)
%r(2) = length of vector 2 (crank)

%r(3) = length of vector 3 (coupler)
%r(4) = length of vector 4 (rocker or slider offset)
%td2 = crank or coupler angular velocity (rad/sec)
%tdd2 = crank or coupler angular acceleration (rad/sec^2)
%sigma = +1 or -1. Identifies assembly mode
%driver = 0 for crank as driver; 1 for coupler as driver
% The results are returned in the vector "values". The answers are
% stored in values according to the following:

%values (1:4,1) = link position
%values (1:4,2) = link angles in degrees
%values (1:4,3) = link angular velocities
%values (1:4,4) = link angular accelerations
%values (1,5) = velocity of point Q
%values (2,5) = velocity of point P
%values (3,5) = acceleration of point Q
%values (4,5) = acceleration of point P
%vakyes (4,6) =absolute values of values(:,4)
%vakyes (4,7) =angles in degrees of values(:,4)
%form = assembly flag. If form = 0, mechanism cannot be
% assembled.

%convert input data

values=zeros(4,7);
% if coupler is the driver, interchange the vetor 3 & 2
% if theta2>180|theta2<0, sigma=-sigma;end
if driver==1,
r=[r(1) r(3) r(2) r(4)];
end
rr=r.*r; fact=pi/180; theta=zeros(4,1); td=zeros(4,1); tdd=zeros(4,1);
theta(1:2)=[theta1 theta2]*fact; t1=theta(1);
tx=theta(2); s1=sin(t1); c1=cos(t1); sx=sin(tx); cx=cos(tx);

% position calculations
A=2*r(1)*r(4)*c1-2*r(2)*r(4)*cx; C=rr(1)+rr(2)+rr(4)-rr(3)-2*r(1)*r(2)*(c1*cx+s1*sx); B=2*r(1)*r(4)*s12*r(2)*r(4)*sx;
arg=B*B-C*C+A*A;
if (arg>=0)
form=1;

% Check for the denominator equal to zero

if abs(C-A)>=eps

t4=2*atan((-B+sigma*sqrt(arg))/(C-A)); s4=sin(t4);
c4=cos(t4);
t3=atan2((r(1)*s1+r(4)*s4-r(2)*sx),(r(1)*c1+r(4)*c4-r(2)*cx)); s3=sin(t3);
c3=cos(t3);
elseif abs(C-A) If the denominator is zero, compute theta(3) first

A=-2*r(1)*r(3)*c1+2*r(2)*r(3)*cx;
B=-2*r(1)*r(3)*s1+2*r(2)*r(3)*sx; C=rr(1)+rr(2)+rr(3)-rr(4)-2*r(1)*r(2)*(c1*cx+s1*sx); arg=B*B-C*C+A*A;
if (arg>=0)
t3=2*atan((-B-sigma*sqrt(arg))/(C-A)); s3=sin(t3);
c3=cos(t3);
t4=atan2((-r(1)*s1+r(3)*s3r(2)*sx),(-r(1)*c1+r(3)*c3+r(2)*cx)); s4=sin(t4);
c4=cos(t4);

end
end

theta(3)=t3; theta(4)=t4;

%velocity calculation
td(2)=td2;
AM=[-r(3)*s3, r(4)*s4; -r(3)*c3, r(4)*c4];
BM=[r(2)*td(2)*sx; r(2)*td(2)*cx]; CM=AM\BM;
td(3)=CM(1);
td(4)=CM(2);

%acceleration calculation

tdd(2)=tdd2;

BM=[r(2)*tdd(2)*sx+r(2)*td(2)*td(2)*cx+r(3)*td(3)*td(3)*c3-r(4)*td(4)*td(4)*c 4;...

r(2)*tdd(2)*cx-r(2)*td(2)*td(2)*sx-r(3)*td(3)*td(3)*s3+r(4)*td(4)*td(4)*s4]; CM=AM\BM;
tdd(3)=CM(1);
tdd(4)=CM(2);
%store results in array values
% coordinates of P and Q
if driver==1,
r=[r(1) r(3) r(2) r(4)];
c2=c3;c3=cx;s2=s3;s3=sx;

td=[td(1) td(3) td(2) td(4)];
tdd=[tdd(1) tdd(3) tdd(2) tdd(4)]; theta=[theta(1) theta(3) theta(2) theta(4)];
else
c2=cx;s2=sx;
end
for j=1:4,
values(j,1)=r(j).*exp(i*theta(j)); values(j,2)=theta(j)/fact; values(j,3)=td(j);
values(j,4)=tdd(j);
end % position vectors values(1,5)=r(2).*exp(i*theta(2));%velocity for point Q values(2,5)=r(4).*exp(i*theta(4));%velocity for point P
values(3,5)=i*r(2).*(tdd(2)-td(2).*td(2)).*exp(i*theta(2));%accel of Q
values(4,5)=i*r(4).*(tdd(4)-td(4).*td(4)).*exp(i*theta(4));%accel of P
for j=1:4,
values(j,6)=abs(values(j,5)); %absolute values for values(:,4)
values(j,7)=angle(values(j,5))/fact; %angles for values(:,4)
end

%find the accelerations
else
form=0;
if driver==1,
r=[r(1) r(3) r(2) r(4)];
for j=1:4, values(j,1)=r(j).*exp(i*theta(j));end % position vectors

end
end

二、drawlinks 函數

有了 f4bar.m 之計算，只要連桿組可以組立(form = 1)，即表示每一個相對應 連桿的角度均可獲得，因而可以進行繪出其相關位置。程式 drawlinks 之目的即 是利用 MATLAB 繪製四連桿之相關位置。故此程式本身會呼叫 f4bar.m 函數以 計算四連桿之向量位置，然後繪圖。在其呼叫 f4bar.m 函數時，其角速度及角加 速度均設為零。其呼叫格式如下：

function drawlinks(r,th1,th2,sigma,driver)

其輸入各式與 f4bar.m 大體相同，茲說明如下：

>>r(1:4) = 各桿之長度，r(1)為固定桿，其餘分別為曲桿、結合桿及被動桿。
>>theta1 = 第一桿之水平角，或為四連桿之架構角，以角度表示。
>>theta2 = 驅動桿之水平夾角，以角度表示。一般為曲桿角，但若為結合桿

驅動，則為結合桿之水平夾角。
>>sigma = +1 or -1. 組合模式，負值表示閉合型，正值為分支型，但有時 需視實際情況而定。
>>driver = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿)

程式中，在呼叫 f4bar.m 過後，僅取其 values 中第一欄中之數據，即 r1、r2、 r3 及 r4 等之複數型式資料。構成基本的 O、R、P、Q 等四點。O 為原點，R 及 Q 為 r2 與 r4 的向量點，直接使用即可，而 P 點則為兩桿組合點，即為 r1+r4 或 r2+r3 的向量和，本程式取其前者。

例三、第二桿為驅動桿

>>drawlinks([3 2 4 2],0,60,-1,0)

2.5

2

P

1.5

1

Q

0.5

0 O R

-0.5

-1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

圖二、四連桿之繪圖

其繪出之四連桿為如圖二。黑色為第一桿，藍色為第二桿，紅色為第三桿， 綠色為第四桿。

例二、第三桿為結合桿(coupler)

>>drawlinks([3 2 4 2],0,60,-1,1)

Q

1.5

1

0.5

0 O R

-0.5

-1

P
-0.5 0 0.5 1 1.5 2 2.5 3 3.5

圖三、以結合桿為驅動桿(r =[3 2 4 2])

圖三即為所得之答案，此時四連桿為分支型(branch)，因為目前之情況無法轉 為閉合型，即使將 sigma 值變號，仍為分支型，如：

>>clf；drawlinks([3 2 4 2],0,60,1,1)

Q

1

0.5

0 O R

-0.5

-1

-1.5

P
-0.5 0 0.5 1 1.5 2 2.5 3 3.5

圖四、當 sigma=1 時並以結合桿為驅動桿(r =[3 2 4 2])

利用 drawlinks 亦可繪出各種角度之圖型，可以作為四連桿運動過程之觀察， 相當方便，例如：將 theta2 改以每 20 度繪一次，由繪至 360 度。此時必須使用 迴圈的繪製方式來達成：

>>clf; for i=10:20:360,drawlinks([1 3 3 4],0,i,-1,0); end

Q Q
Q
Q
3 Q P
P P
Q
Q P P
2

QP P
Q
1

QP P

0 O R
Q
QP P

-1

P P
Q
-2 Q
P P
Q
P P
-3 P

Q Q

Q
-4 -3 -2 -1 0 1 2 3 4 5 6

圖五、多重位置之四連桿運動情形(r =[1 3 3 4])

其結果如圖五，當然若以第三桿(聯結桿)為驅動桿時，亦可獲得同樣的結 果，此時其他數值暫不改變，僅就參數最後一項改為 1，例如：

>>clf;for i10:20:360,drawlinks([1 3 3 4],0,i,-1,1);end

Q Q
Q
Q
3 Q P
P P
P Q
P
Q
P

P Q
P Q
1 P

Q P
P
0 O R
P Q
Q

-1 P
P
Q
P
-2 Q
P
Q P
-3 P

Q Q
Q
-4 -3 -2 -1 0 1 2 3 4 5 6

圖六、以第三桿為驅動桿之情形(r =[1 3 3 4])

圖五與圖六不同之處在於第二桿(曲桿)的角度分佈。前者係由曲桿作均勻分 佈，故可以看出 360 度的均勻分佈；而圖六則為不均均的分佈情況，但第三桿的 角度變化則應為均勻的，雖然無法立即由圖六看得出來。
當 sigma 變號時，亦可看出其不同的轉動方式，如下：
>>clf; for i=10:20:360,drawlinks([1 3 3 4],0,i,1,1); end

Q

Q Q

3 P
Q P P
2 Q
P Q
P
1 P

Q
P Q
0 O R
P
Q P

-1 P
QP Q
P

P
Q P
P Q
P P
-3 Q P
Q Q
Q Q
-4 -3 -2 -1 0 1 2 3 4 5 6

圖七、當 sigma 變號(為正時) (r =[1 3 3 4])

drawlinks 程式內容
function drawlinks(r,th1,th2,sig ma,driver)
%
%function drawlinks(r,th1,th2,sigma,driver)
%draw the positions of four-bar links
%will call f4bar.m funcion
%designed by Din-Sue Fon, NTU
%
%r: row vector for four links
%th1: frame angle
%th2: crank angle or couple angle
%sigma: assembly mode
%driver: 0 for crank, 1 for coupler
%clf;
[r b]=f4bar(r,th1,th2,0,0,sigma,driver); r(3,1)=r(1,1)+r(4,1);
rx=real(r(:,1));rx(4)=0;
ry=imag(r(:,1));ry(4)=0; if b==1
plot([0 rx(1)],[0 ry(1)],'k-','LineWidth',4); hold on;
if driver==0
plot([0 rx(2)][0 ry(2)],'b-','LineWidth',1.5);
plot([rx(2) rx(3)],[ry(2) ry(3)],'r-','LineWidth',2); else
plot([0 rx(2)],[0 ry(2)],'r-','LineWidth',2);
plot([rx(2) rx(3)],[ry(2) ry(3)],'b-','LineWidth',1.5);
end
plot([rx(1) rx(3)],[ry(1) ry(3)],'-g'); plot(rx,ry,'bo');
text(0,0,' O');text(rx(1),ry(1),' R')
text(rx(2),ry(2),' P');text(rx(3),ry(3),' Q'); else
fprintf('Combination of links fail at degrees %6.1f\n',th2);
end
axis equal grid on

三、drawlimits 函數

四連桿之迴轉過程，能完全迴轉的情況仍然很少，有些時候無法獲得完整一 圈的迴轉。亦即依葛列夫定理四連桿之第一或第二類類型決定，前者為完整迴轉 型，後者則有迴轉角度之限制，這些限制因四連桿長度決定之。
四連桿迴轉過程中，有可能其中兩桿會連成一線，或重疊成一線，前者若成 立時即變成三角形，後者若重疊時，亦會構成另一個三角形。理論上連桿構成 三角形應不會有相對運動，故可稱為四連桿之運動極限。由這兩個極端位置，可 以知道四連桿之最終運動限制。

A. 第二桿為驅動桿時

在數學上，表示這兩個狀況之方法可以利用下列二種不等式進行測試：

r1+r2|r3-r4|

而由其不等式之方向，可以構成四種狀況，並進而求得該狀況之角度。下面 為第二桿為驅動桿時之四種情況：

(1) 當 r1+r2≦r3+r4，|r1-r2|≧|r3-r4| 時 θmin =0 ，θmax =2π

r3
r4
r2

r1

(2) 當 r1+r2≧r3+r4，|r1-r2|≧|r3-r4| 時

-1 2 2 2

θmin = - cos
{[r1 +r2 -(r3+r4) ]/(2 r1 r2)}，

-1 2 2 2

θmax = cos
{[r1 +r2 -(r3+r4) ]/(2 r1 r2)}

r3
r2 Amax r4
r1
Amin

(3) 當 r1+r2≧r3+r4，|r1-r2|≦|r3-r4| 時

-1 2 2 2

θmin = cos

{[r1 +r2 -(r3-r4) ]/(2 r1 r2)}，

-1 2 2 2

θmax = cos
{[r1 +r2 -(r3+r4) ]/(2 r1 r2)}

r3
r2 Amax
Amin

r4

r1

(4) 當 r1+r2≦r3+r4，|r1-r2|≦|r3-r4| 時

-1 2 2 2

θmin = cos
{[r1 +r2 -(r3-r4) ]/(2 r1 r2) }，

-1 2 2 2

θmax = 2π- cos
{[r1 +r2 -(r3-r4) ]/(2 r1 r2)}

Amax
Amin

r2 r3

r1
r4

B. 第三桿為驅動桿

第三桿結合桿為驅動桿時，則仍然取決於四連桿屬於葛列斯荷(Grashof)一型或二 型。若屬一型連桿，則當第三桿 r3 為最短桿時，第三桿可以作 360 度迴轉。其 餘之限制條件雖不如以第二桿為驅動桿者，但其極;限狀況是當第二桿與第四桿 相平行時，變成無法繼續迴轉，除非它是處於平行四邊形。將四種情況依下列二 不等式之情況加以分類，在這些分類中，若兩式均為等號時，則應歸屬於第五類：

r1+r3
|r1-r3|>|r2-r4|

(5) 當 r1+r3≦r2+r4，|r1-r3|≧|r2-r4| 時 θmin =0 ，θmax =2π

r3 r4

r2

r1

(6) 當 r1+r3≧r2+r4，|r1-r3|≧|r2-r4| 時

-1 2 2 2

θmin = - cos
{[r1 +r3 -(r2+r4) ]/(2 r1 r3)} ，

-1 2 2 2

θmax = cos
{[r1 +r3 -(r2+r4) ]/(2 r1 r3)}

r2 r4
Amin

r1
r3 Amax

(7) 當 r1+r3≧r2+r4，|r1-r3|≦|r2-r4| 時

-1 2 2 2

θmin = cos
{[r1 +r3 -(r2-r4) ]/(2 r1 r3)} ，

-1 2 2 2

θmax = cos
{[r1 +r3 -(r2+r4) ]/(2 r1 r3)}

r2 r4
Amin

r1
r3 Amax

(8) 當 r1+r3≦r2+r4，|r1-r3|≦|r2-r4| 時

-1 2 2 2

θmin = cos
{[r1 +r3 -(r2-r4) ]/(2 r1 r3) }，

-1 2 2 2

θmax = 2π- cos
{[r1 +r3 -(r2-r4) ]/(2 r1 r3)}

Amin

r1
r3 r4

Amax r2

C. fb_angle_limits 函數

觀察上面討論之四個極限角度，可以寫一組程式進行計算。由於以第三桿驅 動與第二桿驅動，在計算上僅是將其中之 r2 與 r3 之位置對調即可。為尋找上述 極限角度θmin、θmax，可用函數 fb_angle_limits 進行尋找，其格式如下：

function [Qstart, Qstop]=fb_angle_limits(r,Q1,driver)

其中輸入項目有：
r = 四連桿之長度向量，其定義與前函數相同。
Q1 = 第一桿之夾角，角度表示(deg)。
driver = 驅動模式(=0 第二桿驅動； =1 第三桿驅動)。 而輸出項為兩個角度：
Qstart = 驅動桿(第二桿或第三桿)之最低限角度 (deg) Qstop = 驅動桿(第二桿或第三桿)之最高限角度 (deg)

drawlimits 函數則是呼叫 fb_angle_limits.m 函數，然後將其極限位置繪出。其輸 入項目與 drawlinks 函數相同。

D.執行例(第二桿為驅動桿)

(1) drawlimits([1 2 3 4],0,-1,0)

Qstart = 3.6e-005 Qstop = 360

O

sRs12==03.060.0

pP

1.5

1

0.5

0 qQ

-0.5

-1

-1.5

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

圖八、第二桿驅動，屬全迴轉之情形(r = [1 2 3 4])

本例中，無論第二桿或第四桿，均可完全迴轉。

(2) drawlimits([3 5 2 1],0,1,0)
Qstart = -33.557
Qstop = 33.557

p 2.5

2

1.5

s2=33.6

1 q

0.5

0 O R

-0.5

-1 Q

-1.5
s1=-33.6

-2

-2.5


P
-1 0 1 2 3 4 5

圖九、第二類型，右邊有限制，左邊無限(r =[3 5 2 1])

圖九則因 3+5≧2+1，|3-5|≧|2-1|，故屬前述之第二型限制角度。 (3)drawlimits([5 4 1 3],0,1,0)
Qstart = 22.332
Qstop = 51.318

3.5
p 3
2.5
q Q
2

1.5
s2=51.3 P

1
s1=22.3
0.5
0 O R
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

圖十、第三類型，前後角度限制(r =[5 4 1 3])

本例中，因為各桿長度符合第三類型，即 5+4≧1+3，|5-4|≦|1-3|，故其角度 限制如圖十。

(4)drawlimits([4 3 3 5],0,1,0)
Qstart = 28.955
Qstop = 331.04

Q

3

2

P

1
s1=29.0

0 O R

s2=331.0
-1

p

-2

-3

q
-2 -1 0 1 2 3 4 5 6

圖十一、第四類型(r2 驅動) (r =[4 3 3 5])

此類型與第一類型類似，但角度限制在右邊，因為各桿之長度符合該項條 件：4+3≦3+5，|4-3|≦|3-5|。如圖十一。

E.執行例(第三桿為驅動桿)

(5) drawlimits([4 5 3 5],0,-1,1)
Qstart = 3.6e-005
Qstop = 360

pP qQ
4.5
4
3.5
3

2.5
ss12==03.600.0

2
1.5
1
0.5
0 O R
-1 0 1 2 3 4 5

圖十二、第三桿與第一桿平行，但可完全迴轉(r =[4 5 3 5])

(6) drawlimits([4 3 3 3],0,-1,1)

Qstart = -117.28
Qstop = 117.28

1.5
q P

1

0.5
s1=-117.3

0 O R

-0.5


s2=117.3

-1

Q p
-1.5
0 0.5 1 1.5 2 2.5 3 3.5 4

圖十三、 第三桿之限制角度相反，屬第六情況(r =[4 3 3 3 ])

(7) drawlimits([4 4 6 1],0,-1,1)
Qstart = 26.384
Qstop = 55.771

q
0.5
0 O R
-0.5
-1 Q

-1.5
-2
s1=26.4
s2=55.8

-2.5
-3
-3.5 P
p
-2 -1 0 1 2 3 4

圖十四、不同角度但同一象限內之限制，屬第七狀況(r =[4 4 6 1])

(8) drawlimits([3 4 4 6],0,-1,1)

Qstart = 28.955
Qstop = 331.04

Q 5
4 P

3

2 s1=29.0

1

0 O R

-1

-2 s2=331.0

-3

-4 p

-5

q
-4 -2 0 2 4 6 8

圖十五、上下象限之限制角度，兩者和為 360 度(r =[3 4 4 6 ])

E.程式內容

function [Qstart, Qstop]=fb_angle_limits(r,Q1,driver)
%
%function [Qstart, Qstop]=fb_angle_limits(r,Q1,driver)
%
% Function file that calculates the initial and final angle for four-bar
% driver=0 when the crank is the driver.
% driver=1 when coupler is the driver.
% Variables
% r=linkage row vector (cm)
% Q1=frame angle(degree);
% Qstart=initial crank(or coupler) angle (deg)
% Qstop=final crank(or co