畢業(yè)設(shè)計論文基于MATLAB的SMITH圓圖的設(shè)計

1、前 言smith圓圖是p.h.smith于1939年在貝爾實(shí)驗(yàn)室發(fā)明的，它主要用于計算微波網(wǎng)絡(luò)的阻抗、導(dǎo)納及網(wǎng)絡(luò)阻抗匹配設(shè)計，還可用于設(shè)計微波元器件。smith圓圖軟件不僅適用于微波工程設(shè)計，亦可用于電磁場、微波技術(shù)及天線與電波傳播等課程相關(guān)內(nèi)容的教學(xué)。 多年來，smith圓圖作為一個不可缺少的設(shè)計工具，在高頻和微波領(lǐng)域發(fā)揮了重要作用。掌握了smith圓圖的原理和特性，并熟悉其阻抗匹配技術(shù)，將使設(shè)計人員的設(shè)計能力大為加強(qiáng)。盡管目前微波電路cad技術(shù)已進(jìn)入了設(shè)計領(lǐng)域，人們?nèi)源罅渴褂脠A圖作為設(shè)計工具。究其原因，一是因?yàn)楝F(xiàn)有的大部分cad軟件在設(shè)計時要求用戶提供初始電路（包括拓?fù)浜驮担?，因此仍?/p>

2、借助于smith圓圖進(jìn)行圖解設(shè)計；另一方面，對于窄帶情形，利用圓圖通?？芍苯油瓿稍O(shè)計不 必進(jìn)一步進(jìn)行優(yōu)化而且目前不少大型軟件價格昂貴，cad技術(shù)尚未普及，因而一款smith圓圖工具的設(shè)計仍然是非常有必要的。微波網(wǎng)絡(luò)的正弦穩(wěn)態(tài)分析含有復(fù)數(shù)計算，運(yùn)算十分繁瑣和耗時。在計算機(jī)運(yùn)算速度和內(nèi)存不夠發(fā)達(dá)以前，圖解分析法得到長足發(fā)展，其中多年來應(yīng)用最廣的是smith圓圖。在計算微波傳輸線輸入阻抗、導(dǎo)納及阻抗匹配等問題時，它不僅能避開繁瑣的公式及復(fù)數(shù)運(yùn)算，使工程設(shè)計中相關(guān)計算簡單便捷，而且圖解過程物理概念清晰，所得結(jié)果直觀形象。隨著計算機(jī)技術(shù)的飛速發(fā)展，圖解法在計算精度上的固有缺陷日益顯現(xiàn)，因?yàn)?，圓圖的計算精

3、度取決于圓圖中必須有足夠的圓周數(shù)，而過多的圓周會導(dǎo)致圖線過于密集，不便將阻抗、反射系數(shù)、行波系數(shù)及電長度等相關(guān)數(shù)據(jù)從圖上直接讀出。通過對圓圖構(gòu)成的基本原理和應(yīng)用問題的分析，利用現(xiàn)代計算機(jī)技術(shù)可以解決圓圖計算精度等方面存在的問題，為此設(shè)計的smith圓圖軟件既保留圓圖計算直觀、便捷的大眾性，又滿足工程設(shè)計中相關(guān)參數(shù)的計算精度。在計算機(jī)應(yīng)用日益普及的今天，該軟件特別適合電磁場、微波技術(shù)與天線等領(lǐng)域的教學(xué)和工程設(shè)計相關(guān)參數(shù)計算使用。 目 錄摘 要4第章 概述 5第章 smith圓圖軟件構(gòu)成的基本原理 72.1 阻抗圓圖軟件構(gòu)成的基本原理72.3 圓圖軟件的設(shè)計特點(diǎn) 82.3 圓圖的基本應(yīng)用10第章

4、仿真與調(diào)試113.1 調(diào)試113.2 仿真123.3 仿真結(jié)果驗(yàn)算13第章 設(shè)計小結(jié) 14參考文獻(xiàn) 15附錄 源程序代碼17基于matlab語言的smith圓圖軟件設(shè)計摘要：smith圓圖是從事微波工程實(shí)驗(yàn)和天線設(shè)計的重要工具。應(yīng)用matlab作為開發(fā)平臺研制的smith圓圖應(yīng)用軟包，使圓圖的應(yīng)用和計算變得更加方便、快捷，該軟件具有用戶圖形界面，簡單易用，而且計算精度高。關(guān)鍵詞： matlab；gui；smith圓圖； 第1章 概述smith圓圖的基本在于以下的算式: 式中的代表其線路的反射系數(shù)，即阻抗匹配時, s矩陣?yán)锏膕11，zl是歸一化負(fù)載抗，即zl / z0。其中，zl是電路的負(fù)載值z

5、0是傳輸線的特性阻抗值，通常會使用50。圖1-1中的圓形線代表電阻抗力的實(shí)數(shù)值，即電阻值，中間的橫線與向上和向下散出的線則代表電阻抗力的虛數(shù)值，即由電容或電感在高頻下所產(chǎn)生的阻力，當(dāng)中向上的是正數(shù)，向下的是負(fù)數(shù)。圓圖圖最中間的點(diǎn)(1+j0)代表一個已匹配的電阻數(shù)值(zl)，同時其反射系數(shù)的值會是零。圖表的邊緣代表其反射系數(shù)的長度是1，即100%反射。在圖邊的數(shù)字代表反射系數(shù)的角度(0-180度)和波長（由零至半個波長）。有一些圖表是以導(dǎo)納值來表示，把上述的阻抗值旋轉(zhuǎn)180度即可。自從有了計算機(jī)后，此種圖表的使用率隨之而下，但仍常用來表示特定的資料。對于就讀電磁學(xué)及微波電子學(xué)的學(xué)生來說，在解決課

6、本問題仍然很實(shí)用，因此史密斯圖至今仍是重要的教學(xué)用具。 本設(shè)計是運(yùn)用matlab編寫smith圓圖仿真程序，整個圓圖軟件分為用戶圖形界面模塊、圓圖計算模塊、畫圖演示模塊。圓圖計算模塊分為輸入阻抗計算、反射系數(shù)計算、行波系數(shù)、駐波比計算以及整個smith圓圖；畫圖演示模塊分為等歸一化電阻圓、等歸一化電抗圓、反射系數(shù)圓；確定阻抗值在圓圖上的位置、圓圖的基本應(yīng)用、求輸入阻抗及其在圓圖上的位置。具體實(shí)現(xiàn)：1、 本程序能讀出圓圖上的任意一點(diǎn)對應(yīng)的各個值，并能夠根據(jù)輸入的歸一化阻抗畫出相應(yīng)的圓圖、顯示對應(yīng)的值；2、 當(dāng)在matlab環(huán)境下運(yùn)行程序后，會顯示完整的圓圖工具界面，在整個界面的右半部分畫出完整的

7、smith圓圖；3、 當(dāng)在“輸入數(shù)據(jù)”文本框里分別輸入負(fù)載阻抗的實(shí)部及虛部，點(diǎn)擊“確定”按鈕，會在繪圖區(qū)域內(nèi)畫出圓圖，并顯示相應(yīng)的阻抗、反射系數(shù)、駐波比、行波系數(shù)的值。第2章 smith圓圖軟件構(gòu)成的基本原理2.1 阻抗圓圖軟件構(gòu)成的基本原理 圓圖運(yùn)算的基礎(chǔ)是反射系數(shù) ( =u+jv ) 。smith圓圖由反射系數(shù)平面上的等歸一化反射系數(shù)圓族、等歸一化電阻圓族、等歸一化電抗圓族構(gòu)成。利用已知的歸一化阻抗zl=r+jx，實(shí)現(xiàn)等反射系數(shù)圓族程序如下： r=r；x=x；%輸入歸一化阻抗 u= ( r2+x2-1) / ( r2+23r+1+x2) ； v= ( 23x) / ( r2+23r+1+x

8、2) ； tr=23pi3 ( 00.011) ； =sqrt (u2+v2) ；%等反射系數(shù)圓的半徑plot ( r3cos ( tr) ，r3sin ( tr) ，y ) %畫等反射系數(shù)圓實(shí)現(xiàn)等歸一化電阻圓族程序如下：for rr=1/ ( 1+r) ；cr=1-rr；%畫電阻圓 plot ( cr+rr3cos ( tr) ，rr3sin ( tr) ，m )實(shí)現(xiàn)等歸一化電抗圓族程序如下：forx=x；%畫電抗圓 rx=1/x；cx=rx； iftx<pi； plot ( 1-rx3sin ( tx) ，cx-rx3cos ( tx) ，m)else上述作圖在u為橫坐標(biāo)和v為縱坐標(biāo)

9、的復(fù)平面中，r=常數(shù)的曲線是一系列的電阻圓，x=常數(shù)的曲線為一系列的電抗圓。需要指出，圓圖是用極坐標(biāo)來表示反射系數(shù)的。因?yàn)閨=常數(shù)的等反射系數(shù)圓是以原點(diǎn) ( u=0，v=0) 為圓心，以 |為半徑，而相角=常數(shù)的曲線則為經(jīng)過坐標(biāo)原點(diǎn)的一系列直線。為避免混淆，一般只畫出由已知的歸一化阻抗zl=r+jx對應(yīng)的代表大小的同心圓和代表相角的直線，但計算所得的電長度等參數(shù)以菜單形式給出。由上述程序，根據(jù)需要容易畫出相應(yīng)圓圖。 圖2-1 smith圓圖軟件主界面2.2 圓圖軟件的設(shè)計特點(diǎn) 圓圖軟件設(shè)計要求計算結(jié)果以圖形和數(shù)據(jù)并行輸出，處理包括復(fù)數(shù)的矩陣運(yùn)算。為使程序代碼簡單，執(zhí)行運(yùn)算速度快，計算精度高，選

10、擇matlab軟件作為設(shè)計技術(shù)平臺較為理想。該軟件數(shù)學(xué)表示、函數(shù)集豐富且功能強(qiáng)大、良好的用戶界面以及許多函數(shù)本身會繪制圖形且自動選取坐標(biāo)刻度等顯著優(yōu)點(diǎn)，特別適用于大量計算，因此smith圓圖軟件選擇matlab語言來編寫。在matlab中句柄圖形有一定的優(yōu)先級，我們設(shè)計的程序也是按照這樣的順序進(jìn)行的，如圖2-2所示。在這個框架中，最高級的是主界面，它是下面所有界面句柄的父。然后，從主界面上連接到下一層界面中，上一級為下一級的父，每級間通過回調(diào)來實(shí)現(xiàn)連接和數(shù)據(jù)傳遞。圖2-2圓圖軟件程序結(jié)構(gòu)框圖 matalab將圖形窗口 ( figure ) 、標(biāo)題 ( title ) 、文本( text ) 、

11、坐標(biāo)軸 ( axes) 、菜單 ( menu) 、控制框 ( uicontrol )等均作為對象處理，并且有豐富的對象屬性，只要利用get，set，reset，delet，gco，findobj等函數(shù)即可方便地操作這些圖形對象。軟件主頁圖形界面由圓圖和主菜單按鈕組成。界面比較簡潔，窗口中僅有固定的圖形和菜單，按任意菜單鍵即可進(jìn)入相關(guān)計算。菜單的設(shè)計很簡單，只要利用figure，title和text，設(shè)置相應(yīng)的屬性即可。為便于使用，主菜單均使用按鈕，按下按鈕即彈出相應(yīng)的輸入框，主要用到控制框?qū)ο笾械陌粹o類型。輸入框格式與一般windows中的輸入框格式基本相同，上面為各輸入文本編輯條，下邊有“確

12、定”按鈕，只是按鈕相聯(lián)系的功能函數(shù)不同，輸入內(nèi)容也不同。整個輸入框是一個特殊的圖形窗口，其中輸入編輯條由unicontrol對象中的可編輯文本框 ( edit ) 類型實(shí)現(xiàn)。另外，這里的“確認(rèn)”按鈕，由于回調(diào)函數(shù)需要獲得輸入文本的內(nèi)容，實(shí)現(xiàn)比較復(fù)雜， 要利用執(zhí)行字符串命令 ( eval ) 或執(zhí)行字符串函數(shù) ( feval )。 基于上述設(shè)計目標(biāo)，利用matlab強(qiáng)大的作圖功能容易畫出完整的smith圓圖。整個圓圖軟件分為用戶圖形界面模塊、圓圖計算模塊、畫圖演示模塊。上述3大模塊又進(jìn)一步分解，其中用戶圖形界面模塊分為：主頁、主菜單；圓圖計算模塊分為反射系數(shù)計算、單支節(jié)匹配計算、輸入阻抗計算以及

13、整個smith圓圖；畫圖演示模塊分為等歸一化電阻圓、等歸一化電抗圓、反射系數(shù)圓等；確定阻抗值在圓圖上的位置、圓圖的基本應(yīng)用、求輸入阻抗及其在圓圖上的位置以及單支節(jié)匹配等問題。2.3 圓圖的基本應(yīng)用 圓圖是從事天線和微波工程設(shè)計的重要計算工具之一。應(yīng)用圓圖軟件進(jìn)行工程計算，既保留圓圖簡便直觀、物理概念清晰等圖解法固有的優(yōu)點(diǎn)，又能保證設(shè)計計算必要的精確度，主要用于以下幾類計算： 由歸一化阻抗求電壓反射系數(shù)和電壓駐波系數(shù)以及相應(yīng)的電長度等，這類問題可由反射系數(shù)模塊解決； 單支節(jié)匹配等微波網(wǎng)絡(luò)匹配設(shè)計及微波元器件設(shè)計計算等由單支節(jié)匹配模塊實(shí)現(xiàn)； 由負(fù)載阻抗根據(jù)移動的電長度計算輸入阻抗、導(dǎo)納或逆運(yùn)算，計

14、算微波傳輸線上umax 和umin 的位置等，這些問題可由輸入阻抗計算模塊解決。第3章 仿真與調(diào)試3.1 調(diào)試（1） 運(yùn)行matlab軟件，在該環(huán)境下新建編輯器文件，輸入要調(diào)試的smith圓圖工具實(shí)現(xiàn)的程序；（2） 保存好該文件之后，點(diǎn)擊“debug”菜單中的“run”命令對程序進(jìn)行調(diào)試，如果提示有錯誤提示，改正之后重復(fù)調(diào)試，直至沒有錯誤為止；（3） 改正所有的錯誤之后，運(yùn)行，得到預(yù)期設(shè)計的smith圓圖軟件界面，如圖3-1所示：圖3-1 調(diào)試結(jié)果3.2 仿真（1）在用戶界面的“輸入數(shù)據(jù)”欄里分別輸入負(fù)載阻抗的實(shí)部與虛部（仿真時輸入實(shí)部為50，虛部為0），如圖3-2所示：圖3-2 輸入負(fù)載阻抗

15、（2）點(diǎn)擊“確定按鈕”，界面上顯示輸出數(shù)據(jù)：負(fù)載導(dǎo)納、駐波比、反射系數(shù)和行波系數(shù)的數(shù)值，同時在繪圖區(qū)域畫出其歸一化阻抗圓圖。分別如圖3-3及圖3-4所示：圖3-3 輸出數(shù)據(jù) 圖3-4 圓圖3.3 仿真結(jié)果驗(yàn)算Ø 負(fù)載導(dǎo)納y：y=1/r=1/(50+j*0)=0.02Ø 駐波比：=z0*z0/zl （zl=r=50）=50*50/50=50Ø 反射系數(shù)：=（-1）/（+1） =49/50=0.960784Ø 行波系數(shù)k： k=1/=1/50=0.02由上述分析可知，仿真結(jié)果符合定義，結(jié)果正確！第4章 設(shè)計小結(jié)由于我們微波技術(shù)的課程只有短短9周的學(xué)習(xí)時間，因

16、而所學(xué)到的知識也是較為膚淺的，因此，這門課程設(shè)計對我而言既是深入了解微波技術(shù)專業(yè)知識的機(jī)遇，同時也是一個挑戰(zhàn)。通過本次課程設(shè)計，在編程過程中對smith圓圖的有關(guān)知識進(jìn)行了很好的復(fù)習(xí)和掌握。由初步的只能畫出簡單的圓圖，得到對應(yīng)的電抗、導(dǎo)納的值到可以通過電抗、導(dǎo)納的值畫出相應(yīng)的smith 圓圖。而且在程序的編寫及調(diào)試過程中，對matlab的知識也進(jìn)行了很好的運(yùn)用與學(xué)習(xí)。開始時我覺得對smith圓圖知識的掌握還可以，由于沒有專門地去學(xué)，因而在matlab軟件的應(yīng)用上有所困難。最后我是在相關(guān)專業(yè)同學(xué)的指點(diǎn)及自己的摸索下才比較完整的完成該軟件的設(shè)計。該軟件能實(shí)現(xiàn)smith圓圖軟件的基本功能，不足之處也

17、存在著，比如要在分析傳輸線理論時會有所欠缺；另外界面不夠美觀、較為單一，放大鏡功能也沒有能夠加進(jìn)。我覺得以后可以多進(jìn)行類似的課程設(shè)計，給我們以足夠的實(shí)踐空間，從而借以加深我們對所學(xué)知識的理解和掌握。參考文獻(xiàn)1劉志學(xué). 導(dǎo)抗圓圖計算和廣播電視網(wǎng)絡(luò)傳輸匹配設(shè)計 中國有線電視。20042 劉學(xué)觀，郭輝萍. 微波技術(shù)與天線（第二版） 西安：西安電子科技大學(xué)出版社. 20063 喬曉華. 基于matlab語言的smith圓圖軟件設(shè)計. 中國有線電視. 20034 車晴. 電子系統(tǒng)仿真與matlab 北京：北京廣播學(xué)院出版社.2000附錄部分源程序代碼：23function varargout = smi

18、th(varargin)%smith m-file for smith.fig% smith, by itself, creates a new smith or raises the existing% singleton*.% h = smith returns the handle to a new smith or the handle to% the existing singleton*.% smith('property','value',.) creates a new smith using the% given property value

19、pairs. unrecognized properties are passed via% varargin to smith_openingfcn. this calling syntax produces a% warning when there is an existing singleton*.% smith('callback') and smith('callback',hobject,.) call the% local function named callback in smith.m with the given input% argum

20、ents.% *see gui options on guide's tools menu. choose "gui allows only one% instance to run (singleton)".% see also: guide, guidata, guihandles% edit the above text to modify the response to help smith% last modified by guide v2.5 12-jun-2009 14:40:47% begin initialization code - do no

21、t editgui_singleton = 1;gui_state = struct('gui_name', mfilename, . 'gui_singleton', gui_singleton, . 'gui_openingfcn', smith_openingfcn, . 'gui_outputfcn', smith_outputfcn, . 'gui_layoutfcn', , . 'gui_callback', );if nargin && ischar(varargin1

22、) gui_state.gui_callback = str2func(varargin1);endif nargout varargout1:nargout = gui_mainfcn(gui_state, varargin:);else gui_mainfcn(gui_state, varargin:);end% end initialization code - do not edit% - executes just before smith is made visible.function smith_openingfcn(hobject, eventdata, handles, v

23、arargin)% this function has no output args, see outputfcn.% hobject handle to figure% eventdata reserved - to be defined in a future version of matlab% handles structure with handles and user data (see guidata)% varargin unrecognized propertyname/propertyvalue pairs from the% command line (see varar

24、gin)global phaseangle max grconstant = linspace(0,10,5);phaseangle = linspace(0,2*pi,50);unitgamma = exp(j*phaseangle);%plot the unit circle in the complex planehold on;plot(real(unitgamma),imag(unitgamma),'r');%set(gcf,'position',0 0 1280 990);axis squarezoom onaxis(-1.1 1.1 -1.1 1.

25、1);max=2001;bound2=0;bound3=0;min_bound1=0;min_bound2=0;max_bound2=0;h=0;word=0;gr = linspace(-1,1,max);hold on;interval = .01:.01:.2,.22:.02:.5,.55:.05:1,1.1:.1:2,2.2:.2:5,6:10,12:2:20,30:10:50;interval2= .01:.01:.5,.55:.05:1,1.1:.1:2,2.2:.2:5,6:10,12:2:20,30:10:50;%plot real axisplot(gr, zeros(1,l

26、ength(gr),'r');%equations were derived using the symbolic toolbox as follows%solve('r=(1-gr2-gi2)/(1-gr)2+gi2)','gi')%bound was derived as follows%solve('1/(r+1)*(-(r+1)*(r-2*r*gr+r.*gr2-1+gr2)(1/2)=0','gr')for r = interval2, min_bound1 = (r-1)/(r+1); if(r<

27、.2) if(mod(r,.1)=0) max_bound = (-1+22+r2)/(22+r2+2*r+1); elseif(mod(r,.02)=0) max_bound = (-1+.52+r2)/(.52+r2+2*r+1); else max_bound = (-1+.22+r2)/(.22+r2+2*r+1); if(r=.05 | (r<.151 & r>.149) min_bound2 = (-1+.52+r2)/(.52+r2+2*r+1); max_bound2 = (-1+12+r2)/(12+r2+2*r+1); end end elseif(r&

28、lt;1) if(mod(r,.2)=0) max_bound = (-1+52+r2)/(52+r2+2*r+1); elseif(mod(r,.1)=0) max_bound = (-1+22+r2)/(22+r2+2*r+1); elseif(r=.25 | r=.35 | r=.45) temp = (-1+.52+r2)/(.52+r2+2*r+1); min_bound2 = max(min_bound1, temp); max_bound = (-1+12+r2)/(12+r2+2*r+1); elseif(r<.5) max_bound = (-1+.52+r2)/(.5

29、2+r2+2*r+1); else max_bound = (-1+12+r2)/(12+r2+2*r+1); end elseif(r<5) if(mod(r,2)=0) max_bound = (-1+202+r2)/(202+r2+2*r+1); elseif(mod(r,1)=0) max_bound = (-1+102+r2)/(102+r2+2*r+1); elseif(r>2) max_bound = (-1+52+r2)/(52+r2+2*r+1); else if(mod(r,.2)=0) max_bound = (-1+52+r2)/(52+r2+2*r+1);

30、 else max_bound = (-1+22+r2)/(22+r2+2*r+1); end end elseif(r<10) if(mod(r,2)=0) max_bound = (-1+202+r2)/(202+r2+2*r+1); else max_bound = (-1+102+r2)/(102+r2+2*r+1); end else if(r=10|r=20) max_bound = (-1+502+r2)/(502+r2+2*r+1); elseif(r=50) max_bound = 1; elseif(r<20) max_bound = (-1+202+r2)/(

31、202+r2+2*r+1); else max_bound = (-1+502+r2)/(502+r2+2*r+1); end end index = ceil(min_bound1+1)*(max-1)/2+1); actual_value = gr(index); if(actual_value<min_bound1) index = index + 1; end min=index; index = ceil(max_bound+1)*(max-1)/2+1); actual_value = gr(index); if(actual_value>max_bound) inde

32、x = index - 1; end min2 = ceil(min_bound2+1)*(max-1)/2+1); actual_value = gr(min2); if(actual_value<min_bound2) min2 = min2 + 1; end max2 = ceil(max_bound2+1)*(max-1)/2+1); actual_value = gr(max2); if(actual_value<max_bound2) max2 = max2 + 1;end r_l_a=1/(r+1)*(-(r+1)*(r-2*r.*gr(min:index)+r.*g

33、r(min:index).2-1+gr(min:index).2).(1/2); r_l_b=-1/(r+1)*(-(r+1)*(r-2*r.*gr(min:index)+r.*gr(min:index).2-1+gr(min:index).2).(1/2); r_l_b(1)=0; r_l_a(1)=0; r_l_a2=1/(r+1)*(-(r+1)*(r-2*r.*gr(min2:max2)+r.*gr(min2:max2).2-1+gr(min2:max2).2).(1/2); r_l_b2=-1/(r+1)*(-(r+1)*(r-2*r.*gr(min2:max2)+r.*gr(min

34、2:max2).2-1+gr(min2:max2).2).(1/2); r_l_a3=1/(r+1)*(-(r+1)*(r-2*r.*gr(min2:index)+r.*gr(min2:index).2-1+gr(min2:index).2).(1/2); r_l_b3=-1/(r+1)*(-(r+1)*(r-2*r.*gr(min2:index)+r.*gr(min2:index).2-1+gr(min2:index).2).(1/2); %fix resolution issues in .2-.5 range if(r>.2 & r<.5 & (mod(r,.

35、02)=0)if(r=1) color = 'r' else color ='b' end plot(gr(min:index),r_l_a(1:index-min+1),color,gr(min:index), r_l_b(1:index-min+1),color); if(r<=1) if(mod(r,1)=0) word = num2str(r) '.0' else word = num2str(r); end if(mod(r,.1)=0) set(text(gr(min),0,word),'rotation',90

36、,'horizontalalignment','left','verticalalignment','bottom'); end elseif(r<=2) if(mod(r,.2)=0) if(mod(r,1)=0) word = num2str(r) '.0' else word = num2str(r); end set(text(gr(min),0,word),'rotation',90,'horizontalalignment','left','

37、verticalalignment','bottom'); end elseif(r<=5) if(mod(r,1)=0) set(text(gr(min),0,num2str(r) '.0'),'rotation',90,'horizontalalignment','left','verticalalignment','bottom'); end else if(mod(r,10)=0) set(text(gr(min),0,num2str(r),'rotat

38、ion',90,'horizontalalignment','left','verticalalignment','bottom'); end end elseif(r=.25 | r=.35 | r=.45) plot(gr(min2:index),r_l_a3,'b'); plot(gr(min2:index),r_l_b3,'b'); end if(r=.05 | (r>.149 & r<.151) plot(gr(min2:max2),r_l_a2(length(

39、gr(min2:max2)-length(r_l_a2)+1:length(r_l_a2),'b'); plot(gr(min2:max2),r_l_b2(length(gr(min2:max2)-length(r_l_b2)+1:length(r_l_b2),'b'); end end%equations were derived using the symbolic toolbox as follows%solve('2*gi/(1-gr)2+gi2)=x','gi')%bound was derived as follows

40、%solve('1-x2+2*x2*gr-x2*gr2=0','gr')%solve('1/2/x*(2+2*(1-x2+2*x2*gr-x2*gr2)(1/2)=(1-gr2)(1/2)','gr')for x = interval, inter_bound = (-1+x2)/(x2+1); %intersection with unit circle: all values must be less than this imag_bound = (-1+x)/x; %boundary of imagination: all

41、values must be greater than this angle_point = 0; if(inter_bound = 0) angle_point = sqrt(1-inter_bound2)/inter_bound; end imag_bound_y = 1/2/x*(-2+2*(1-x2+2*x2.*inter_bound-x2.*inter_bound.2).(1/2); imag_rad = (imag_bound2 + imag_bound_y2)(1/2); condition = imag_rad < 1; if(inter_bound > 1) in

42、ter_bound = 1; elseif(inter_bound < -1) imag_bound=-1; end if(imag_bound > 1) imag_bound = 1; elseif(imag_bound < -1) imag_bound=-1; end %used solve function to find intersection of appropriate circle with corresponding hyperbolics %solve('-1/(r+1)*(-(r+1)*(r-2*r*gr+r*gr2-1+gr2)(1/2)=1/

43、2/x*(-2+2*(1-x2+2*x2*gr-x2*gr2)(1/2)','gr') %the following conditional tree creates the internal bounding between the two types of curves for variable resolution if(x<.2) if(mod(x,.1)=0) max_bound = (-1+x2+22)/(x2+22+2*2+1); elseif(mod(x,.02)=0) max_bound = (-1+x2+.52)/(x2+.52+2*.5+1)

44、; else max_bound = (-1+x2+.22)/(x2+.22+2*.2+1); end elseif(x<1) if(mod(x,.2)=0) max_bound = (-1+x2+52)/(x2+52+2*5+1); elseif(mod(x,.1)=0) max_bound = (-1+x2+22)/(x2+22+2*2+1); elseif(x<.5) max_bound = (-1+x2+.52)/(x2+.52+2*.5+1); else max_bound = (-1+x2+12)/(x2+12+2*1+1); end elseif(x<5) if

45、(mod(x,2)=0) max_bound = (-1+x2+202)/(x2+202+2*20+1); elseif(mod(x,1)=0) max_bound = (-1+x2+102)/(x2+102+2*10+1); elseif(x>2) max_bound = (-1+x2+52)/(x2+52+2*5+1); else if(mod(x,.2)=0) max_bound = (-1+x2+52)/(x2+52+2*5+1); else max_bound = (-1+x2+22)/(x2+22+2*2+1); end end elseif(x<10) if(mod(

46、x,2)=0) max_bound = (-1+x2+202)/(x2+202+2*20+1); else max_bound = (-1+x2+102)/(x2+102+2*10+1); end else if(x=10|x=20) max_bound = (-1+x2+502)/(x2+502+2*50+1); elseif(x=50) max_bound = 1; elseif(x<20) max_bound = (-1+x2+202)/(x2+202+2*20+1); else max_bound = (-1+x2+502)/(x2+502+2*50+1); end end in

47、ter_index = ceil(inter_bound+1)*(max-1)/2+1); imag_index = ceil(imag_bound+1)*(max-1)/2+1); index4 = ceil(max_bound+1)*(max-1)/2+1); index1 = max(inter_index,imag_index); %maximum index for c,d index2 = min(imag_index,inter_index); %minimum index for c,d if(condition) index3=imag_index; else index3=

48、inter_index; end actual_value1 = gr(index1); actual_value2 = gr(index2); actual_value3 = gr(index3); actual_value4 = gr(index4); if(actual_value1 > inter_bound & index1 = inter_index)|(actual_value1 > imag_bound & index1 = imag_index) index1 = index1 - 1; end if(actual_value2 < inte

49、r_bound & index2 = inter_index)|(actual_value2 < imag_bound & index2 = imag_index) index2 = index2 + 1; end if(actual_value3 < inter_bound & index3 = inter_index)|(actual_value3 < imag_bound & index3 = imag_index) index3 = index3 + 1; end if(actual_value4 > max_bound) ind

50、ex4 = index4 - 1; end min=index2; max2=index1; max3=index4; min2 = index3; % actual_value1 = gr(min); % actual_value2 = gr(max2); % min=1; % max2=max; % min2=1; x_l_a = real(1/2/x*(-2+2*(1-x2+2*x2.*gr(min2:max3)-x2.*gr(min2:max3).2).(1/2); x_l_b = real(1/2/x*(2-2*(1-x2+2*x2.*gr(min2:max3)-x2.*gr(min2:max3).2).(1/2); x_l_c= real(1/2/x*(2+2*(1-x2+2*x2.*gr(min:max2)-x2.*gr(min:max2).2).(1/2); x_l_d= real(1/2/x*(-2-2*(1-x2+2*x2.*gr(min:max2)-x2.*gr(min:max2).2).(1/2); if(min2<max3) x_l_c(1)=x_l_b(1); x_l_d(1)=x_l_a(1); end check1 = abs(round(10000*1