2007年6月11日 星期一

hw13

1.試設計一組複式齒輪,使其轉速比為125(請說明思考步驟及結果)。
ans:
先用最簡單的想法來做這題
題目要求齒輪比為125,是5的立方,所以只要串聯三個1:5的齒輪組就可達成目的

齒輪一 1  齒輪二 5  齒輪三  25  齒輪四 125  


如果減少一個齒輪呢?
√125=11.1803 
所以可以設計成

齒輪一 1  齒輪二 11.1803  齒輪三  125  
這樣的好處是可以省空間及成本,並且減少動力的損耗

2.請指出本學期中你自己最感得意的一次作業(請說明其原因,且該作業必須在自己的部落格內)。
老實說,機動學的每個作業都讓我很累很痛苦
學期剛開始完全摸不著頭緒,一半是因為對matlab程式使用不熟悉,一半是因為自己不夠用心去研讀,常常看著題目思考很久還不曉得如何下手,最令我頭痛的是,居然連請教同學都還聽不太懂!讓我感到信心全失!
後來跟著同學一同討論,上課也更加認真聽講,我漸漸進入狀況,之前還常常跟同學討論就是整夜沒睡的,也很謝謝大家肯提示我幫忙我,尤其時林詠舜同學張志鵬同學及蔡智鴻同學,非常謝謝他們!
我檢視自己的作業,覺得最滿意的就是作業六了
可能是因為作業六用到程式的部份比較少吧,當我看到題目時就比較不惶恐
自己細心的把那部分課本及講義讀通,然後開始計算
求出答案來時我非常有成就感!

還有一次作業令我印象深刻,作業七
這份作業我相當用心的跟同學研究程式碼
花了很多很多時間
剛開始有一小部份圖畫出來,高興的要命!讓我更有動力繼續下去了
後來我們遇到瓶頸,試過很多方法都不曉得如何解決,當時已經是凌晨兩點了!
於是我們就先關電腦休息,隔天再去學校請教已經完成的高手
最後當然是有順利完成該項作業阿,真的很謝謝同學的幫助

當我越來越了解老師循序漸進的教導時
便抓到作業操作的感覺了
我本來就不是特別聰明的學生,那時我已經落後大家很多很多了
但是慢慢清楚老師撰寫的程式用法
做作業也變得輕鬆許多,也不用那麼常問東問西的,同學說個大概,我就可以了解
比起當初,感覺真的很好
雖然我並沒有把程式摸的熟透,但是這堂課讓我學到很多
matlab的一些程式我都比較會運用了

謝謝老師!

2007年6月1日 星期五

hw12

1.5/31日全程來上課。
2.一組標準全齒輪齒輪之徑節為8,齒數分別為30T與48T,其工作壓力角為25度試求其接觸線長度,與接觸比。
兩齒輪之節圓、基圓直徑各為如何?請列式計算其結果。
此組齒輪是否會產生干涉現象?試列式證明之。
可否利用draw_gear.m繪出其接合情形,並繪出其動畫效果。

接觸比: (Contact ratio),表在同一時間內有多少對牙齒同時相互嚙合。
接觸比mc應為接觸路徑長度與基周節之比。可以用MATLAB函式
contact_ratio()計算
其中輸入參數:

Pd:徑節 n2, n3:兩齒輪之齒數 phi:壓力角function [c_ratio,c_length,ad,pc,pb,d2,d3,ag]=contact_ratio(pd,n2,n3, phi)
d2g=pi/180;
pangle=phi*d2g;
cosx=cos(pangle);sinx=sin(pangle);
ad=1./pd;pc=pi./pd;
pb=pc.*cosx;
r2=n2./(2*pd);r3=n3./(2*pd);d2=2*r2;d3=2*r3;
rb2=r2.*cosx;rb3=r3.*cosx;
ax=sqrt((r3+ad).^2-(r3.*cosx).^2)-r3.*sinx;
xb=sqrt((r2+ad).^2-(r2.*cosx).^2)-r2.*sinx;
c_length=ax+xb;
c_ratio=c_length./pb;
ag1=[ax./rb2 xb./rb2 c_length./rb2]/d2g;
ag2=[ax./rb3 xb./rb3 c_length./rb3]/d2g;
ag=[ag1;ag2];

輸入:
[c_ratio, c_length,ad,pc,pb,r2,r3,ag]=contact_ratio(8,30,48,25)
得到下列資訊:
c_ratio(接觸比) =
1.5028(無單位)
c_length (接觸長度)=
0.5349(吋)
ad (齒冠)=
0.1250(吋)
Pc(周節) =
0.3927
pb (基周節)=
0.3559
r2(齒輪節圓直徑) =
3.7500(吋)
r3 (齒輪節圓直徑)=
6(吋)
ag =(度)
(齒輪2之接近角) 9.1921 (遠退角)8.8419 (作用角)18.0340
(齒輪3之接近角)5.7450 (遠退角)5.5262 (作用角)11.2712
Ans:
接觸長度: 0.5349吋
接 觸 比: 1.5028
節圓直徑: 節圓直徑=齒數÷徑節 D=N/Pd
3.75吋 / 6吋
基圓直徑: 基圓直徑=節圓直徑×cos(壓力角) D'=D×cos(ψ)
3.75cos25=3.398吋
6cos25 = 5.438吋
根據公式9.47,測試干涉之條件為: (N2²+2N2 x N3)sin²>= 4 + 4N3
(30²+2*30*48)sin²25>=4(1+48)
可知不會干涉
亦可使用function [x]=isinterf(phi,N1,N2)
輸入: isinterf(25,30,48)
輸出: ans = 0
可知不會干涉



draw_gear.m
輸入參數分別定義如下:Dp: 節矩N: 齒數phi: 壓力角range: 繪出之部份x0,y0: 齒輪中心座標

輸入move2_gear(10 ,30,48,25,10)
輸入項分別是節矩,兩齒輪齒數,壓力角,齒輪一的角速度
得到動畫

http://www.youtube.com/watch?v=zFhtuYXliho

2007年5月25日 星期五

hw11

本人范詠晴(5/24)有來上課

1.某凸輪開始時先在0-100∘區間滯留,然後提升後在200至260∘區間滯留,其高度(衝程)為 5公分,其餘l由260∘至360∘則為返程。升程採用等加速度運動,返程之運動型式自定。設刻度區間為10∘,試繪出其高度、速度及加速度與凸輪迴轉角度間之關係。

ans:
由題目知道凸輪在0-100度之間是圓形,沒有凸出,因此沒有衝程
而當100度提升到200度時,此段有提升,非圓形因此有衝程
程式輸入升程的地方要輸入100度到200度衝程是5公分
返程也就是此段下降,輸入260度
題目又要求升程是等加速度上升 可用拋物線的形式,代號2

使用函式function plot_dwell(ctheta,s,pattern,range)來繪出高度、速度及加速度與凸輪迴轉角度間之關係圖

輸入plot_dwell(0:10:360,5,[2 1],[100 200 260])
得到下圖



2.設凸輪之半徑為15公分,以順時針方向旋轉,其從動件為梢型,垂直接觸,長為10公分,從動件之運動係依照第二項之運動型式。試繪出此凸輪之工作曲線。
ans:
本題使用函式function [x,y]=pincam(cth,r0,s,e,L,range,pattern,cw)
參數解釋:
cth 為上題的ctheta 代入[0:10:360]
r0 為桿長15
s 是衝程5
e 偏置量0(若升程和返程都是用等加速度運動,就沒有偏置量)
L 桿長10
range 使用[100 200 260]
pattern 代入[2 1]
cw 是凸輪轉動方向,逆時鐘為1,順時鐘為-1
輸入[x,y]=pincam([0:10:360],15,5,0,10,[100 200 260],[2 2],-1)


3.你能讓此凸輪迴轉嗎?
ans:
可以
參考講義上的pincam程式再做修改
function [x,y]=pincam_2(cth,r0,s,e,L,range,pattern,cw)
figure(1);
n=0;
for mm=1:10:360;
n=n+1;
m=mm*pi/180;
clf;
th=cth*pi/180;
s0=sqrt(r0*r0-e*e);f
or i=1:length(cth)
t=th(i)*cw;
A=[cos(t-m) -sin(t-m);sin(t-m) cos(t-m)];
[ym,yy,yyy]=dwell(cth(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);
end
[yw,yyw,yyyw]=dwell(cth,range,pattern)
y1=yw*s+r0;
y2=yw*s+r0+L;
line([0 0],[y1(n) y2(n)],'linewidth',3,'color','red')
hold on;
plot([0 x],[0 y],'ro',x,y,'k-')
axis ([-70 70 -70 70])
pause(0.03)
end

動畫 http://www.youtube.com/watch?v=PxB-ncghhrM

2007年5月24日 星期四

Meeting紀錄與心得

針對5/17早上的與教授小組討論
我歸納了幾項在機動學作業上可再改進的部份


作業重點
在做作業時,最重要的是一個程式的輸入輸出參數及其單位
再來是解釋程式碼還有使用的方法介紹呼叫的目的,這樣才有確實學到東西,才能真的瞭解程式的撰寫及執行
最後,老師表示其實有沒有放function是無所謂的,因為同學所使用的都是修改自老師的教學講義提供的function
但我認為可以附註在每次作業的最後面,以便未來要參考時馬上就有function可以研究
更可以提供其他有興趣的網友這些function資料

blog編輯
不少人覺得部落格的編輯方法有些地方不順手
1.貼上圖片後,文字行距需重新調整
解決辦法是:先貼圖,再放文字,如此一來可省下很多時間

2.在kkman裡頭開啟,無法跳行,也就是編輯文章時按enter沒有用
解決方法是:用ie開啟來寫文章

3.調整文章排版,有時設定好還是會跑掉,要花兩倍時間整理,相當不便
解決方法:在文章前加上<pre>
以及在文章後加上</pre>
謝謝老師提供此法

4.頁面色彩
盡量挑乾淨簡單的樣式,也不要選用亮黃色等讓人看了不舒適的顏色

2007年5月22日 星期二

hw10

本週(5/17)有來上課
請思考速度與加速度的問題,當一桿以某特定點M等角速度迴轉時,其端點P之速度方向如何?其加速度方向如何?若該特定點M復以等速水平運動,則同一端點P之速度與加速度方向會變為如何?若M點同時也有加速度,則點P會有何變化?若以此推理四連桿的運動,則點P與Q之速度與加速度方向會與桿一(固定桿)之兩端點之關係如何?與我們前面的作業分析結果有無共通之處?(參看第六章之四連桿機構之運動分析)
設有一運動之曲柄滑塊連桿組合,設滑塊之偏置量為零,且在水平方向移動,試以此機構之曲桿長度及角度,以及連結桿之長度為輸入項,利用matlab寫出一程式計算在不同曲柄角度時,六點瞬心之對應位置。可順便探討六點瞬心與曲柄角間之關係。


ans:分析
一.
設M為圓心,PM為其半徑r(m),具角速度ω(rad/s)。
從動力學中公式v=ωxr可以知道,此桿所得到之速度為ωr。再由右手定則得知其速度之方向為ω指向r之大拇指方向(即垂直r)
因為是等角速度轉動,角加速度為0
但有向心加速度r*ω^2 (m/s^2)

二.
若以點M以等速V(m/s)水平移動,可由向量式表現其位置:PM桿向量為s = rcos(θ)i+rsin(θ)j,
則此時速度向量可以表示為V(m/s)= r*ω*sin(θ)(m/s) i - r*ω*cos(θ)(m/s) j,加速度同上題

三.
若此時M點也具有加速度a(m/s^2)時,則M點所具有的速度應為原本的再加上瞬時加速度所造成的速度,
其應為v = r*ω*sin(θ)(m/s) i - r*ω*cos(θ)(m/s) j + a對θ之積分
其加速度為V(m/s)+r*ω^2*cos(θ)(m/s) i + r*ω^2*sin(θ)(m/s) j

四.
綜合前面結果,P點與Q點之速度與加速度方向將與固定桿之兩端點關係呈現出一個特定個關係,
以第一桿為固定桿,P點與第二桿作連結,若將P點限制在以第二桿為半徑的圓周上,
即為將P點與Q點之速度與加速度以及M點為圓心、PQ兩點為連結之分析過程一樣



程式:
function slider_draw(R,L,e)
%R=桿一長
%L=桿二長
%e=偏置量
th1=0; %R與L在同一平面上
if R>L
th2=asind(L/R);
else
th2=90;
end
th=linspace(th1,th2,100);
d=slider_solve(th,R,L,e,1);
x=R*cosd(th);
y=R*sind(th);
for n=1:100
hold on;
%畫出桿子與滑塊
line([0,x(n),d(n)],[0,y(n),e]);
line([d(n)-4,d(n)+4,d(n)+4,d(n)-4,d(n)-4],[e-3,e-3,e+3,e+3,e-3]);
%畫出瞬心位置
plot(0,0,'r*');
plot(x(n),y(n),'r*');
plot(d(n),e,'r*');
plot([0,0],[0,e-d(n)*(y(n)-e)/(x(n)-d(n))],'r*:');
plot([x(n),0],[y(n),e-d(n)*(y(n)-e)/(x(n)-d(n))],'r*:');
plot([x(n),d(n)],[y(n),y(n)*d(n)/x(n)],'r*:');
plot([d(n),d(n)],[0,y(n)*d(n)/x(n)],'r*:');
axis equal;
axis ([-70 70 -70 70]);
pause(0.05);
clf;
end

動畫
http://www.youtube.com/watch?v=N1c2Ls_eNGI

2007年5月6日 星期日

hw9

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

R=10+43=53; L=53+5=58曲桿 R=53cm連桿 L=58cm偏置量 e=10cm滑塊機構的位置分析:滑塊能移動,d的長度會變利用課本上提供的slider_solve程式,算出在幾種曲柄角度之下之對應d的長度和聯結桿角度

用slider_limit程式
[s,theta1,theta2]=slider_limit(53,58,10)
得到
s = 1.1757e+002 -8.6603e+000itheta1 = 4.9031theta2 = 90.0000 -75.4561i
(s為衝程,theta1和theta2為界限角)
>>abs(s)
衝程
s=117.8932 (cm)
分析由曲桿驅動和滑塊驅動的運動狀況:
動畫網址
http://www.youtube.com/watch?v=2_gKqNVY-Jw
@曲桿驅動畫四連桿,輸入
drawldlinks([75,53,58,10],0,60,1,0)
使用程式sld_angle_limits、sldbox及drawsldlimits計算並畫運動的界限驅動桿為第二桿(driver=0),又符合講義上第三項的條件"r3>-r4,且r2>r3-r4",sigma=1
>>drawsldlimits([75,53,58,10],0,1,0)
輸出:
Qstart =3.6000e-006Qstop =360.0000
四連桿若用曲桿驅動,可從0.0000036度到360度(曲桿與水平的夾角)轉一圈。
@滑塊驅動畫出四連桿的位置
>>drawldlinks([75,53,58,10],0,60,1,2)

再畫出其運動的界限:此四連桿符合|r2-r3|<r4 的條件,sigma=1;由滑塊驅動,driver=2得到上述結果後,便可輸入
drawldlimits([75,53,58,10],0,1,2)
得:
Qstart =-116.5719Qstop =116.5719

當滑塊驅動時,限制是從-116.5719度到116.5719度



hw8

本人(4/26)有上課
有一組四連桿,其桿長分別為r=[4 3 3 5],由桿2驅動,設第一固定桿角度theta1=0度;
角速度 td2=10rad/s; 角加速度tdd2=0 rad/s^2。

問題一:設桿2角度theta2=45度時,求各點之位置、速度與加速度為何?

角度:
>> abs(val(:,2))'

ans =

0 45.0000 166.4283 163.5328

角速度:
>> abs(val(:,3))'

ans =

0 10.0000 16.2681 4.9677
角加速度
>> abs(val(:,4))'

ans =

0 0 491.4428 383.6120
問題二:繪出此四連桿之相關位置及標明各點之速度方向及大小(以程式為之)。






問題三:當桿2迴轉時,求出此組四連桿之限制角度,並繪出其位置(以程式為之)。
>> [Ang1, Ang2]=fb_angle_limits([4 3 3 5],0,0)

Ang1 =

28.9550


Ang2 =

331.0450



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


問題四:設theta2=[0:20:360],試繪出此組四連桿之重疊影像,解釋為何有些沒有值。

>> for i=10:20:360,drawlinks([4 3 3 5],0,i,-1,0); end
Combination of links fail at degrees 10.0
Combination of links fail at degrees 350.0

因為由程式求出限制角度範圍,所以會有兩項沒有值

問題五:若將問題三考慮在內,只在可迴轉的範圍內迴轉,請問你能讓此組四連桿作成動畫方式迴轉嗎?

for i=10:20:360,drawlinks([4 3 3 5],0,i,-1,0);
pause(0.5)
end
動畫網址:http://www.youtube.com/watch?v=S0ffsSi0lu4

本作業需用到的function
function [data,form] = f4bar(r,theta1,theta2,td2,tdd2,mode,linkdrive)
%function [data,form] = f4bar(r,theta1,theta2,td2,tdd2,mode,linkdrive)
% This function analyzes a four-bar linkage when the crank is the
% driving link. The input data are:
% theta1,theta2 are angles in degrees
% r=[r1 r2 r3 r4]= lengths of links(1=frame)
%td2 = crank or coupler angular velocity (rad/sec)
%tdd2 = crank or coupler angular acceleration (rad/sec^2)
%mode = +1 or -1. Identifies assembly mode
%linkdrive = 0 for crank as driver; 1 for coupler as driver
%data (1:4,1) = link positions for 4 links
%data (1:4,2) = link angles in degrees
%data (1:4,3) = link angular velocities
%data (1:4,4) = link angular accelerations
%data (1,5) = velocity of point Q
%data (2,5) = velocity of point P
%data (3,5) = acceleration of point Q
%data (4,5) = acceleration of point P
%data (1,6) = position of Q
%data (2,6) = position of P
%form = assembly status. form = 0, mechanism fails to assemble.
if nargin<7,linkdrive=0;end mode="1;end" data="zeros(4,6);" linkdrive="=" r="[r(1)" rr="r.*r;d2g=" t1="theta(1);tx=" s1="sin(t1);c1=" sx="sin(tx);cx=" 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)*s1-2*r(2)*r(4)*sx;" pos="B*B-C*C+A*A;">=0,
form=1;
% Check for the denominator equal to zero
if abs(C-A)>=1e-5
t4=2*atan((-B+mode*sqrt(pos))/(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);
else
% 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);
pos=B*B-C*C+A*A;
if pos>=0,
t3=2*atan((-B-mode*sqrt(pos))/(C-A));
s3=sin(t3); c3=cos(t3);
t4=atan2((-r(1)*s1+r(3)*s3+r(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)*c4;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 data
% coordinates of P and Q
if linkdrive==1,
c2=c3;c3=cx;s2=s3;s3=sx;
r(2:3)=[r(3) r(2)];theta(2:3)=[theta(3) theta(2)];
td(2:3)=[td(3) td(2)];tdd(2:3)=[tdd(3) tdd(2)];
else
c2=cx;s2=sx;
end
for j=1:4,
data(j,1:4)=[r(j)*exp(i*theta(j)) theta(j)/d2g td(j) tdd(j)] ;
end % position vectors
data(1,5)=r(2)*td(2)*exp(i*theta(2));%velocity for point Q
data(2,5)=r(4)*td(4)*exp(i*theta(4));%velocity for point P
data(3,5)=r(2)*(i*tdd(2)-td(2)*td(2))*exp(i*theta(2));%acc of Q
data(4,5)=r(4)*(i*tdd(4)-td(4)*td(4))*exp(i*theta(4));%acc of P
data(1,6)=data(2,1);%position of Q, again
data(2,6)=data(1,1)+data(4,1);% position of P

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


須輸入的參數說明:

r(1:4) = 各桿之長度,r(1)為固定桿,其餘分別為曲桿、結合桿及被動桿,即r=[r1 r2 r3 r4]。
theta1 = 第一桿之水平角,或為四連桿之架構角,以角度表示。
theta2 = 驅動桿之水平夾角,以角度表示。一般為曲桿角,但若為結合桿驅動,則為結合桿之水平夾角。
td2 = 驅動桿(第二桿或第三桿)之角速度(rad/sec)。
tdd2 = 驅動桿(第二桿或第三桿)之角加速度(rad/sec^2)。
mode = +1 or -1. 組合模數,負值表示閉合型,正值為分支型,但有時需視實際情況而定。
linkdrive = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿) 。




畫出連桿
function [values]=drawlinks(r,th1,th2,mode,linkdrive)
if nargin<5,linkdrive=0;end mode="1;end" rr="values(:,1);" rx="real(rr(:,1));rx(4)=" ry="imag(rr(:,1));ry(4)=" b="=" linkdrive="=">

function [Ang1, Ang2]=fb_angle_limits(r,theta1,linkdrive)
if nargin<3,linkdrive=0;end theta1="0;end" linkdrive="=" r="[r(1)" r1="r(1);r2=" r3="r(3);r4=" rmin="min(r);rmax=" rtotal="sum(r);" ang1="0;" ang2="2*pi;">(r3+r4)& abs(r1-r2)(r3+r4)& abs(r1-r2)>=abs(r3-r4)
Ang1=acos((r2^2-(r4+r3)^2+r1^2)/(2*r1*r2));
Ang2=-Ang1;
end
if (r1+r2)<=(r3+r4)&abs(r1-r2)< ang1="acos((r2^2-(r4-r3)^2+r1^2)/(2*r1*r2));" ang2="2*pi-Ang1;" adjst =" (Ang2-Ang1)*1e-8;" ang1="theta1+(Ang1+adjst)*180/pi;" ang2="theta1+(Ang2-adjst)*180/pi;" tt="ang1;ang1=" ang2="tt;end

function drawlimits(r,th1,sigma,driver)
[Qstart, Qstop]=fb_angle_limits(r,th1,driver)
clf;
[values b]=f4bar(r,th1,Qstart,0,0,sigma,driver);
rr=values(:,1);rr(3)=rr(1)+rr(4);
rx=real(rr);rx(4)=0;
ry=imag(rr);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');
text(rx(2)/2,ry(2)/2,['s1=' num2str(Qstart,'%6.1f')]);
else
fprintf('Combination of links fails at degrees %6.1f\n',Qstart);
end
[values b]=f4bar(r,th1,Qstop,0,0,sigma,driver);rr=values(:,1);
rr(3)=rr(1)+rr(4);
rx=real(rr);rx(4)=0;
ry=imag(rr);ry(4)=0;
if b==1
if driver==0
plot([0 rx(2)],[0 ry(2)],'b-','LineWidth',1);
plot([rx(2) rx(3)],[ry(2) ry(3)],'r-','LineWidth',1.5);
else
plot([0 rx(2)],[0 ry(2)],'r-','LineWidth',1.5);
plot([rx(2) rx(3)],[ry(2) ry(3)],'b-','LineWidth',1);
end
plot([rx(1) rx(3)],[ry(1) ry(3)],'g-');
plot(rx,ry,'bo');
text(rx(2),ry(2),' p');text(rx(3),ry(3),' q');
text(rx(2)/2,ry(2)/2,[' s2=' num2str(Qstop,'%6.1f')]);
else
fprintf('Combination of links fail at degrees %6.1f\n',Qstop);
end
axis equal
grid on