作業十一

我本週有上課。

我將本作業的三個小題合併成一個大程式，因為老師沒有在題目內明講偏置量e為多少，故我的程式設計成可自行輸入偏差量；由於返程的運動形式為自訂，故我的程式設計成去程及返程的運動形式都可以自行輸入。以下為程式碼以及解釋：

-------------------------------------------------------------------------------------------------

function homework9(r0,s,e,L,range,pattern,cw)
% 把本作業的三個小題整合在一個程式內，分別以三個figure呈現出來。
% uses plot_dwell.m
% uses pincam.m

% Inputs:
% cth:angle of cam, degrees
% r0:radius of base circle
% e:offset
% s:stroke
% L:length of pin
% cw:rotation direction of cam(-counterclockwise,+clockwise
% pattern = denote the type of motion used(a 3 element-row matrix)
% 1:uniform 2:parabolic 3:simple harmonic 4: cycloidal
% 5:polynomial motion
% example [4 3]
% range =the degrees the specific motion starts, eg.[90 180 240]

% 我一律設定所有圖表都是每度畫一次，使用者不得自行更改。
ctheta=0:1:359;

% r0必小於e，否則會有凸輪轉到一半從動件跟凸輪失去接觸的情況。
if r0<=e
error('Please try another set of r0 & e, since r0 cannot be smaller than e!')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%% 第一小題 %%%%%%%%%%%%%%%%%%%%%%%%%%%
% 以下是將老師的plot_dwell稍做修改後得到：

figure(1); % 第一小題，使用figure1。
clf;
[y,yy,yyy]=dwell(ctheta,range,pattern);
plot(ctheta,y*s,'b-',ctheta,yy*s,'k-',ctheta,yyy*s,'r-'); % 我僅把此行的handle拿掉，其他都跟原程式碼一樣。
legend('Displacement','Velocity','Acceleration',3)
xlabel('Elapsed Angle, degrees')
grid

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 第二小題 %%%%%%%%%%%%%%%%%%%%%%%%%%%
% 以下是將老師的pincam修改後得到：

figure(2); % 第二小題，使用figure2。
clf;
th=ctheta*pi/180;
s0=sqrt(r0*r0-e*e);

for i=1:length(ctheta)

t=th(i)*cw;
A=[cos(t) -sin(t);sin(t) cos(t)];
[ym,yy,yyy]=dwell(ctheta(i),range,pattern);
x0=s0+ym*s;
Sx=[0 x0 x0+L;e e e];
X=A\Sx;
x(i)=X(1,2);y(i)=X(2,2);

% 每十度畫一次從動件的延伸和從洞件。
% 若是每一度都畫一次，圖形會密到看不出來。
r= mod(i,10); % 計數器i可以被10整除時就畫一次。
if r == 0
line(X(1,1:2),X(2,1:2)); % 畫從動件的延伸
line(X(1,2:3),X(2,2:3),'linewidth',3,'color','red'); % 畫從動件
end

end

hold on;
plot(0,0,'ro',x,y,'k-') % 畫工作曲線
axis([-(r0+s+L+5) (r0+s+L+5) -(r0+s+L+5) (r0+s+L+5)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 第三小題 %%%%%%%%%%%%%%%%%%%%%%%%%%%
% 以下是將第二小題的程式再進行修改後所得到：

figure(3); % 第三小題，使用figure3。
clf;

% 分成有沒有偏移量e來處理
if e==0 % 無偏移量e

while 1 % 讓動畫可以無限循環
for i=1:length(ctheta)

clf; % 製作動畫所必須的更新畫面。
% 將所有點繞原點旋轉1度後畫工作曲線。
x = x.*cosd(1) - y.*sind(cw) ;
y = x.*sind(cw) + y.*cosd(1) ;
plot(0,0,'ro',x,y,'k-');
hold on;

% 畫從動件的延伸以及從動件。
plot([x(i) x(i)+L],[0 0],'red','LineWidth',3);
plot([0 x(i)],[0 0],'blue');
axis([-(r0+s+L+5) (r0+s+L+5) -(r0+s+L+5) (r0+s+L+5)]);
pause(0.01);

end
end

else % 有偏移量e

while 1 % 讓動畫可以無限循環
for i=1:length(ctheta)

clf; % 製作動畫所必須的更新畫面。

% 當cw=-1，也就是順時針旋轉時：
% i等於360時x(i+1)會讀到第361項，需修正為第1項。
if i == 360 && cw==-1
i = 0;
end
% 當cw=-1，也就是順時針旋轉時：
% i等於1時x(i-1)會讀到第0項，需修正為第360項。
if i == 1 && cw==1
i = 360;
end
% 雖然在這裡有可能會改到計數器i，但是經過測試後發現：計數器本身並不會因為i被暫時改變而計算錯！
% 例：i=1且cw=1時i被改為360，但for迴圈下次再跑進來的時候i會是2而不是361！

% 以下寫法特殊，另行說明。
R2 = ( (x(i-cw))^2 + (y(i-cw))^2 )^0.5 ;
phi = asind(y(i-cw)/R2) - asind(e/R2) ;

x = x.*cosd(phi) - y.*sind(-phi) ;
y = x.*sind(-phi) + y.*cosd(phi) ;
plot(0,0,'ro',x,y,'k-');
hold on;
% 以上寫法特殊，另行說明。

% 畫從動件的延伸線所切出的包絡圓。
p=linspace(0,2*pi,360)';
pp=e*(sin(p)+j*cos(p));
plot(pp);

% 畫從動件的延伸以及從動件。
plot([x(i-cw) x(i-cw)+L],[e e],'red','LineWidth',3);
plot([x(i-cw) 0],[e e],'blue');
axis([-(r0+s+L+5) (r0+s+L+5) -(r0+s+L+5) (r0+s+L+5)]);
pause(0.01);

end
end

end

-------------------------------------------------------------------------------------------------
homework9(15,5,5,10,[100 200 260],[2 5],-1)

（選擇偏置量=5、返程運動型式=5）

-------------------------------------------------------------------------------------------------

另行說明部份：
我在那裡用了不太直觀的寫法。假設偏差量是e，然後我想要把凸輪順時針轉動theta角，這時候不能先用pincam把圖畫出來再把圖轉theta角！理由將配合以下兩圖說明：
上面是錯誤的轉法：假設凸輪是上面的橢圓形，從動件為紅線，左圖是一開始的樣子。我想得到凸輪順時針轉了theta角（上圖theta為45度，程式內是1度）後的圖形，然後我就直接把整張凸輪的圖都轉了theta角後，再把從動件放上去，結果就會發生上圖所畫出來的有誤差的情況。

這張是正確的轉法：因為我不能直接轉theta角，所以我必須繼續把圖旋轉，直到錯誤的那張圖的誤差量等於零，總共必須旋轉phi角。

注意！若偏差量e為0，則不會有此問題。這也是為什麼e==0那部份的程式碼較簡單的原因！

phi角與theta角的關係如下：

          R2 = ( (x(i-cw))^2 + (y(i-cw))^2 )^0.5 ;
   phi = asind(y(i-cw)/R2) - asind(e/R2) ;

作業十

1.我本週有上課。

2.
由於本題使用相當多數學式，無法在blog上顯示，故我在此po上pdf連結：
http://homepage.ntu.edu.tw/~b94202029/mechanism_HW10_1.pdf
另外我也將內容抓成圖檔貼在下方：

2.
本題也使用相當多數學式，無法在blog上顯示，故我在此po上pdf連結：
http://homepage.ntu.edu.tw/~b94202029/mechanism_HW10_2.pdf
另外我也將內容抓成圖檔貼在下方：

以下是程式碼，內附解釋：
function [ICposition]=instantaneousCenterPosition(theta,r2,r3)
% ICposition代表六點瞬心位置：
% ICposition=[ X12(必為0) X23 X34 X41(設為0) X13 X24(必為0) ; Y12(必為0) Y23 Y34(必為0) Y41(設為inf，正負號與sind(theta)同號) Y13 Y24 ]
% theta請輸入degree

% 不接受360度以外角度，不接受曲桿長大於等於連結桿長。
%（曲桿長等於連結桿長時，一旦theta到達90度，則滑塊就有可能停留在原地不動，但曲桿仍然在旋轉。）
%（曲桿長大於連結長時，一旦theta到達90度，則滑塊會被抬離地面。）
if (theta>=360) || (theta<0)>=r3) || (r2<=0) || (r3<=0) ; error('Not acceptable input(s).'); end; % 先製造出ICposition以便使用。 % ICposition=[ X12(必為0) X23 X34 X41(設為0) X13 X24(必為0) ; Y12(必為0) Y23 Y34(必為0) Y41(設為inf，正負號與sind(theta)同號) Y13 Y24 ] % (14瞬心的位置使用(0,-inf)代表之) ICposition=[ 0 0 0 0 0 0 ; 0 0 0 -inf 0 0 ]; % 以下計算方法的推導過程請見pdf檔 phi = asind( (r2/r3)*sind(theta) ) ; L = r2*cosd(theta) + r3*cosd(phi) ; ICposition(1,2)=r2*cosd(theta); ICposition(2,2)=r2*sind(theta); ICposition(1,3)=L; ICposition(2,4)=-inf; ICposition(1,5)=L; ICposition(2,5)=L*tand(theta); ICposition(2,6)=L*tand(phi); % 畫出滑塊機構。閉合型或分支型的區別在這邊只會造成左右的差異而已。 % 也就是說，分支型的150度等於閉合型的30度，以此類推。 % 故我們只挑其中一種來做就可以。 [values]=drawsldlinks([0 r2 r3 0],0,theta,1,0); % 畫出尋找瞬心13和24的輔助線 line([ICposition(1,1);ICposition(1,6)],[ICposition(2,1);ICposition(2,6)],'color','m','linestyle',':','marker','p'); line([ICposition(1,2);ICposition(1,6)],[ICposition(2,2);ICposition(2,6)],'color','m','linestyle',':','marker','p'); if ICposition(2,5)~=inf & ICposition(2,5)~=-inf line([ICposition(1,2);ICposition(1,5)],[ICposition(2,2);ICposition(2,5)],'color','m','linestyle',':','marker','p'); line([ICposition(1,3);ICposition(1,5)],[ICposition(2,3);ICposition(2,5)],'color','m','linestyle',':','marker','p'); end ---------------------------------------------------------------------------------------------------------- 執行結果：

瞬心位置數據(theta=30/90/150;r2=2;r3=5)

圖形(theta=30;r2=2;r3=5)

圖形(theta=90;r2=2;r3=5)

圖形(theta=150;r2=2;r3=5)

瞬心位置數據(theta=210/270/330;r2=2;r3=5)

圖形(theta=210;r2=2;r3=5)

圖形(theta=270;r2=2;r3=5)

圖形(theta=330;r2=2;r3=5)

作業九

我本週有上課。

請就教科書中第四章第五節之偏置機構作另類分析，分析過程可採你所知的方式（包括講義中所列的方法）。運動中分以曲桿驅動及滑塊驅動的方式，並說明運動的 界限或範圍。設此機構之曲桿長Rcm ， 連桿Lcm，滑塊之偏置量為10cm等數據作分析。其中，R=10+(學號末二碼），Ｌ＝Ｒ+5 。

我的學號：b94202029
曲桿長R=39;
連桿長L=44;
滑塊偏置量e=10;

首先網路講義對於sldlink的解說有誤： r(1)為固定桿長，其長度會變化，故可設為0。 r(4)為偏置量，不是可有可無，當有偏置量時必須輸入。

另外網路講義中sld_angle_limits的程式碼沒有把所有狀況都寫進去：
第二桿為驅動桿的第2種情況沒有寫入程式碼中
應該加入以下程式碼（請插在case 0 中，if和第1個else if 中間）
------------------------------------------------------------------------------------------
elseif r3+r4>=r2 & r4>=r3 % 網路講義中的(2) 此部份是網路的程式碼中遺漏的部份
Qstart=asin((r4-r3)/r2);
Qstop=pi-asin((r4-r3)/r2);
------------------------------------------------------------------------------------------

網路講義介紹drawsldpaths的部份，其呼叫方式、輸入參數與程式內容部份，引數都不一樣。正確的應該是：
drawsldpaths(r6,th6,r,th1,td2,tdd2,sigma,npts,driver,mode)
r6:外加固定桿長度
th6:外加固定桿角度
r:各桿件長度
th1:地桿水平角度
td2:驅動桿角速度
tdd2:驅動桿角加速度
sigma：分支型or閉合型
npts：分割點數目
driver：=0，曲桿驅動；=1，結合桿驅動；=2，滑塊驅動
mode：=0，畫簡單位置圖；=1 畫所有圖表；=2畫所有圖表，但用簡單位置圖。
而固定桿就只能在曲桿上，無法挑選。

最後，move_sldpaths程式碼內有很多個"符號，造成無法執行。但現在已經修改過了。
但這個程式也無法做出連續的連桿動作。關於這部份我會在底下討論運動情形時詳述。

------------------------------------------------------------------------------------------------
我們先來檢視這個滑塊機構的極限運動情形。現在設：
r=[0 39 44 10];
theta1=0;

使用sld_angle_limits，找出驅動桿極限運動情形，並用drawsldlimits畫出其位置：
注意：當驅動桿為第二及第三桿時，其Qstart與Qstop分別表示驅動桿之最小及最大角度；但當驅動件為滑塊時，其不再是角度，而是其對應之桿一最小及最大長度，亦即r1min與r1max。
-----------------------------------------------------------------------------------------------
驅動桿為曲桿情形：

r=[0 39 44 10];
theta1=0;
driver=0;
[Qstart, Qstop]=sld_angle_limits(r,theta1,driver)

Qstart =
-60.6679
Qstop =
240.6679

sigma=1;
drawsldlimits(r,theta1,sigma,driver);

sigma=-1;
drawsldlimits(r,theta1,sigma,driver);



驅動桿為結合桿情形：

driver=1;
[Qstart, Qstop]=sld_angle_limits(r,theta1,driver)

Qstart =
-41.2306
Qstop =
221.2306

sigma=1;
drawsldlimits(r,theta1,sigma,driver);

sigma=-1;
drawsldlimits(r,theta1,sigma,driver);




驅動桿為滑塊之情形：

r=[0 39 44 10];
theta1=0;
driver=2;
[r1min, r1max]=sld_angle_limits(r,theta1,driver)

r1min =
-82.3954
r1max =
82.3954

sigma=1;
drawsldlimits(r,theta1,sigma,driver);

sigma=-1;
drawsldlimits(r,theta1,sigma,driver);


-----------------------------------------------------------------------------------------------
由以上觀察可得知：在極限角度下的滑塊桿件的狀態不會跟滑塊桿件是分支型或閉合型有關！其實sld_angle_limits的input本來就沒有sigma，所以只是把極限狀況畫出來的drawsldlimits，在正常的滑塊分析上，可能可以被改寫成不需要輸入sigma。
-----------------------------------------------------------------------------------------------


到目前為止，我們已經知道了在曲桿、結合桿、滑塊為驅動桿時，驅動桿的極限角度（或位置）以及整組桿件的狀態。
接下來我們要分析的是這個桿件如何運動。首先要釐清一點：不管是用曲桿、結合桿、滑塊當作驅動桿，整個桿件組做完一次完整的運動的歷程是不會改變的。但是這個運動歷程會和這個桿件是分支型還是閉合型有關。有關這部份，我的作業八裡面就有提過。所以以下我只挑選用曲桿作為驅動桿，討論分支型與閉合型下的桿件運動。
(此處的運動是指桿件行進的軌跡，不是桿件的速度與加速度。不同的驅動桿和不同的輸入都會造成整個桿件的速度與加速度不同。)

我先討論sigma=1的情形。畫出曲桿在各個角度下時，整個桿件的狀態：
-----------------------------------------------------------------------------------------------
以下有呼叫drawsldlinks程式者，此程式內容有更動：我將axis([-90,90,-30,30])加在axis equal之前一行，以固定整個座標軸的大小。但座標軸隨著桿件而平移是無可避免之事。
------------------------------------------------------------------------------------------------
r=[0 39 44 10];
theta1=0;
driver=0;
sigma=1;

[Qstart, Qstop]=sld_angle_limits(r,theta1,driver)
[values]=drawsldlinks(r,theta1,Qstart,sigma,driver);
pause;
pause(0.05);
clf;

for theta2=(-60):1:240;
[values]=drawsldlinks(r,theta1,theta2,sigma,driver);
pause(0,05);
clf;
end

[values]=drawsldlinks(r,theta1,Qstop,sigma,driver);



------------------------------------------------------------------------------------------------

這是從Qstart開始到Qstop停止的影片。注意到滑塊一開始是往右邊移動，後來就開始往左移動，不斷往左移動直到影片結束。那如果曲桿的的運動是從Qstart開始到Qstop、再從Qstop開始到Qstart、再從Qstart開始到Qstop......一直下去，這樣的情況下整個桿件會怎麼移動？
上一段提到滑塊在曲桿角度接近Qstop時，一直都是往左移動。所以到了Qstop的瞬間，滑塊的速度方向是往左，此時結合桿施加給滑快的力完全沒有平行方向的分量，所以滑塊在水平方向的速度在這一瞬間不會改變，也就是說滑塊在曲桿角度經過Qstop的時候會繼續往左移動。這正代表當整個桿件從Qstop開始回到Qstart的時候，桿件的狀態是sigma=-1而不是sigma=1！！老師的move_sldpaths程式的錯誤就是在這裡。

sigma=-1的情況就不用討論，只是和sigma=1的情況左右相反罷了。

所以我只要把從Qstop開始回到Qstart的運動寫下，再用while使他做無限迴圈，就可以製造出這個連桿的運動動畫。
------------------------------------------------------------------------------------------------
以下有呼叫drawsldlinks程式者，程式有更動，更動內容在上面已經提過
------------------------------------------------------------------------------------------------
r=[0 39 44 10];
theta1=0;
driver=0;

[Qstart, Qstop]=sld_angle_limits(r,theta1,driver)

% k矩陣的row1：從k(1,1)到k(1,303)是儲存逆時針旋轉的角度值，從k(1,304)到k(1,606)是儲存逆時針旋轉的角度值。
% k矩陣的row2：從k(1,1)到k(1,303)是儲存逆時針旋轉的mode，從k(1,1)到k(1,303)是儲存逆時針旋轉的mode。
k=[Qstart linspace(-60,240,301) Qstop linspace(240,-60,301) ; linspace(1,1,302) linspace(-1,-1,302)];

while 1
for l=1:604;

clf;
theta2=k(1,l);
sigma=k(2,l);
[values]=drawsldlinks(r,theta1,theta2,sigma,driver);
pause(0.01);

end
end



------------------------------------------------------------------------------------------------
作業八和作業九的最後一題，我都得到「連桿在極限角度時將在閉合與分支之間切換型態」的結論，故我大膽推測：可能所有的四連桿都有這樣的性質！！

另外，可用drawsldpaths程式，輸出速度、加速度、角速度、角加速度隨角度變化圖，另外還可附上桿件運動路徑。由於運動路徑已經可以很清楚的從上面的影片看出，所以我就只把速度、加速度、角速度、角加速度隨角度變化圖附上。
------------------------------------------------------------------------------------------------
設角速度=180/pi，角加速度=0，npts=50。

driver=0，sigma=1：
drawsldpaths(0,0,[0 39 44 10],0,180/pi,0,1,50,0,2);

driver=0，sigma=-1：
drawsldpaths(0,0,[0 39 44 10],0,180/pi,0,-1,50,0,2);

driver=1，sigma=1：
drawsldpaths(0,0,[0 39 44 10],0,180/pi,0,1,50,1,2);

driver=1，sigma=-1：
drawsldpaths(0,0,[0 39 44 10],0,180/pi,0,-1,50,1,2);

driver=2，sigma=1：
drawsldpaths(0,0,[0 39 44 10],0,180/pi,0,1,50,2,2);

driver=2，sigma=-1：
drawsldpaths(0,0,[0 39 44 10],0,180/pi,0,-1,50,2,2);

作業八

我4/26有上課。

8.1
四連桿桿長r=[4 3 3 5]，分別為固定桿、曲桿、結合桿、被動桿。 桿二為驅動桿，故linkdrive=0。 固定桿角度theta1=0(degree)；角速度td2=10(rad/s)；角加速度tdd2=0(rad/s^2)。 桿二角度假設為theta2=45(degree)。 題目沒有給定此連桿為閉合型或分支型，故mode=+1 or -1。

使用網路講義中第六章的f4bar程式找出各點的位置、速度、加速度：
------------------------------------------------------------------------------------------------
r=[4 3 3 5]; td2=10; tdd2=0; theta1=0; theta2=45; linkdrive=0;

mode=1;
[data,form] = f4bar(r,theta1,theta2,td2,tdd2,mode,linkdrive)

mode=-1;
[data,form] = f4bar(r,theta1,theta2,td2,tdd2,mode,linkdrive)
------------------------------------------------------------------------------------------------
執行結果：
data = 輸出矩陣，其大小為 4 X 7，各行之資料分配如下：

　　　 1 　   2(deg)  3(rad/s)  4(rad/s^2)    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點之速度與加速度量，第五行為向量，第六行為絕對量，第七行為夾角，以度數表示。

教授的講義在此處有誤：
1.程式輸出中沒有column7，經查看後f4bar程式的設計中也沒有設計column7，故上圖中只有6個column。
2.data(1,6)和data(2,6)為Q&的位置(以複數表示)，data(3,6)和data(4,6)恆為0，皆不是column5的絕對值。

有了r2、r3、theta2、theta3、omega2、omega3、alpfa2、alpha3，我們就可以求出q點和p點的速度、加速度（位置已求出，就是data(1,6)和data(2,6)）。
注意theta2&theta3是角度，故要使用cosd&sind函數。
-------------------------------------------------------------------------------------------------
r2=r(2);
r3=r(3);
theta2=data(2,2);
theta3=data(3,2);
omega2=data(2,3);
omega3=data(3,3);
alpha2=data(2,4);
alpha3=data(3,4);

Vq=j*r2*(omega2)*(cosd(theta2)+j*sind(theta2))
Vp=j*r2*(omega2)*(cosd(theta2)+j*sind(theta2))
+j*r3*(omega2)*(cosd(theta3)+j*sind(theta3))

Aq=(j*r2*alpha2-r2*omega2*omega2)*(cosd(theta2)+j*sind(theta2))
Ap=(j*r2*alpha2-r2*omega2*omega2)*(cosd(theta2)+j*sind(theta2))
+(j*r3*alpha3-r3*omega3*omega3)*(cosd(theta3)+j*sind(theta3))
-------------------------------------------------------------------------------------------------
執行結果：

mode=1

mode=-1

8.2
使用drawlinks畫出四連桿形狀，再利用之前f4bar所得到的資料畫出速度線、加速度線：
-------------------------------------------------------------------------------------------------
r=[4 3 3 5]; td2=10; tdd2=0; theta1=0; theta2=45; linkdrive=0;
mode=1;
[data,form] = f4bar(r,theta1,theta2,td2,tdd2,mode,linkdrive)
r2=r(2);
r3=r(3);
theta2=data(2,2);
theta3=data(3,2);
omega2=data(2,3);
omega3=data(3,3);
alpha2=data(2,4);
alpha3=data(3,4);

% 速度、加速度公式
Vp=j*r2*(omega2)*(cosd(theta2)+j*sind(theta2));
Vq=j*r2*(omega2)*(cosd(theta2)+j*sind(theta2))
+j*r3*(omega2)*(cosd(theta3)+j*sind(theta3));
Ap=(j*r2*alpha2-r2*omega2*omega2)*(cosd(theta2)+j*sind(theta2));
Aq=(j*r2*alpha2-r2*omega2*omega2)*(cosd(theta2)+j*sind(theta2))
+(j*r3*alpha3-r3*omega3*omega3)*(cosd(theta3)+j*sind(theta3));


Vp=Vp/10 % 速度線過長 調成十分之一
Vq=Vq/10 % 速度線過長 調成十分之一
Ap=Ap/100 % 加速度線過長 調成十分之一
Aq=Aq/100 % 加速度線過長 調成十分之一

drawlinks(r,theta1,theta2,mode,linkdrive);

% 把起點設在PQ兩個旋轉結上，再把速度向量、加速度向量加上去
VpVector=[real(data(1,6)) imag(data(1,6)) ; real(data(1,6)+Vp) imag(data(1,6)+Vp)];
VqVector=[real(data(2,6)) imag(data(2,6)) ; real(data(2,6)+Vq) imag(data(2,6)+Vq)];
ApVector=[real(data(1,6)) imag(data(1,6)) ; real(data(1,6)+Ap) imag(data(1,6)+Ap)];
AqVector=[real(data(2,6)) imag(data(2,6)) ; real(data(2,6)+Aq) imag(data(2,6)+Aq)];


line([VpVector(1,1);VpVector(2,1)],[VpVector(1,2);VpVector(2,2)],'color','c'); %Vp 青色線
line([VqVector(1,1);VqVector(2,1)],[VqVector(1,2);VqVector(2,2)],'color','b'); %Vq 藍色線
line([ApVector(1,1);ApVector(2,1)],[ApVector(1,2);ApVector(2,2)],'color','g'); %Ap 綠色線（和黑色連桿重疊會變青）
line([AqVector(1,1);AqVector(2,1)],[AqVector(1,2);AqVector(2,2)],'color','r'); %Aq 紅色線
-------------------------------------------------------------------------------------------------
執行結果：

mode=1
mode=-1

8.3
使用fb_angle_limits畫出四連桿形狀，再利用之前drawlinks畫出限制角度時的四連桿形狀：
-------------------------------------------------------------------------------------------------
r=[4 3 3 5]; td2=10; tdd2=0; theta1=0; theta2=45; linkdrive=0;
mode=1;

% 求出限制角度
[Ang1, Ang2]=fb_angle_limits(r,theta1,linkdrive)

% 將限制角度代入drawlink程式中，畫出限制角度時的四連桿形狀
drawlinks(r,theta1,Ang1,mode,linkdrive);
drawlinks(r,theta1,Ang2,mode,linkdrive);
-------------------------------------------------------------------------------------------------
執行結果：

mode=1

mode=-1

在x軸上的圖，是Ang1=28.9550度的圖；在x軸下的圖，是Ang2=331.0450度的圖。這個四連桿，不論是分支型或閉合型，極限角度都一樣，故兩張圖一樣。
值得注意的是，這個configuration就是從分支型轉換到閉合型那一瞬間configuration。拿Ang1=28.9550度的來看：若此時Q點往QR連線的右半平面移動，則四連桿會變成閉合型；若此時Q點往QR連線左半平面移動，則四連桿會變成分支型。

8.4
利用for loop將不同theta2時的四連桿形狀畫出。
-------------------------------------------------------------------------------------------------
r=[4 3 3 5]; theta1=0; linkdrive=0;
mode=1;

for k=1:18;

theta2=20*k;
drawlinks(r,theta1,theta2,mode,linkdrive);

end
-------------------------------------------------------------------------------------------------
執行結果：

mode=1

mode=-1

不論mode=1或mode=-1，我們都可以發現在theta2=20/340/360處畫不出來。老師的drawlinks是設計成接收到無法實際存在的configuration時，會顯示出"Combination of links fail at degrees xxx.x"的字樣，告訴使用者該四連桿無法作xxx.x度的輸入。
原本在8.3小題我們只知道兩極限角度，但不知道可以移動的範圍是優弧(在此四連桿中為從28.9550度順時針繞到331.0450度)還是劣弧(在此四連桿中為從28.9550度逆時針繞回331.0450度)，但從8.4小題中，我們就可以判斷出這個四連桿可以移動的範圍是從28.9550度順時針繞到331.0450度。
另外，仔細觀察mode=1和mode=-1的兩張圖，可以發現兩張圖其實根本就是上下顛倒而已。這代表此四連桿有一個特性：若該連桿在x軸以上是分支型，則該連桿在x軸以下就會轉成閉合型，反之亦然。

8.5
先利用8.4小題的程式碼來判斷該連桿到底作怎麼樣的運動：
將8.4小題的程式碼中的theta2改成從極限角度30&330度代換，並畫圖觀察Q點的運動情形

-------------------------------------------------------------------------------------------------
r=[4 3 3 5]; theta1=0; linkdrive=0;
mode=1;

for k=1:11;

theta2=30*k;
drawlinks(r,theta1,theta2,mode,linkdrive);

end
-------------------------------------------------------------------------------------------------
執行結果：

mode=1

-------------------------------------------------------------------------------------------------

看mode=1的情況：
在theta逆時針接近331.0450度時，可看出Q會往左上方移動，故到達極限位置後，Q往左上方移動且P往右下方移動。這告訴我們這個四連桿在x軸上方為分支型，逆時針轉到x軸下方後變成閉合型，在通過331.0450度後又會變回分支型，再回到x軸上方後又變回閉合型，在通過28.9550度後又變回分支型......一直下去！mode=-1的情況就只是上下顛倒，不需要再討論。

所以在繪圖時：當驅動桿逆時針旋轉時，必須將mode設定為1；當驅動桿順時針轉回時，必須將mode設定-1。以上設定適用於起始位置在x軸上方且為分支型的情況，若上下顛倒或是型態對調都會使mode改變對調。

由於老師的drawlinks程式無法把axis鎖定在固定大小，所以我做了一點修改，命名為drawlinksNew：
-------------------------------------------------------------------------------------------------
因原程式碼貼上blog時會變成亂碼，故我以圖檔代替之：


基本上我只是在所有plot前面都加了axis([-5 5 -5 5])，把整張圖鎖在10*10的座標裡面。另外在最後面的axis equal也刪除，因為已經沒有必要。

以下是動畫的程式碼：
-------------------------------------------------------------------------------------------------
r=[4 3 3 5]; theta1=0; linkdrive=0;
mode=1;

axis([-5 5 -5 5]);

% k矩陣的row1：從k(1,1)到k(1,303)是儲存逆時針旋轉的角度值，從k(1,304)到k(1,606)是儲存逆時針旋轉的角度值。
% k矩陣的row2：從k(1,1)到k(1,303)是儲存逆時針旋轉的mode，從k(1,1)到k(1,303)是儲存逆時針旋轉的mode。

k=[linspace(29,331,303) linspace(331,29,303); linspace(1,1,303) linspace(-1,-1,303)];

while 1
for l=1:604;

clf;
theta2=k(1,l);
mode=k(2,l);
drawlinksNew(r,theta1,theta2,mode,linkdrive);
pause(0.01);

end
end
-------------------------------------------------------------------------------------------------
執行結果：由於我的電腦同時執行Matlab跟CamStudio會變慢，所以動畫有點不順，請見諒！

Grashof Inequality之證明

最近從書上看到了葛拉索公式的證明，覺得很有趣，就自己整理了一下打成檔案，放在網路上給大家參考。

我做的pdf檔可以從http://homepage.ntu.edu.tw/~b94202029/GrashofInequality.pdf下載。

Grashof Inequality

考慮如圖1的四連桿。若要完成一完整的旋轉，則該四連桿必須通過圖2和圖3的位置。令a為AB桿長、b為BC桿長、c為CD桿長、d為AD桿（接地桿）長。

先假設AB桿長小於接地桿桿長：

a＜d

運用三角不等式（兩邊之和大於第三邊）於圖2，可以得到

a+d＜b+c

b＜c+a+d

c＜b+a+d

再運用於圖3，可以得到

d-a＜b+c

b＜c+d-a

c＜b+d-a

再來將六個不等式之中較弱的論證去除：

1. 由於(e)式為真則保證(b)式為真，故(b)式可去除。把(e)式寫成

a+b＜c+d

2. 由於(f)式為真則保證(c)式為真，故(c)式可去除。把(c)式寫成

a+c＜b+d

3. 由於(a)式為真則保證(d)式為真，故(d)式可去除。

最後將(e’)與(a)相加、(f’)與(a)相加，可以得到

2a+b+d＜2c+b+d

2a+c+d＜2b+c+d

又因一開始的假設為a，故a為最短桿桿長。再回頭看(a)、(e’)、(f’)：三個不等式的左手邊皆為(a+另一桿桿長)，其中一式必定代表(最短桿桿長+最長桿桿長)< (另兩桿桿長相加)。故最後得到

s+l＜p+q

其中s=a為最短桿桿長、l為最長桿桿長、p和q為另兩桿桿長。

　　注意上述論證有AB桿長小於接地桿桿長的假設，且在此假設下所得到的最後結果顯示最短桿必接在接地桿上。若要將此公式推廣（最短桿可能為接地桿或連結桿），只要將觀察者改設成與AB桿或BC桿或CD桿一起移動（也就是將座標系改架在AB桿或BC桿或CD桿上）即可，且將此公式推廣的同時，亦可得知最短桿兩端的旋轉結必可作360度旋轉的訊息。