平板傳輸線電磁仿真

1、電磁仿真綜合實(shí)踐Numerical Simulation of Electromagnetic Field作者 張泳 Author Yong Zhang指導(dǎo)老師 王秉中Instructor Bing-Zhong Wang單位 電子科大物理電子學(xué)院99777Class 99777 CPE of UESTC設(shè)計(jì)時(shí)間 2002年5月8日6月13日【摘要】仿真是通過建立系統(tǒng)模型，并且對(duì)模型進(jìn)行實(shí)際研究，以代替對(duì)于實(shí)際系統(tǒng)的研究。Matlab 5.3是一種功能強(qiáng)大工程數(shù)學(xué)軟件。本文闡明了怎樣利用Matlab 5.3實(shí)現(xiàn)一維以及二維電磁仿真實(shí)例的建模分析、編程調(diào)試、結(jié)果分析?！続bstract】Simul

2、ation is an available means in studying real system by Model-building.This paper building a numerical simulation of electromagnetic field model based on the Matlab 5.3,the design is discussed in detail.【關(guān)鍵詞】電磁仿真，計(jì)算電磁學(xué)，電磁場(chǎng)【Key Words】 Electromagnetic Simulation,Computational Electromagnetics,Electroma

3、gnetic Field【設(shè)計(jì)目的】1 課程性質(zhì)：屬于非計(jì)算機(jī)專業(yè)計(jì)算機(jī)基礎(chǔ)教育第三層次計(jì)算機(jī)綜合應(yīng)用層次的課程2 任務(wù)： 21 加強(qiáng)學(xué)生的計(jì)算機(jī)綜合應(yīng)用能力、尤其是運(yùn)用計(jì)算機(jī)分析和解決專業(yè)問題的能力培養(yǎng)，使學(xué)生對(duì)獨(dú)立進(jìn)行科學(xué)研究有初步實(shí)踐。 22 初步掌握一種純數(shù)值電磁仿真方法時(shí)域有限差分法。 23 初步學(xué)會(huì)綜合應(yīng)用一種程序設(shè)計(jì)語(yǔ)言進(jìn)行科學(xué)與工程計(jì)算。 24 增強(qiáng)科技報(bào)告的寫作能力，學(xué)會(huì)相關(guān)軟件使用。正文第一部分：電磁仿真的意義現(xiàn)代科學(xué)研究的基本模式：“科學(xué)實(shí)驗(yàn)，理論分析，數(shù)值計(jì)算”三位一體。計(jì)算電磁學(xué)：以電磁場(chǎng)理論為基礎(chǔ)，以計(jì)算機(jī)及計(jì)算機(jī)技術(shù)為工具和手段，運(yùn)用計(jì)算科學(xué)所提供的各種方法，解決

4、復(fù)雜電磁問題的一門應(yīng)用科學(xué)。電磁場(chǎng)工程電磁場(chǎng)理論 電磁規(guī)律 實(shí)驗(yàn)驗(yàn)證 數(shù)學(xué)方程計(jì)算電磁學(xué) 計(jì)算數(shù)據(jù) 建模與仿真 運(yùn)算手段和方法 優(yōu)化與設(shè)計(jì)【電場(chǎng)問題求解方法】解析法電磁場(chǎng)問題的傳統(tǒng)解法。強(qiáng)調(diào)電磁分析和數(shù)學(xué)分析通常給出緊湊高效的計(jì)算程序。程序的最終擁護(hù)只具有很少的彈性，只能改變很少的參數(shù)，結(jié)構(gòu)的主要特性已被編進(jìn)程序，適用于專用程序開發(fā)。純數(shù)值法直接以數(shù)值的、程序的形式描述電磁場(chǎng)問題。普適性強(qiáng)，用戶彈性大。特定問題的邊界條件、電氣結(jié)構(gòu)、激勵(lì)等特性可以不編入基本程序，而由用戶輸入，更好的情況是通過圖形界面輸入。用戶不必具備高深的電磁理論、數(shù)學(xué)及數(shù)值技術(shù)方面的知識(shí)。適用于通用商業(yè)程序開發(fā)。受到硬件限制

5、大。 FDTD TLM FEM 惠更斯原理 變分原理麥克斯韋方程數(shù)學(xué)描述 及 離散化處理結(jié)果計(jì)算計(jì)算程序數(shù)值模型 邊界條件材料特性 后處理 用戶界面 第二部分 時(shí)域有限差分法基本原理【差分原理】前項(xiàng)差分 誤差 O(h)后項(xiàng)差分 誤差 O(h)中心差分 誤差 O(h*h) 用中心差分代替對(duì)空間、時(shí)間坐標(biāo)的微分 【Yee的差分網(wǎng)格】每個(gè)磁場(chǎng)分量被四個(gè)電場(chǎng)分量環(huán)繞，每個(gè)電場(chǎng)分量被四個(gè)磁場(chǎng)分量環(huán)繞。E和H的空間取樣位置相差半個(gè)空間步長(zhǎng)，E和H時(shí)間取樣位置相差半個(gè)時(shí)間步長(zhǎng)。Yee的差分算法考慮空間一個(gè)無(wú)緣區(qū)域,其媒質(zhì)的參數(shù)不隨時(shí)間變化且各項(xiàng)同性,則麥克斯韋旋度方程在直角坐標(biāo)系中,寫成分量式:= (1)

6、(2) (3) (4) (5) (6)【蛙跳格式】H0 H1/2 H3/2E0 E1 E2【解的穩(wěn)定條件】：； 【數(shù)值色散】 在時(shí)域有限差分網(wǎng)格中，數(shù)值波模的傳播速度將隨頻率改變，即有色散。這種色散由數(shù)值網(wǎng)格引起，而非物理上客觀存在。為了減小色散，實(shí)用中通常取空間步長(zhǎng)滿足 。第三部分 均勻平行板傳輸線傳輸特性仿真 FDTD基本過程【平行板傳輸線中的主模TEM模】根據(jù)平行板傳輸線中TEM波的特性對(duì)MAXWELL旋度方程簡(jiǎn)化得： 其中 若取時(shí)域波形為高斯波形， ，如下圖所示。 【一維差分格式】E【網(wǎng)格劃分】Ex Hy Ex Hy Ex Hy Ex Hy Ex k=0 k=1/2 k=1 k=3/2

7、 k=2 k=Nz-1/2 k=Nz T=0.5ns，則=1GHz，d=0.18m,L=6m, =0.3m, =0.015m, 400, 【時(shí)間步進(jìn)過程】在初始（n=0）時(shí)令所有場(chǎng)兩為零。將高斯脈沖激勵(lì)加在輸入口 由FDTD方程算出n+1/2時(shí)刻所有場(chǎng)點(diǎn)的磁場(chǎng) 由FDTD方程算出n+1時(shí)刻所有場(chǎng)點(diǎn)的電場(chǎng) 令金屬邊界上的切向電場(chǎng)分量為零，并在截?cái)噙吔缟鲜褂媒財(cái)噙吔鐥l件 記錄并輸出兩個(gè)觀測(cè)面處的電壓值 n n+1重復(fù)上一步驟，直到（n=nt）脈沖波全部傳出計(jì)算區(qū)域。 【截?cái)噙吔鐥l件】 終端短路： 終端匹配：（吸收邊界條件） 【一維仿真實(shí)例】終端匹配：入射波大約經(jīng)過1000步后到達(dá)邊界，由于采用吸收

8、邊界條件，高斯波幾乎被完全吸收，通過圖可以觀測(cè)到吸收后只有振幅小于的波形存在，這在誤差的允許范圍之內(nèi)。第一條、第二條分別為第一、二觀測(cè)面的電壓。時(shí)間步長(zhǎng)sZ方向網(wǎng)格數(shù)（0015m）Z方向網(wǎng)格數(shù)（0015m）電壓強(qiáng)度V電場(chǎng)強(qiáng)度V/M電場(chǎng)強(qiáng)度V/M第二條第一條 入射波 被吸收后的波形 觀測(cè)面電壓終端短路：入射波到達(dá)終端后由于截?cái)噙吔鐥l件采取的是短路，所以波形反轉(zhuǎn)，向入射口傳輸，由于Nt取得足夠大（Nt=2000），波形會(huì)到達(dá)入射口，所以應(yīng)該在波形到達(dá)入射口之前，激勵(lì)脈沖已完全進(jìn)入入射口后，將入射口邊界條件切換為吸收邊界條件。第一條第二條分別為第一、二觀測(cè)面的電壓。在橫軸下方的是反射電壓，上方為入射

9、電壓。第二條電場(chǎng)強(qiáng)度V/M時(shí)間步長(zhǎng)sZ方向網(wǎng)格數(shù)（0015m）Z方向網(wǎng)格數(shù)（0015m）第一條第一條電壓強(qiáng)度V電場(chǎng)強(qiáng)度V/M 被反射的波形 被吸收后的波形 觀測(cè)面電壓第四部分 帶擋板的平行板傳輸線傳輸特性仿真X L/2 d s=d/3Z L【二維FDTD方程】由于不連續(xù)性的出現(xiàn)，沿X方向場(chǎng)的均勻性遭到破壞，沿Y方向場(chǎng)的均勻性仍然保持，此時(shí)平行板傳輸線中的場(chǎng)分量有Ex,Ez,Hy.【參數(shù)設(shè)置】：數(shù)組EX(NX,NZ+1),EZ(NX+1,NZ),HY(NX,NZ),L/2應(yīng)選足夠長(zhǎng)，使得在擋板反射波返回入射口之前，激勵(lì)脈沖已完全進(jìn)入入射口，并且入射口邊界條件已切換為吸收邊界條件?！居^測(cè)面電壓波形

10、】【入射、反射、透射電壓頻譜】【網(wǎng)絡(luò)的S參數(shù)】現(xiàn)在我們?cè)儆涗浵驴傠妷骸⒕涂梢苑蛛x出入射電壓、反射、透射電壓。=- = 終端短路的反射系數(shù)：終端匹配的S參數(shù)： 【終端條件】終端匹配終端短路=0【觀測(cè)面電壓波形】【入射、反射、透射電壓頻譜】【二維仿真實(shí)例】 終端匹配：可以通過對(duì)Ex波形的觀測(cè)知道，波形到達(dá)擋板時(shí)產(chǎn)生尖端效應(yīng)，高斯波有一部分被“擠”過擋板，有一部分被擋板反射。此時(shí)應(yīng)該將入射口和終端都設(shè)為吸收邊界條件。 如下圖所示：Z方向網(wǎng)格數(shù)(0.015m)X方向網(wǎng)格數(shù)(0.015m)電場(chǎng)強(qiáng)度Ex V/MZ方向網(wǎng)格數(shù)(0.015m)X方向網(wǎng)格數(shù)(0.015m)電場(chǎng)強(qiáng)度Ex V/M X方向網(wǎng)格數(shù)(0.

11、015m) 入射波 到達(dá)擋板 電場(chǎng)強(qiáng)度Ex V/MZ方向網(wǎng)格數(shù)(0.015m)Z方向網(wǎng)格數(shù)(0.015m)X方向網(wǎng)格數(shù)(0.015m)X方向網(wǎng)格數(shù)(0.015m)電場(chǎng)強(qiáng)度Ex V/M 通過擋板 兩端都吸收電壓幅度(V)時(shí)間步長(zhǎng)（s）第一個(gè)觀測(cè)面反射電壓第二個(gè)觀測(cè)面透射電壓第一個(gè)觀測(cè)面入射電壓 參考面電壓 S參數(shù)頻率f (0.1GHz)終端匹配時(shí)S參數(shù)終端短路：由于擋板與短路板之間構(gòu)成一種諧振腔結(jié)構(gòu)，時(shí)域脈沖波在其間來(lái)回反射，入射口仍然是吸收邊界條件，波完全消失的時(shí)間會(huì)比終端匹配情況要長(zhǎng)。Nt足夠大（Nt=20000）如下圖所示：電場(chǎng)強(qiáng)度ExV/MX方向網(wǎng)格數(shù)(0.015m)z方向網(wǎng)格數(shù)(0.0

12、15m) 終端短路、入射口吸收，Ex被終端反射的波形電壓幅度V 頻率f (0.1GHz) 終端短路時(shí)第一參考面電壓 反射系數(shù)頻率f (0.1GHz) 終端短路時(shí)反射系數(shù) 第五部分 課程總結(jié)與討論短短幾周的電磁仿真綜合實(shí)踐很快就結(jié)束了，這幾周緊張的學(xué)習(xí)給我留下了深刻的印象。我們已經(jīng)學(xué)習(xí)了電磁場(chǎng)與電磁波、微波技術(shù)基礎(chǔ)等課程，對(duì)電磁理論有一定的認(rèn)識(shí)和了解，但那只是完全理論性質(zhì)的。通過這門課程的學(xué)習(xí)，我對(duì)電磁場(chǎng)等相關(guān)知識(shí)有了進(jìn)一步的了解，有了一定感性的認(rèn)識(shí)，感覺電磁場(chǎng)也不是那么和天書一樣的難以理解。通過用計(jì)算機(jī)模擬將電磁現(xiàn)象展示在我們面前，使得我不但加深了對(duì)電磁理論的理解，同時(shí)也鍛煉了我們計(jì)算機(jī)應(yīng)用方

13、面的能力，對(duì)科學(xué)研究有了一定的了解。我有如下幾點(diǎn)體會(huì)：首先，這次實(shí)踐使得我對(duì)計(jì)算機(jī)的應(yīng)用有了進(jìn)一步的提高。以往我對(duì)計(jì)算機(jī)的認(rèn)識(shí)只是處于一個(gè)很簡(jiǎn)單的層次上，對(duì)其應(yīng)用、原理也只是略知皮毛。但通過這段時(shí)間的學(xué)習(xí)，自己用計(jì)算機(jī)以及相關(guān)軟件完成了這次實(shí)踐活動(dòng)。雖然在開始的時(shí)候?qū)浖约霸砝斫獠粔颍峭ㄟ^老師的多次講解，同學(xué)之間的互相探討，漸漸地自己掌握了一定的方法?？粗约壕幊龅某绦蛘{(diào)試出令人滿意的結(jié)果，還是相當(dāng)喜悅的。其次，這次的實(shí)踐活動(dòng)讓我加深了對(duì)專業(yè)知識(shí)的了解，對(duì)電磁理論也不象以前那樣見到就躲了。當(dāng)我用自己編的程序捕捉到了它們，對(duì)它們的產(chǎn)生、傳播、反射、透射等現(xiàn)象有了更為直接的認(rèn)識(shí)和了解。自己

14、本來(lái)對(duì)電磁理論有一定的興趣，苦于沒有很好的方法幫助自己更好的理解和學(xué)習(xí)，或許通過這次實(shí)踐活動(dòng)，自己以后能夠找到適合自己的學(xué)習(xí)電磁理論的方法。第三，通過對(duì)原理、程序的不斷學(xué)習(xí)、操作，使得我認(rèn)識(shí)到做任何事情都要認(rèn)真，持之以恒。王老師科學(xué)嚴(yán)謹(jǐn)?shù)淖黠L(fēng)給我留下了深刻的印象，不論是做人還是做科研，王老師都是那么的一絲不茍。這一點(diǎn)讓我也不敢有所懈怠。同時(shí)我也感覺到自己能力的欠缺，分析問題、解決問題的能力還有待提高。相信在以后我會(huì)通過自己的不斷努力，提高自己的能力，學(xué)習(xí)和工作上更加細(xì)致認(rèn)真。第四，這次實(shí)踐活動(dòng)，我感覺是我上大學(xué)以來(lái)收獲最大的一門課程。王老師講課不是傳統(tǒng)的灌輸式的教育，學(xué)生是真正意義上學(xué)習(xí)的主體

15、。課前需要我們對(duì)理論方面的知識(shí)多加復(fù)習(xí)和學(xué)習(xí)，上課的時(shí)候需要我們積極的開動(dòng)腦筋，如何把理論轉(zhuǎn)化為實(shí)際可見，課后也需要我們不斷的思考，如何才能使得我們做的軟件更加趨于完善。在掌握計(jì)算電磁學(xué)有關(guān)理論知識(shí)的同時(shí)，更加重要的是我們學(xué)到了方法，這一點(diǎn)我覺得才是最為重要的。下面說(shuō)一下我對(duì)這門課程的一些建議吧。我覺得要消除一些學(xué)生對(duì)這門課程的緊張情緒（開始我就有這樣的感覺），適當(dāng)?shù)脑黾诱n時(shí)數(shù)，畢竟我們目前涉及的真的只是皮毛而已。增加課時(shí)數(shù)，可以讓大家對(duì)理論的理解加深，同時(shí)也會(huì)有更多的同學(xué)以后投入到電磁場(chǎng)領(lǐng)域方面的學(xué)習(xí)和研究。還有上機(jī)時(shí)間最好不要靠得太近了，畢竟學(xué)習(xí)有一個(gè)消化吸收的過程。這次的科研實(shí)踐活動(dòng)是令

16、人愉快的。它讓我體會(huì)到了快樂，當(dāng)然也有抓耳撓腮的痛苦；它讓我學(xué)習(xí)到了更多新的知識(shí)，新的方法；它讓我看到了自己的不足之處，在以后的學(xué)習(xí)過程中，我一定會(huì)加以改正的。最后我還要特別感謝王老師，以及研究生師姐在這次實(shí)踐活動(dòng)中給予我的幫助。正是他們一絲不茍的工作態(tài)度和嚴(yán)謹(jǐn)?shù)淖黠L(fēng)深深的教育了我，以后我也會(huì)盡我最大的努力去做好每一件事情。第六部分 附錄一、參考文獻(xiàn)1、計(jì)算電磁學(xué)（節(jié)選），王秉中，電子科技大學(xué)，2000年版2、電磁仿真綜合實(shí)踐，王秉中，電子科技大學(xué)3、MATLAB基礎(chǔ)與應(yīng)用簡(jiǎn)明教程，張平 ，北京航空航天大學(xué)出版社，2001年1月第一版二、相關(guān)程序1、一維終端匹配%定義變量t=0.5E-9;c=

17、3E8;fmax=1/(2*t);lammin=c/fmax;epslon=8.854E-12;mu=4E-7*p;y0=120*pi;dz=lammin/20;dt=dz/(2*c);t0=3*t;l=6;d=0.18;nz=400;nt=1000;k1=101;k2=301;cm=dt/(mu*dz);ce=dt/(epslon*dz);cmur=(dt*c-dz)/(dt*c+dz);ex(1:nz+1)=0;hy(1:nz)=0;vref1(1:nt)=0;vref2(1:nt)=0;%程序for n=1:nt ex(1)=exp(-(n*dt-t0)2/t2); for i=1:nz

18、 hy(i)=hy(i)+cm*(ex(i)-ex(i+1); end exm=ex(nz); for j=2:nz ex(j)=ex(j)+ce*(hy(j-1)-hy(j); end ex(nz+1)=exm+cmur*(ex(nz)-ex(nz+1);%一維匹配終端吸收條件 vref1(n)=ex(k1)*d; vref2(n)=ex(k2)*d; %plot(1:nz+1,ex);pause(0.001);%電場(chǎng)強(qiáng)度監(jiān)控 x=1:nt; plot(x,vref1(x),'r',x,vref2(x),'g');pause(0.001);%觀測(cè)面電壓圖 %g

19、ridend%一維終端匹配2、一維終端短路%定義變量y0=120*pi;dz=lammin/20;dt=dz/(2*c);t0=3*t;l=6;d=0.18;nz=400;nt=2000;k1=101;k2=301;cm=dt/(mu*dz);ce=dt/(epslon*dz);cmur=(dt*c-dz)/(dt*c+dz);ex(1:nz+1)=0;hy(1:nz)=0;vref1(1:nt)=0;vref2(1:nt)=0;for n=1:nt for i=1:nz hy(i)=hy(i)+cm*(ex(i)-ex(i+1); end if n<=800 ex(1)=exp(-(n

20、*dt-t0)2/t2); for j=2:nz ex(j)=ex(j)+ce*(hy(j-1)-hy(j); endelse exm=ex(2); for j=2:nz ex(j)=ex(j)+ce*(hy(j-1)-hy(j); end ex(nz+1)=0; %終端短路邊界條件 ex(1)=exm+cmur*(ex(2)-ex(1); %高斯波完全進(jìn)入后，入口改為吸收邊界條件end vref1(n)=ex(k1)*d; vref2(n)=ex(k2)*d; x=1:nt;plot(x,vref1(x),'r',x,vref2(x),'g');%兩個(gè)觀測(cè)面的

21、電壓波形 pause(0.01);% plot(1:nz+1,ex);pause(0.001);%實(shí)時(shí)監(jiān)控Ex波形 %gridend%一維終端短路3、二維終端匹配%定義變量t=0.5E-9;c=3E8;fmax=1/(2*t);lammin=c/fmax;epslon=8.854E-12;mu=4E-7*pi;y0=120*pi;dz=lammin/20;dx=lammin/20;dt=dz/(2*c);t0=3*t;l=6;d=0.18;nz=400;nt=1000;k1=101;k2=301;cm=dt/(mu*dz);ce=dt/(epslon*dz);cf=dt/(epslon*dx)

22、;cmur=(dt*c-dz)/(dt*c+dz);nx=12;ez(1:nx+1,1:nz)=0;ex(1:nx,1:nz+1)=0;hy(1:nx,1:nz)=0;ver1in(1:nt)=0;ver1re(1:nt)=0;ver2tr(1:nt)=0;exb(1:nx,1)=0;exe(1:nx,1)=0;nf=10;df=1e8;jj=sqrt(-1);gver1in(1:nf)=0;gver1re(1:nf)=0;gver2tr(1:nf)=0;s11(1:nf)=0;s21(1:nf)=0;for n=1:nt %時(shí)間循環(huán) for i=1:nx %二維FDTD過程 for k=1:

23、nz hy(i,k)=hy(i,k)+cm*(ez(i+1,k)-ez(i,k)*dz/dx+ex(i,k)-ex(i,k+1); end end if n<=300 ex(1:12,1)=exp(-(n*dt-t0)2/t2); for i=1:nx for k=2:nz ex(i,k)=ex(i,k)+ce*(hy(i,k-1)-hy(i,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; %設(shè)置擋板 for i=2:nx for k=1:nz ez(i,k)=ez(i,k)+cf*(hy

24、(i,k)-hy(i-1,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; else %高斯波完全進(jìn)入，入射端改為吸收邊界條件 exb(1:nx,1)=ex(1:nx,2);exe(1:nx,1)=ex(1:nx,nz); for i=1:nx for k=2:nz ex(i,k)=ex(i,k)+ce*(hy(i,k-1)-hy(i,k); end end for i=2:nx for k=1:nz ez(i,k)=ez(i,k)+cf*(hy(i,k)-hy(i-1,k); end endex

25、(1:nx,1)=exb(1:nx,1)+cmur*(ex(1:nx,2)-ex(1:nx,1);%入口邊界條件ex(1:nx,nz+1)=exe(1:nx,1)+cmur*(ex(1:nx,nz)-ex(1:nx,nz+1);%吸收邊界條件 ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0;endfor i=1:nx %分離出入射電壓、反射電壓、透射電壓 ver2tr(n)=ver2tr(n)+ex(i,301)*d/12; end if n<=500 for i=1:nx ver1in(n)=ver1in(n)+

26、ex(i,101)*d/12; end else for i=1:nx ver1re(n)=ver1re(n)+ex(i,101)*d/12; end end %x=1:nt%plot(x,ver1in(x),'r',x,ver1re(x),'b',x,ver2tr(x),'m');%各個(gè)電壓波形surf(ex); %監(jiān)控Ex波形pause(0.01)endfor f=1:nf for n=1:nt %傅氏變換 gver1in(f)=gver1in(f)+dt*ver1in(n)*exp(-2*pi*f*df*n*dt*jj); gver1re(

27、f)=gver1re(f)+dt*ver1re(n)*exp(-2*pi*f*df*n*dt*jj); gver2tr(f)=gver2tr(f)+dt*ver2tr(n)*exp(-2*pi*f*df*n*dt*jj); endendfor f=1:nf %S參數(shù) s11(f)=gver1re(f)/gver1in(f); s21(f)=gver2tr(f)/gver1in(f);endfor i=1:nf s1(i)=abs(s11(i); s2(i)=abs(s21(i);end%x=1:nf;%plot(x,s1,'g',x,s2,'r');%S參數(shù)%g

28、rid %二維終端匹配4、二維終端短路%定義變量t=0.5E-9;c=3E8;fmax=1/(2*t);lammin=c/fmax;epslon=8.854E-12;mu=4E-7*pi; y0=120*pi;dz=lammin/20;dx=lammin/20;dt=dz/(2*c); t0=3*t;l=6;d=0.18;nz=400; nt=20000;k1=101;k2=301; cm=dt/(mu*dz);ce=dt/(epslon*dz); cf=dt/(epslon*dx);cmur=(dt*c-dz)/(dt*c+dz);nx=12; ez(1:nx+1,1:nz)=0;ex(1:

29、nx,1:nz+1)=0;hy(1:nx,1:nz)=0; vref1(1:nt)=0;vref2(1:nt)=0;ver1in(1:nt)=0; ver1re(1:nt)=0;ver2tr(1:nt)=0;exb(1:nx,1)=0;exe(1:nx,1)=0; nf=10;df=1e8;jj=sqrt(-1); gver1in(1:nf)=0;gver1re(1:nf)=0;gver2tr(1:nf)=0; s11(1:nf)=0;s21(1:nf)=0;for n=1:nt for i=1:nx %二維FDTD過程 for k=1:nz hy(i,k)=hy(i,k)+cm*(ez(i+

30、1,k)-ez(i,k)*dz/dx+ex(i,k)-ex(i,k+1); end end if n<=300 ex(1:12,1)=exp(-(n*dt-t0)2/t2); for i=1:nx for k=2:nz ex(i,k)=ex(i,k)+ce*(hy(i,k-1)-hy(i,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; %設(shè)置擋板 for i=2:nx for k=1:nz ez(i,k)=ez(i,k)+cf*(hy(i,k)-hy(i-1,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; else %高斯波完全進(jìn)入，入射端改為吸收邊界條件 ex