1523090027唐江梅电磁场数值分析论文
- 格式:doc
- 大小:114.00 KB
- 文档页数:8
东南大学
硕士学位论文
插值与拟合技术在电磁场频域和谱域问题中的应用
姓名:王彩芹
申请学位级别:硕士
专业:电磁场与微波技术
指导教师:周后型
20070101
东南大学硕士学位论文
分为三个步骤。
第一步,确定指数函数的项数肘值.
假设在Sommeffeld积分的被积函数厂(,)半周期的采样点个数m,贝塞尔函数的近似半周期为q=n'Ip。
.当对(4.10)式左端,(f)进行均匀采样时,取采样步长为△丁=qlm,则有
Ⅳ
乃=,(p△丁)=∑R矽,p=0,1,…,N一1(4.Is)
l-l
式中极点毛=P枷(f=1,2,---,M),复指数s一般具有负的实部;Ⅳ为采样点总个数,其确定公式N=mK,K为有限整数。
依采样值‘(p=o,1,…,N-))定义两个矩阵fYl】和【X】如下:
【Y2】=
【Yl】=石Z
Z五
::
(4.19)
“.20)
三称为罚参数,取值范围为耐≤三≤Ⅳ一M。
虽然£<Ⅳ,但是,Ⅳ个采样值在构造矩阵【YI】和【Y2】时都会被用到。
将这两个矩阵按以下方式进行分解
fY2】=瞄】陋儿Zo】f五】(4.21)
fY】=[ZI]IR]IZ,】一’(4.22)式中
【ZI】=1
毛
zy也4)
1
乃
;
z,。
1’
1]
20I
;f
,‘Ⅳ.¨)I
‘lJ‘Ⅳ一Lh材
L
£
卜
砖
止
也
肌
1L—li■。
第33卷第4期电子科技大学学报V ol.33 No.4 2004年8月Journal of UEST of China Aug. 2004 电磁波在等离子体层中衰减的数值分析王华东,鄢 扬,郭 军(电子科技大学物理电子学院成都 610054)【摘要】采用分层介质方法处理非均匀等离子体层,研究了电磁波射向覆盖磁化等离子体层的金属平板时电磁波的衰减特性。
着重讨论Epstein密度分布的等离子体层,分析了等离子体电子密度、电子碰撞频率、外加磁场和入射角度等因素对电磁波衰减的影响,得出在一定等离子体密度和电子碰撞频率下,等离子体对电磁波的吸收特性。
关键词电磁波; 磁化等离子体; 反射; 衰减中图分类号TN011; O539 文献标识码 ANumerical Analysis of Attenuation of Electromagnetic Wavein Non-Uniform Magnetized Plasma SlabWang Huadong,Yan Yang,Guo Jun(School of Physical Electronics, UEST of China Chengdu 610054)Abstract The attenuation of the electromagnetic wave in the non-uniform plasma slab cling to a metal plane is studying in this paper. The plasma slab is divided into a number of sub-slabs each of which is dealt with as a uniform slab. Detailed numerical calculations are conducted on Epstein-distributions.The effects of electron number density, momentum transfer collision rate, magnetic field and angle of incidence on the attenuation are taken into account, and get the character of plasma’s absorbing electromagnetic wave.Key words electromagnetic wave; magnetized plasma; reflection; attenuation电磁波与等离子体之间相互作用是等离子体物理学中重要的研究方向之一,近年来研究表明等离子体在飞行器隐身方面具有良好前景。
电磁轨道炮导轨电磁场数值分析本页仅作为文档封面,使用时可以删除This document is for reference only-rar21year.March电磁轨道炮导轨电磁场数值分析刘旭洋,李凡,刘平成,杨安波(中科院高能物理研究所,北京,100049)Numerical Analysis of Electromagnetic RailgunLiu-Xuyang,Li-Fan,Liu-Pingcheng,Yang-Anbo(Institute of High Energy Physics Chinese Academy of Sciences)摘要:电磁轨道炮作为一种新概念动能武器,由于具有发射速度高,加速时间短,能量释放易于控制等优点,使其在军事领域中有着巨大的应用潜力。
电磁轨道炮采用了电磁发射技术,本文在矩形电枢与简单轨道模型的基础上分析其磁扩散方程,以及在线电流模型下电枢的受力情况。
利用有限元分析软件Ansys-Maxwell 3D进行了轨道电枢上的电磁场计算。
关键词:电磁轨道炮;电磁场有限元计算ABSTRACT: Electromagnetic railgun as a new concept of kinetic energy weapons, has since the launch of high speed, acceleration time is short, easy to control the release of energy, etc., so that it has great potential in the military field. Electromagnetic rail gun uses electromagnetic emitting technology, on the basis of a simple rectangular armature orbit model using finite element analysis software Ansys-Maxwell 3D was the calculation of electromagnetic field on the track of the armature.KEY WORDS:electromagnetic railgun; finite element analysis of electromagnetic field作者简介刘旭洋,(1988-),男,辽宁籍,超导磁体与精密机械李凡(1993-),男,山西籍,光学方向刘平成(1993-),男,新疆籍,研究方向为核技术应用杨安波(1990-),男,湖北籍,粒子物理方向1 引言电磁发射装置是一种利用电磁力将弹丸加速到超高速的发射装置,可用于摧毁空间低轨卫星,拦截来袭导弹以及发射小型卫星等。
基于MATLAB的电磁场数值分析***(江南大学物联网工程学院,江苏无锡214036)摘要:使用计算机进行电磁场数值分析已成为电磁场的工程开发、科研和教学的重要手段。
编程实现从电磁场微分方程到有限元求解全过程需要很好的理论基础和编程技巧,难度较高。
该文介绍了电磁场数值分析的基本理论并通过两个实例介绍了使用MATLAB工具箱实现电磁场偏微分方程的有限元解法。
实验结果表明这一方法具有操作简单明了,运算速度快,计算误差可控制等优点。
关键词:电磁场数值分析有限元Electromagnetic Field Numerical Analysis Based on MATLAB***(College of Institute of Things, Jiangnan University, Wuxi, Jiangsu 214036,China) Abstract: Numerical analysis by computer is now an important way to development, research and teaching of electromagnetic field. Good theory foundation and programming skill are necessary for programming to solve the electromagnetic partial differential equation by FEM method. This article introduces the basic theory of Electromagnetic Field Numerical analysis and the fem method to solve electromagnetic partial differential equation by MATLAB PDE Toolbox. The results show that this method is easy to handle and easy to control the error with faster operation speed.Keywords: Electromagnetic Field Numerical Analysis, MATLAB, finite element method1 前言MATLAB语言是美国MathWorks公司开发的计算机软件,经过多年的发展与完善,现已成为国际上最流行的科学与工程计算的软件工具。
1.5工程电磁场数值分析的概况1.5 工程电磁场数值分析概况本章讲述了电磁场计算基础问题,内容涉及到电磁场的基本理论,各种类型电磁场的特性,所求场域中媒质的影响、场源的作用和域外源作用反映的边界条件,电磁场的规范和场解答的唯一性问题等。
作为实际电磁场问题,我们有可能通过定性分析,基于电磁场的基本方程逐步建立起相应的场的控制方程,逐步分析和建立场域边界条件和初值条件,形成定解问题,这是计算场的一种数学模型,这是求解工程电磁场的基础,我们需要努力掌握。
求解电磁场的数学模型,要想用解析的方法来确定这种定解问题的解答,那是十分困难甚至是不可能的。
因此自1864年Maxwell方程组诞生后相当长一段时间,电磁场计算问题解决极为有限。
随着电子计算机和计算技术的发展,数值计算方法逐步产生和发展用电磁场数值分析的方法解决各类电磁场问题计算才成为可能。
1.5.1电磁场数值分析的任务电磁场数值分析的基本任务是根据Maxwell方程组和电磁场的基本理论,进行:(1)建立逼近实际工程电磁场问题的连续型的数学模型;(2)采用相应的数值计算方法,经过离散化处理而转化为等价的离散型数学模型;——由离散数值构成的联立代数方程组,应用有效的代数方程组解法,计算出待求离散数学模型的离散解(即数值解);(3)依据计算结果求得所需点处的场强,或某一区域的能量、损耗分布,或某个电磁参数的值,为工程判断和优化设计提供依据。
进行电磁场数值分析,必须具备一定的数学、物理基础和有关电磁场的专门知识,特别是建立数学模型,在很大程度上还有赖于实际的工程知识和运用数值计算方法的经验的积累,使得有可能提出恰当地理想化假设,较为准确地给定问题的定解条件,建立起把握了问题实质的数学模型。
1.5.2 电磁场数值分析主要运用方法电磁场数值分析的主要内容是讲述各种实用的数值计算方法,主要有:(1)应用于微分方程数学模型的有:有限差分法(时域有限差分法)、有限元法等;(2)应用于积分方程数学模型的有:边界元法、模拟电荷法,矩量法等。
场冷和零场冷过程中磁化强度的模型化与数值模拟‘诊寝诊姘《奄低温物理V o1.21,No.5的模型化与数值模拟李世亮闻海赵忠贤中辩晓两秀丽丽釉襄蟊再蚕蚕虿聚态中心,北京1000801999年4月22日收到2斗f通过采用蒙特卡罗模拟方法和磁通运动的热激活模型(TAFM),车文计算了在场冷(Fc)和零场冷(zFc)过程中的磁场位形及磁化强度.基于简单的Ande肌模型的U(j)依赖关系.主要的实验结果都可以较好的被模拟.计算表明,在零场冷中,磁场位形非常接近Bean临界态模型.而在场冷中则相反,表现的是磁通冻结的图象零场冷和场冷中的磁化强度都很强的依赖于钉扎势和临界电流密度,而只微弱的依赖于温度变化的速度.1引言众所周知.磁通排斥是超导体的一个主要特性,因此无论是在基础研究还是高能应用方面.测量并研究超导体的磁化强度都是必要的.在各种方法中,广泛采用的是被称为零场冷(ZFC)和场冷(Fc)的方法_1].零场冷和场冷中的磁矩通常被用来决定超导转变温度J同时,也被用于研究磁通动力学.例如.对不可逆点(7,B.在其以下磁通线被钉扎的)的研究.在高能应用领域.如超导磁悬浮方面.悬浮力被认为直接和零场冷的磁矩相连J关于所谓超导永磁体的想法也是基于在场冷过程中的磁通冻结现象_5J.综上所述,对场冷和零场冷过程的详细研究和更好的理解就非常必要了.本文中,利用蒙特卡罗模拟方法并基于热激活磁通运动模型(TAFM).我们计算了磁场位形和磁化强度.我们将模拟主要的实验现象.如磁滞回线.磁化强度对温度的依赖曲线.及其与钉扎势和临界电流密度的关系.2模拟计算考虑一个无限大平面样品,沿z轴厚2n.边界为=一n及z=a.由于在穿透深度(以内的Meissner屏蔽电流的影响.B很快的沿z从外场0H降到o(H)一1).考虑到《4.因此令B(z=±4)=o[H()一1(丁)】.其中.为下临界场,1(T)=l(0)[1一(丁)].:一另外.假设在样品内部存在着N个沿z轴周期分布的钉扎中心.邻近的两个钉扎中心距离为2.z0对一个给定的钉扎中心.磁通线可以跳跃到其邻近位置.磁通线的运动采用蒙特卡罗方法来模拟.采用一个一维数组代表着有N个钉扎中心的上述样品.每个钉扎中心处的磁感应强度为B(z).在每步计算中,从随机选定一个位置选低温物理21卷出固定大小的一束磁通线,将它跳出其位置的概率与一个在0到1之间的随机数比较.如果该随机敦小于其概率,则认为磁通线跳出自己的位置;反之,磁通线将留在原来的位置.这里采用的方法和Schmck【引与V and日Beek[等采用的方法是相同的.一根磁通线发生跳跃的机率为=expPexpf一KBT1(1)=l_l(1)其中u(j)是激活能,采用Kim-Anderson模型,有u(T,,)=uo(T)-[1±(J/J(T)],其中‘+’表示反向跳跃,’一’表示正向跳跃.电流密度J由Maxwell方程得到.,(z)=一1/o-dB/d..Tc在所有的计算中,我们都将采用如下关系:U口(T)=uo(0)[1一(T/)].L(T)=L(0)[1一丁)】[.磁矩M采用下式计算1;MJ[,uoHo—B(z)】出(2)图1计算机得到的磁滞回线.其中参数选取如下:(0)/=165.(0)=10”A/m2口=0lnn1.=25×10一图1给出了采用这种方法得到的磁滞回线,与实验符合的很好.该图的参数选取为:”(o)/kBT=16.5,(0)=10”A/m2.如果没有…特别说明的情况下.以下计算中的参数将取.为:n=0.1mm.口=2.5×10m(即有40个.钉扎中心),微观振动频率.=加,ol(O)=0.1Tesla,外场oHo=0.05Tesla.’’4结果与讨论图2给出了计算得到的零场冷和场冷曲线,其中U口(0)=lOSK,(0)=10”A/m2,020*********K图2磁化率在场玲和零场过程中对温度的依赖关系Uo(O)/KBT=105K,Jf(0)=10”A/空心表示零场拎而实心表示场冷.=90K.可以看到,二者在低温区都是平坦的.下面将要说明该平坦是有不同的起源的.这里需要说明的是.uo(o)比实验值大很多(一般(0)=1o0~1000K[9).这是因为对高温超导体.实际的相邻钉扎中心距离为10纳米量级.比较一下可知,这相当于在我们假定的钉扎5期李世亮等:场冷和零场冷过程中磁化强度的模型化与数值模拟365中心问存在上千个真正的钉扎中心下面我们称真正的钉扎中心为实钉扎,而称计算中的钉扎中心为虚钉扎.如果在两个虚钉扎间存在n个实钉扎,则一根磁通线必须通过一个实钉扎才能到达下一个虚钉扎.若假设只有正向跳跃,则在两个虚钉扎中心间跳跃的机率是:立:血e-()’=l|=InU【,『I=一市lJ(3)与(1)式比较,虚钉扎势要比实钉扎势大n 倍.若考虑反向跳跃,(3)式将变得很复杂, 但结论却不会改变,即计算中的钉扎势要比实验得到的钉扎势大很多.图3(a)~(c)给出了计算得到的零场冷,场冷升温和场冷降温过程中的磁通密度位形.图形的上面部分代表可逆区,其内部的磁场是均匀的,并且始终和表面磁场保持一致.因此,可逆区的磁通动力学与过程是譬妄蚕:0.048;;;;三0o47(c)场冷升温010-0050O00050】0圈3三种过程中的磁场位形:(a)零场冷(b)场冷降温(c)场冷升温.无关的.零场冷和场冷的区别主要在低温区,零场冷的磁场位形即使在部分穿透状态下也非常接近Bean临界态,而场冷则表现出曲线形状.因此可以知道,M(T)曲线中场冷与零场冷在低温区的平坦形状有非常不同起因:对零场冷,是因为Meisslaer屏蔽效应;对场冷,是因为磁通冻结.可以注意到场冷升温与降温的磁场位形很相似:磁通线在低温时被冻结.高温时进入可逆区.实际上,在将要进入可逆区时,二者表现是不同的:在样品表面附近,场冷降温的t3()很接近直线,而场冷升温为v字形.下面将讨论一下不同参数对零场冷和场冷曲线的影响.首先,考虑参数Uo(0),如图4.图4(a)为零场冷,在低温下不同曲线圈4(a)零场冷(b)场冷过程中不同Uo(O)的磁化重合在一起i图4(b)为场冷过程,在低温下率.场冷过程中的空心表示场冷升温,实心表曲线高度与u0(0)成单调变化关系.这是因示场冷降温.接近不可逆温度处可以看到由于为对零场冷,起因为完全的Meissner屏蔽,磁场位形不同造成的场冷升温与降温的区别.与钉扎势无关;而对场冷.起因为磁通冻结,366低温物理21卷很强的依赖于钉扎强度.在零场冷中,可以看到M(T)曲线随(0)增大向高温区移动,这是因为大(0)导致运动损耗的增加.同样原因,场冷中磁化率的绝对值随U口(0)增大而减小.改变J,(0)可以得到类似的结果,如图5.还有一个会影响M(T)曲线的是变温速度如图6(a)和图6(b),若速度降低,磁弛豫增加.从我们数据可以发现,尽管温度变化速度dT变化了两个量级.它对M(T)曲线影响却不大,这说明当达到一个适合的磁场状态后,温度变化速度的影响很小.需要说明的是图6中的温度变化值是相对数值,必须和微观振动频率一起加以考虑. 81848790K848790K薹2/K10-0.60.1’01降温速度sr4{1(场{L—..篓匮i0A027j变温速度,KfI●r1:蝴持,1’l.01■——■——暇.03J——,————.————————.————,————.————一81848790K图5(a)零场冷(b)场冷过程中不同』(0)的磁化图6(a)零场冷(b)场冷过程中不同温度变化率.场冷过程中的空心表示场冷升温而实心速度的磁化率.场冷过程中的空心表示表示场玲降温.场冷升温而实心表示场冷降温.总之,本文基于热激活磁通运动模型.通过蒙特卡罗方法模拟了零场冷和场冷过程的磁化行为.主要的实验现象都可以较好的模拟.通过计算表明,在零场冷过程中,磁场位形很接近Bean临界态模型.而在场冷中则呈现出磁通冻结的花样.场冷和零场冷中的磁化行为很强的依赖于钉扎势和临界电流密度,而只微弱的依赖于温度变化速度.【1JForarevie~vsee,YYeshurtm,A.PMskeem~,A.Smulov,RevMad.,醴(1996),911.【2]U.Welt,.WK..G,W.Cazbuee,K.GvanderV ocft.J.Z,‰Reo..62(1989).12 [3]IOx:ng,Z},~LImmi,J瑚.Y.№,KKishio.Te嘟瞄,Y.Bando,MTakano.Sdo~.276{1997).硎.[4]Ylwam,H.GLee,o蛳,37(1997),807[5】HHP.Watson,I.Ycutms.s曲棚晰&-/.,10(1997),944【6]HG..~..Imack,Rie瞄en,J.G.撼 C.JvanderBe,ek,PHKes.曲C,l”(1992),337.【7]cJvanderBeek.G.J.M眦啪rn私,PHKes.HG.sdm,R.Griesm*,.讳衄C.1卯(1992),320.f8]J.R.H印∞乃_Y.Rs岫.etal,跏.Reo,B46(1993).14440[9]H.H.wm.H.G.~..Imack,R(髑锄,BDan,JRe,zror,0妇C,241(199S),353.5期李世亮等:场冷和零场冷过程中磁化强度的模型化与数值模拟367M0DELINGANDSIMULA T10N0NTHEMAGNETIZA T10NINFIELD-COOLINGANDZER0.FⅢLD.COOLINGPROCESSESS.L.”H.H.WeI1Z.X.Zhao/’,htmna/工正哪血&棚嘞.1ma~weRO衄andc.o~/o-Muzo-Phys/cs,aAScL-n0es,PO.603,&衄100080(Reoelved22Apdt.1999)ByusingM∞teCarlosh!ulafionandthelmallyactivatedfluxmotion(fM)taxi. d,wecalculatethemag~ticfieldprofileandthusthemagnetizationintheridd-cooling(Fc)and0er6eld—co~ling(ZFC)pzocesses.Based∞thesbnpleKim-Anders ~U()dependence,themajoqserimentalobservaticaq.sc∞benicelympmduced.I tisfoundthatthecalculated6eldpm~ileB(z)ksclose∞theBeancriticalsterem:xt elinthe℃pmcess,in~-acp∞mfasttOthatiI1theFCpmcesschshsw sacurvedfrozen fluxpatIem.1l1emagnefzationinboththe℃andFCprodessesissm:wglydependent onthepinningpotentialandthecfiti二alo.!rrentdensity.hutweaklydependentonthetE∞mTamre-sweepingrate.。
“电磁场与电磁波”和“微波技术”课内实验大纲及实验指导书唐万春,车文荃编制陈如山审定南京理工大学通信工程系2006年12月目录1.“电磁场与电磁波”课内实验大纲2.“电磁场与电磁波”课内实验指导说明书实验一电磁波参量的测定实验二电磁波的极化3.“微波技术”课内实验大纲4.“微波技术”课内实验指导说明书实验一传输线的工作状态及驻波比测量实验二微波网络散射参量测试5.“电磁场与电磁波”和“微波技术”课内实验评分标准南京理工大学实验教学大纲课程名称:电磁场与电磁波开课实验室:电磁场与微波技术实验室执笔人:唐万春审定人:陈如山修(制)订日期: 2005年4月*由学校出版、印刷的实验教材(或指导书),统一写作“南京理工大学出版”。
“电磁场与电磁波”课内实验指导书唐万春编写南京理工大学通信工程系二00六年十二月实验一电磁波参量的测定实验1.实验目的a)观察电磁波的传播特性。
b)通过测定自由空间中电磁波的波长,来确定电磁波传播的相位常数k和传播速度v。
c)了解用相干波的原理测量波长的方法。
2.实验内容a)了解并熟悉电磁波综合测试仪的工作特点、线路结构、使用方法。
b)测量信号源的工作波长(或频率)。
3.实验原理与说明a)所使用的实验仪器分度转台晶体检波器可变衰减器喇叭天线反射板固态信号源微安表实验仪器布置图如下:体检波器图1 实验仪器布置图参阅图1。
固态信号源所产生的信号经可变衰减器至矩形喇叭天线,由喇叭天线辐射出去,在接收端用矩形喇叭天线接收,接收到的信号经晶体检波器后通过微安表指示。
b) 原理本实验利用相干波原理,通过测得的电磁波的波长,再由关系式2,k v f kπωλλ===得到电磁波的主要参量k ,v 等。
实验示意图如图2所示。
图中0r P 、1r P 、2r P 和3r P 分别表示辐射喇叭、固定反射板、可动反射板和接收喇叭,图中介质板是一23030()mm ⨯的玻璃板,它对电磁波进行反射、折射后,可实现相干波测试。
二维TE波FDTD法PML吸收边界
摘要:由于计算空间有限,一段时间后波面会被反射回来,这是我们不愿看到的。
为了避免反射波的存在,需要在计算边界设置适当的吸收边界。
目前最为有效的吸收边界是PML层。
关键词:二维;TE波;PML层
1引言
时域有限差分(F D T D )方法已逐渐成为解决电磁散射和电磁波传输等问题的有力工具。
研究F D T D 方法的核心问题是寻求一种理想的吸收边界, 使截断面反射最小。
1 9 9 4 年J.P.Berenger 提出了“完全匹配层(P M L ) ”这种新边界完全匹配层技术在计算电磁学领域特别是在FEM中得到广泛的应用是在Sacks995年提出的单轴各向异性介质PLM以后的事,此后的发展证明,作为一种新型的边界截断策略,不管是与FDTD的结合,还是用在FEM中,PLM技术都获得了巨大成功。
2计算模型
有一个长和宽都为ke=80的矩形介质,在它的周围围上厚度为kp=8的PLM层吸收边界(如图一所示),请编写程序验证电磁波在经过PLM层时完全被吸收(即电磁波没有反射)。
图一模型图
图二左为程序运行结果Ey图右为Hzy图
3结论:本文的算例证明了PLM作为一种边界处理的有效性,PLM技术也在有关
键参数的选取优化问题,这需要研究大量的算例来推出其规律,因而还有很多的工作有待研究。
参考文献:
[1]z S Sacks,D M Kingsland,R Lee and Jin-Fa Lee.A
perfectly matched anisotropic absorber for use as an
absorbing boundary condition[J].IEEE Trans.Antennas
Propagat.,1995,43(12):1460~1463.
[2]A Mitchell,T Aberle,D M Kokotoff and M WAustim
An anisotropic PML for use with biaxial media
D].IEEE Trans.Microwave Theory Tech.,1999,
47(3):374~377.
[3]班永灵.填充各向异性介质的背腔式微带天线的矢量有限元法分析[D].北大学硕士论文,2003年5月.
[4]班永灵(1978一),男,河南人,电子科技大学博士研究生,主要研究方向有限元法(特剐是高阶矢量有限元)、快速扫频技术以及并行计算。
[5]周乐柱(1944一),男,贵州人,北京大学教授、博士生导师,CIE会士,IEEE 高级会员,主要研究方向:计算电磁学及其应用(散射、天线、微波器件),通信中的电磁场问题等。
[6]聂在平(1946一),男,陕西人,电子科技大学教授,博士生导师,现任电子科技大学副校长,主要研究方向:计算电磁学、电磁散射与逆散射、非均匀
介质中的场与渡、移动通信中智能天线技术等。
附录:二维TE波FDTD法PML吸收边界程序
clc,clear
ke=80;
kp=8;
ex=ones(ke+2*kp+1,ke+2*kp+1);
e=ones(ke+2*kp+1,ke+2*kp+1);
ey=ones(ke+2*kp+1,ke+2*kp+1);
hz=ones(ke+2*kp+1,ke+2*kp+1);
hzx=ones(ke+2*kp+1,ke+2*kp+1);
hzy=ones(ke+2*kp+1,ke+2*kp+1);
kc=(ke+2*kp)/2;
e0=8.85*10-12;
mu0=4*3.14*10-7;
pi=3.14;
c0=3e8;
f=1e9;
lambda=c0/f;
dx=lambda/10;
dy=lambda/10;
dt=(dx/2)/c0;
sinxm=0.333;
sinmxm=sigxm*(mu0/e0);
nn=3;
T=0;
for n=1:1000;
T=T+1;
hz=hzx+hzy;
e=ex+ey;
hz(kc,kc)=4*sin(2*pi*f*T*dt);
for i=1:(2*kp+ke+1);
for j=(kp+1):(ke+kp+1);
ex(i,j)=ex(i,j)+dt/(e0*dx)*(hz(i,j)-hz(i,j-1));
end
end
for i=1:(ke+2*kp+1);
for j=2:kp;
sinx=sinxm*((kp+2-j)/kp)^nn*2*e0/dt;
ex(i,j)=ex(i,j)*exp(-sigx*dt/e0)+(1-exp(-sigx*dt/e0))/(dy*sigx)*sqrt(
e0/mu0)*(hz(i,j)-hz(i,j-1));
end
end
for i=1:(ke+2*kp+1);
for j=(ke+kp+2):(ke+2*kp+1);
sinx=sinxm*((j-ke-kp)/kp)^nn*2*e0/dt;
ex(i,j)=ex(i,j)*exp(-sigx*dt/e0)+(1-exp(-sigx*dt/e0))/(dy*sigx)*sqrt( e0/mu0)*(hz(i,j)-hz(i,j-1));
end
end
for i=(kp+1):(ke+kp+1);
for j=1:(2*kp+ke+1);
ey(i,j)=ey(i,j)-dt/(e0*dx)*(hz(i,j)-hz(i-1,j));
end
end
for i=2:kp;
for j=1:(ke+2*kp+1);
sinx=sinxm*((kp+1-i)/kp)^nn*2*e0/dt;
ey(i,j)=ey(i,j)*exp(-sigx*dt/e0)-(1-exp(-sigx*dt/e0))/(dx*sigx)*sqrt( e0/mu0)*(hz(i,j)-hz(i-1,j));
end
end
for i=(ke+kp+2):(ke+2*kp+1);
for j=1:(ke+2*kp+1);
sinx=sinxm*((i-ke-kp-1)/kp)^nn*2*e0/dt;
ey(i,j)=ey(i,j)*exp(-sinx*dt/e0)-(1-exp(-sigx*dt/e0))/(dx*sigx)*sqrt( e0/mu0)*(hz(i,j)-hz(i-1,j));
end
end
%%%%%%%%%%%%%%%磁场循环开始
for i=1:(ke+2*kp+1);
for j=1:kp;
hzy(i,j)=hzy(i,j)*exp(-sigmx*dt/mu0)+(1-exp(-sigmx*dt/mu0))/(dx*sigmx )*sqrt(mu0/e0)*(ex(i,j+1)-ex(i,j));
end
end
for i=1:(2*kp+ke+1);
for j=(kp+1):(kp+ke);
hzy(i,j)=hzy(i,j)+dt/(mu0*dy)*(ex(i,j+1)-ex(i,j));
end
end
for i=1:(ke+2*kp+1);
for j=(kp+ke+1):(ke+2*kp);
sigmx=sigmxm*((j-kp-ke)/kp)^nn*2*mu0/dt;
hzy(i,j)=hzy(i,j)*exp(-sigmx*dt/mu0)+(1-exp(-sigmx*dt/mu0))/(dx*sigmx )*sqrt(mu0/e0)*(ex(i,j+1)-ex(i,j));
end
end
for i=1:kp;
for j=1:(ke+2*kp+1);
sigmx=sigmxm*((kp-i+1)/kp)^nn*2*mu0/dt;
hzx(i,j)=hzx(i,j)*exp(-sigmx*dt/mu0)-(1-exp(-sigmx*dt/mu0))/(dx*sigmx )*sqrt(mu0/e0)*(ey(i+1,j)-ey(i,j));
end
end
for i=(kp+1):(kp+ke);
for j=1:(2*kp+ke+1);
hzx(i,j)=hzx(i,j)+dt/(mu0*dy)*(ey(i+1,j)-ey(i,j));
end
end
for i=(ke+kp+1):(ke+2*kp);
for j=1:(ke+2*kp+1);
hzx(i,j)=hzx(i,j)*exp(-sigmx*dt/mu0)-(1-exp(-sigmx*dt/mu0))/(dx*sigmx )*sqrt(mu0/e0)*(ey(i+1,j)-ey(i,j));
end
end
ey(1,:)=0;
ex(:,1)=0;
hzx(2*kp+ke+1,:)=0;
hzy(:,2*kp+ke+1)=0;
plot(hz(:,kc))
axis([46,ke+2*kp,45,50])
figure(1)
subplot(1,2,1)
mesh(ey)
drawnow
subplot(1,2,2)
mesh(hzy)
drawnow
figure(2)
mesh(hz)
drawnow
end
hz
%pcolor(hz)
figure(3)
contour3(hz,5)
axis([kp,ke+kp,kp,kp+ke])
shading interp
%print-djpeg D2pml相位图
%print-djepg D2pml 等高线图。