二维TM波讨论平面波源(使用直接算方法)的加入
平面波公式
平面波公式
平面波公式是物理学中的一个概念，表达了质点在空间中随时间的变化。
平面波公式可以被用来描述振动，波浪等物理变化情况。
它为我们探索宇宙奥秘提供了重要的指导意义。
平面波公式的本质
平面波公式是由线性液动力学理论建立起来的，这里只考虑了液体在物体表面上的摩擦力，忽略了液体内部的变化现象，所以它可以简化为形如：
F=Acos（k（x-ut））
其中A是振幅，k是波数，u是速度，t是时间，x是空间变量。
该公式描述的是质点随时间的变化受空间变量影响的力学运动规律。
应用
以上的公式可以用来解释许多现象，比如波浪在大海里的传播；地震波在地底传播；声音在空气中传播等等。
平面波公式具有广泛的应用，可以应用于地面建筑和海洋结构物的动力设计，军事和航空航天技术中，将其用于测量和导航，以及电子技术中，将其用于电磁波的控制。
另外，平面波公式还可以用于液力学的研究。
液力学是研究液体力学运动的理论，其关键观点是，当液体运动时，液体中心会产生一个力，并且受空间变量影响，利用平面波公式，可以详细描述液体运动情况，为液力学学习者提供了重要的参考。
最后，平面波公式也可以用来模拟物理实验，比如微波加热实验
中，可以利用平面波公式来模拟微波在物质表面的传播，这从很大程度上提高了实验的可信度和准确性。
总结
总而言之，平面波公式是一个重要的概念，它是描述质点在不同空间位置随时间的变化的表达式，它为我们解开物理世界的奥秘提供了有价值的指导。
平面波公式的本质是线性液动力学，其应用非常广泛，可以用在液力学，航空航天，电子技术，建筑和海洋结构物等，还可以用来模拟实验。
二维波动方程
二维波动方程二维波动方程是物理学和数学中一个重要的概念。
它是一类非线性的偏微分方程，用来描述物理系统中的波动现象。
二维波动方程的应用非常广泛，可以用来研究电磁波、热传导、流体力学等领域中的波动现象。
其中最常用的二维波动方程之一是Korteweg-de Vries方程(KdV方程)，它用来描述一维流体中的纵波。
它可以用来描述潮汐、河流涨落、海浪等现象。
它的形式是：u_t + uu_x + u_{xxx} = 0。
其中u(x,t)表示流体的速度场，x和t分别表示空间坐标和时间。
另一个重要的二维波动方程是Sine-Gordon方程，它描述了非线性玄学领域中的波动现象，例如在线性电力线的电磁波传播，在非线性声学领域中的声波传播。
这个方程的形式是：u_{tt}-u_{xx}+sin(u)=0。
还有一种重要的二维波动方程是Nonlinear Schrödinger方程(NLS方程) ,它是研究光学波动、超流体波动等领域中的重要方程。
NLS方程的形式是：i∂u/∂t+1/2∂²u/∂x²+|u|²u=0。
在这些方程中，通常都采用数值方法和解析方法来解决它们。
数值方法包括有限差分法、有限元法、有限体积法等,而解析方法则有传统的小扰动分析和现代的非线性分析等。
总的来说，二维波动方程是物理学和数学领域中重要的概念，它们被广泛应用于研究物理系统中的波动现象。
使用数值方法和解析方法来解决这些方程是目前常用的方法。
综上所述，二维波动方程是一类重要的非线性偏微分方程，应用广泛，用来描述物理系统中的波动现象。
主要有Korteweg-de Vries方程，Sine-Gordon方程，Nonlinear Schrödinger方程等。
在研究中采用数值方法和解析方法来解决这些方程是常见的方法。
TM波FDTD讨论平面波源的加入
! TM波FDTD讨论平面波源的加入module data_moduleimplicit noneinteger,parameter::nx0=0,nx1=360,ny0=0,ny1=360,nz0=-100,nz1=1200integer,parameter::nxl1=nx0+80,nxl2=nx1-80,nyl1=ny0+80,nyl2=ny1-80 !连接边界real,parameter::f=2.0e8,c=3.0e8,delt=0.0177,deltt=delt/6.0e8,eps0=8.85e-12,miu0=1.2566e-6,pi= 3.14159real,parameter::w=2*pi*f,s=-0.477369real,parameter::p=-1.0/3.0,q=-miu0*c/6,r=-miu0*c/2,p1=1/(2*miu0*c),p2=1/(2*eps0*c)real,parameter::tal=2e-9,t0=0.8*tal,fai=pi/3.0real cez,chx,chyinteger,parameter::nt=2000,m0=200integer ncomplex Ez3(nx0:nx1,ny0:ny1)real Ez4(nx0:nx1,ny0:ny1),Ez2(nx0:nx1,ny0:ny1) !记录幅值提取时的实部和虚部real sita(nx0:nx1,ny0:ny1),Ez0(nx0:nx1,ny0:ny1) !记录幅值和相位real Ez(nx0:nx1,ny0:ny1),Hx(nx0:nx1,ny0:ny1),Hy(nx0:nx1,ny0:ny1),Ez1(nx0:nx1,ny0:ny1) real Ei0(nz0:nz1),Hi0(nz0:nz1),Ei1(nz0:nz1)real Ezi(nx0:nx1,ny0:ny1),Hxi(nx0:nx1,ny0:ny1),Hyi(nx0:nx1,ny0:ny1)end module data_module!/////////////////////////////////////////////////////////////////////////////////////////////////subroutine inc()use data_moduleimplicit noneinteger i,j,kreal t,dt=n*delttEi1=Ei0do k=nz0,nz1-1Hi0(k)=Hi0(k)-p1*(Ei1(k+1)-Ei1(k))end do!Ezido i=nxl1,nxl2do j=nyl1,nyl2d=real(i-nxl1)*cos(fai)+real(j-nyl1)*sin(fai)Ezi(i,j)=(d-int(d))*Ei0(m0+int(d)+1)+(1-(d-int(d)))*Ei0(m0+int(d))end doend dodo k=nz0+1,nz1-1Ei0(k)=Ei1(k)-p2*(Hi0(k)-Hi0(k-1)) ! 入射波的场量end doEi0(m0-2)=sin(w*t) !exp(-(4*pi*(t-t0)**2/Tal/tal))k=nz1Ei0(k)=Ei1(k-1)-1.0/3.0*(Ei0(k-1)-Ei1(k)) !入射波右吸收边界K=nz0Ei0(k)=Ei1(k+1)-1.0/3.0*(Ei0(k+1)-Ei1(k)) !入射波左吸收边界!Hxij=nyl1do i=nxl1,nxl2d=real(i-nxl1)*cos(fai)+(j-0.5-nyl1)*sin(fai)+0.5Hxi(i,j-1)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*sin(fai) end doj=nyl2do i=nxl1,nxl2d=real(i-nxl1)*cos(fai)+(j+0.5-nyl1)*sin(fai)+0.5Hxi(i,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*sin(fai) end do!Hyii=nxl1do j=nyl1,nyl2d=(i-0.5-nyl1)*cos(fai)+real(j-nyl1)*sin(fai)+0.5Hyi(i-1,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*(-cos(fai)) enddoi=nxl2do j=nyl1,nyl2d=(i+0.5-nyl1)*cos(fai)+real(j-nyl1)*sin(fai)+0.5Hyi(i,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*(-cos(fai)) enddowrite(13,*)n,Ei0(0),Ei0(20),Ei0(30),Ei0(100),sin(w*t)end subroutine inc!/////////////////////////////////////////////////////////////////////////////////////////////////program mainuse data_moduleimplicit noneinteger i,j,kreal topen(11,file='A.txt')open(12,file='fai.txt')open(13,file='Ez.txt',recl=300)Ez=0;Hx=0;Hy=0;Ez1=0Ei0=0;Hi0=0Ez4=0;Ez2=0;Ez3=cmplx(0.0,0.0)cez=1.0;chx=sin(fai);chy=-cos(fai)do n=0,ntcall inc()call FDTD()print *,'step',n,Ez(nx1/2,ny1/2)!write(13,*)n,Ei0(0),Ei0(20),Ei0(30),Ei0(100),sin(w*t)Ez1=Ezend doEz4=Ezdo n=nt+1,nt+42call inc()call FDTD()Ez1=Ezend doEz2=Ezcall AP()do j=ny0,ny1do i=nx0,nx1write(11,*)i,j,Ez0(i,j)write(12,*)i,j,sita(i,j)end doend doclose(11)close(12)close(13)end program main!///////////////////////////////////////////////////////////////////////////////////////////////// subroutine FDTD()use data_moduleimplicit noneinteger i,jcall cal_H()call cal_E()end subroutine FDTD!///////////////////////////////////////////////////////////////////////////////////////////////// subroutine cal_H()use data_moduleimplicit noneinteger i,jreal,external::Eido j=ny0,ny1-1Hx(i,j)=Hx(i,j)-p1*(Ez(i,j+1)-Ez(i,j))end doend dodo j=ny0,ny1do i=nx0,nx1-1Hy(i,j)=Hy(i,j)+p1*(Ez(i+1,j)-Ez(i,j))end doend do!连接边界i=nxl1-1do j=nyl1,nyl2Hy(i,j)=Hy(i,j)-p1*Ezi(i+1,j)end doi=nxl2do j=nyl1,nyl2Hy(i,j)=Hy(i,j)+p1*Ezi(i,j)end doj=nyl1-1do i=nxl1,nxl2Hx(i,j)=Hx(i,j)+p1*Ezi(i,j+1)end doj=nyl2do i=nxl1,nxl2Hx(i,j)=Hx(i,j)-p1*Ezi(i,j)end doend subroutine cal_H!///////////////////////////////////////////////////////////////////////////////////////////////// subroutine cal_E()use data_moduleimplicit noneinteger i,jreal,external::Hido j=ny0+1,ny1-1do i=nx0+1,nx1-1Ez(i,j)=Ez(i,j)+p2*(Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1))end doend do!Ez(nx1/2,ny1/2)=Ez(nx1/2,ny1/2)-p2/delt*sin(w*(n+0.5)*deltt) call mur()! 连接边界!四条边do j=nyl1+1,nyl2-1Ez(i,j)=Ez(i,j)-p2*Hyi(i-1,j)end doi=nxl2do j=nyl1+1,nyl2-1Ez(i,j)=Ez(i,j)+p2*Hyi(i,j)end doj=nyl1do i=nxl1+1,nxl2-1Ez(i,j)=Ez(i,j)+p2*Hxi(i,j-1)end doj=nyl2do i=nxl1+1,nxl2-1Ez(i,j)=Ez(i,j)-p2*Hxi(i,j)end do!四个角i=nxl1j=nyl1Ez(i,j)=Ez(i,j)+p2*(Hxi(i,j-1)-Hyi(i-1,j))i=nxl2Ez(i,j)=Ez(i,j)+p2*(Hxi(i,j-1)+Hyi(i,j))j=nyl2Ez(i,j)=Ez(i,j)+p2*(-Hxi(i,j)+Hyi(i,j))i=nxl1Ez(i,j)=Ez(i,j)+p2*(-Hxi(i,j)-Hyi(i-1,j))end subroutine cal_E!/////////////////////////////////////////////////////////////////////////////////////////////////subroutine mur()use data_moduleimplicit noneinteger i,j!四条边i=nx0do j=ny0+1,ny1-1Ez(i,j)=Ez1(i+1,j)+p*(Ez(i+1,j)-Ez1(i,j))+q*(Hx(i+1,j)+Hx(i,j)-Hx(i+1,j-1)-Hx(i,j-1)) end doi=nx1-1do j=ny0+1,ny1-1Ez(i+1,j)=Ez1(i,j)+p*(Ez(i,j)-Ez1(i+1,j))+q*(Hx(i+1,j)+Hx(i,j)-Hx(i+1,j-1)-Hx(i,j-1)) end doj=ny0do i=nx0+1,nx1-1Ez(i,j)=Ez1(i,j+1)+p*(Ez(i,j+1)-Ez1(i,j))-q*(Hy(i,j+1)+Hy(i,j)-Hy(i-1,j+1)-Hy(i-1,j)) end doj=ny1-1do i=nx0+1,nx1-1Ez(i,j+1)=Ez1(i,j)+p*(Ez(i,j)-Ez1(i,j+1))-q*(Hy(i,j+1)+Hy(i,j)-Hy(i-1,j+1)-Hy(i-1,j)) end do!角点i=nx0j=ny0Ez(i,j)=Ez1(i+1,j+1)+s*(Ez(i+1,j+1)-Ez1(i,j))i=nx1Ez(i,j)=Ez1(i-1,j+1)+s*(Ez(i-1,j+1)-Ez1(i,j))j=ny1Ez(i,j)=Ez1(i-1,j-1)+s*(Ez(i-1,j-1)-Ez1(i,j))i=nx0Ez(i,j)=Ez1(i+1,j-1)+s*(Ez(i+1,j-1)-Ez1(i,j))end subroutine mur!/////////////////////////////////////////////////////////////////////////////////////////////////////// subroutine AP()use data_moduleimplicit noneinteger i,jreal mreal tmp1,tmp2,tmp3complex tmptmp=cmplx(Ez2(nx1/2,ny1/2),Ez4(nx1/2,ny1/2))do j=ny0,ny1do i=nx0,nx1Ez3(i,j)=cmplx(Ez2(i,j),Ez4(i,j))Ez0(i,j)=sqrt(Ez2(i,j)**2+Ez4(i,j)**2)m=Imag(Ez3(i,j)/tmp)/real(Ez3(i,j)/tmp)if(real(Ez3(i,j)/tmp)==0 .and. Imag(Ez3(i,j)/tmp)>0) thensita(i,j)=pi/2elseif(real(Ez3(i,j)/tmp)==0 .and. Imag(Ez3(i,j)/tmp)<0) thensita(i,j)=-pi/2elsesita(I,j)=atan(m)end ifif (abs(sita(i,j))<=1e-6) thensita(i,j)=0end ifend doend doend subroutine AP!///////////////////////////////////////////////////////////////////////////////////////////////////////。
平面波的特性与波长的计算
平面波的特性与波长的计算波动现象是自然界中普遍存在的一种物理现象，而波长是用来描述波的特性的定量指标之一。
在物理学中，平面波是一种理想化的波模型，其特性与波长之间存在一定的关系。
本文将探讨平面波的特性以及如何计算波长。
1. 平面波的特性平面波是一种在无限远处传播、波前呈平面的波动。
与球面波和柱面波不同，平面波的波阵面是一个平面，不存在波源和波纹的聚焦或发散现象。
平面波特性的关键在于其波矢和传播方向。
平面波的波矢k定义为波矢的模长与传播方向的乘积，即k = |k| * n。
其中，|k|表示波矢的大小，n表示波的传播方向，即单位矢量。
平面波的表达式可以表示为Ae^(ik·r)，其中A为振幅，k为波矢，r为位置矢量。
平面波的传播速度与波矢的大小成正比，即v = ω/|k|，其中v表示传播速度，ω表示角频率。
2. 波长的计算波长是用来描述波动的空间周期性的物理量，指在一个完整波动周期内，波动传播的距离。
对于平面波来说，波长的计算可以通过波矢的大小求倒数得到。
波矢的模长|k|与波长λ之间存在以下关系：|k| = 2π/λ。
由此可得，波长λ = 2π/|k|。
这个公式说明了波矢的大小与波长之间的倒数关系，即波矢越大，波长越短，反之亦然。
以光波为例，光是一种电磁波，其波长范围很宽，从红外到紫外都有不同波长的光。
对于可见光来说，其波长范围约为380nm到780nm之间。
3. 示例计算现在我们以一个具体的例子来计算平面波的波长。
假设我们有一个平面波，其振幅A为2，波矢k为1.5。
根据前述公式，波长λ = 2π/|k| = 2π/1.5 ≈ 4.188。
根据计算结果可知，该平面波的波长约为4.188个长度单位，可以是米、厘米、纳米等等，具体的单位需根据具体的物理系统来确定。
4. 总结平面波是一种在无限远处传播且波前呈平面的波动，其特性由波矢和传播方向决定。
波长是描述波动空间周期性的物理量，对于平面波来说，波长可以通过波矢的大小求倒数得到。
PML吸收边界条件下TM波的时域有限差分法分析
会带来较强的数值反射，选取线性分布可以改善结果。通常，选取 n=2 或 3，可进一步改进 吸收效果。
4. TM 波实例分析
4. 1 问题提出
二 维 TM 波 在 自 由 空 间 传 播 ， 采 用 PML 吸 收 边 界 ， 激 励 源 使 用 脉 冲 源
为了让实际的计算空间与实际的无限的物理空间等效须对有限空间的周围边界做特殊处理使得向边界面行进的波在边界处保持向外行进特点无反射现象即好像被吸收了一样并且不会使内部空间的场产生畸变具有这种功能的边界条件称为吸收边界条件吸收边界条件处理成为了fdtd的一大重点和发展方向pml是近年来发展起来的吸收边界条件45二维情况下pml吸收层区分为图3所示的边和角顶两种区域
∆t ε (m)
− σ (m) 2
+ σ (m)
=
1− 1+
σ (m)∆t 2ε (m) σ (m)∆t
(5)
∆t 2
2ε (m)
∆t
CB(m)
=
ε(m)
1 +σ(m)
=
ε(m) 1+σ(m)∆t
(6)
∆t 2
2ε(m)
CP(m)
=
µ(m)
∆t µ(m)
− σm(m) 2
+ σm(m)
1− σm(m)∆t
外行波在四边的 PML 界面上是无反射的。
图 3 完全匹配层的设置
在四个角顶区域，介质参数为 PML( σx , σmx , σy , σmy ),而于之相邻的边上的参数为 PML(σx ,σmx ,0,0,)和 PML(0,0，σy ,σmy )。根据上面的讨论，满足匹配条件时，在角顶边和
二维超声速普朗特-迈耶系数波流场的数值解
课程设计题目：二维超声速普朗特-迈耶系数波流场的数值解学院：飞行器工程学院专业名称：飞行器设计与工程班级学号： 07034211学生姓名：李桂平指导教师：刘勇二O一O年十一月第一部分1．物理问题简介：普朗特——迈耶稀疏波的解析解图-1中，超声速流围绕着一个尖的扩张角膨胀，无数个无限弱的马赫波组成了稀疏波，在尖角处展开成扇形。
扇形稀疏波的波头与來流方向的夹角1，而2是其波尾与下游方向的夹角。
1和2称为马赫角，定义为：和Ma1和Ma2分别为上下游的马赫数。
通过稀疏波的流动是等熵流动。
当流体通过稀疏波后，马赫数增加，压力、温度和密度降低；图-1中标明了这些变化趋势。
在中心稀疏波前的流动是均匀的，马赫数为Ma1，而且流动平行与波前的壁面。
稀疏波后的流动是均匀的，马赫数为Ma2，并且平行于下游的壁面。
在稀疏波内，流动参数光滑变化，流线弯曲，如图-1所示。
稀疏波内的流动是二维的，唯一的例外是折角的定点，它是一个奇点，壁面流线的方向在此处有一个突然的变化，而且此处的流动参数也是不连续的。
这个奇点对流动的数值解会产生影响。
给定超声速来流条件和拐角处的偏转角，下游参数是唯一确定的。
对于完全气体，在稀疏波后的流动有精确的解析解，下面给出这个解。
流过中心稀疏波的流动，其解析解取决于简单的关系式……1 式中，f是普朗特——迈耶函数；是流动偏转角。
对于完全气体，普朗特——迈耶函数是Ma和γ的函数，定义为……2 解析解中如下依次得到。
对给定的Ma1，从式(2)计算函数f1。
然后，对给定的偏转角θ，从式(1)得到f2。
用这样得到f2的值，通过求解式(2)求出Ma2。
式(2)是关于Ma2的隐式关系式，需要用试凑法求解。
波后的压力温度和密度都可以由等熵流动关系式： (3) (4)和状态方程： (5)得到。
借助是式(1)～式(5)，中心稀疏波后的流动就完全确定了。
2.问题的提法考虑图-2所示的物理平面。
来流马赫数为2，来流的压力、密度、温度分别为：1.01x105N/m2、1.23kg/m3、286.1K。
平面波法计算能带的基本原理
平面波法计算能带的基本原理能带理论是描述固体中电子能量分布的理论，它对于理解固体的电子性质和材料的导电、绝缘以及半导体特性具有重要意义。
平面波法是计算能带的一种常用方法，本文将介绍平面波法计算能带的基本原理。
平面波法是一种基于量子力学的计算方法，它假设电子波函数能够用平面波表示。
平面波是一种特殊的波函数形式，它具有平直的波前和无限延伸的波长。
根据波动方程，平面波的波函数可以表示为Ae^(ik•r)的形式，其中A是振幅，k是波矢，r是位置矢量。
在能带计算中，我们通常采用周期性边界条件，即假设固体是无限重复的晶体结构。
根据布洛赫定理，电子的波函数可以表示为平面波和周期函数的乘积形式。
这样，我们可以通过计算晶体中的平面波的行为来研究电子的能量分布。
平面波法的基本思路是，首先选择一组适当的平面波作为基函数，然后将它们用作电子波函数的展开系数。
通常，我们选择的平面波基函数是满足布洛赫定理的，即满足晶格平移对称性的平面波。
这样，我们就可以通过求解薛定谔方程来确定电子的能量和波函数。
在平面波法中，薛定谔方程可以写为哈密顿算符作用在波函数上等于能量的形式。
由于平面波是哈密顿算符的本征函数，所以我们可以将波函数用平面波展开，并代入薛定谔方程中，从而得到一个简化的形式。
通过对展开系数的求解，我们可以得到电子的能量和波函数。
在实际计算中，我们通常采用周期性边界条件，将晶体划分为一个个无限重复的晶胞。
每个晶胞内的电子波函数可以用平面波展开，而不同晶胞之间的平面波之间存在相位差。
通过求解薛定谔方程，我们可以得到不同晶胞中的平面波的展开系数，从而确定电子的能量和波函数。
平面波法计算能带的基本原理就是通过选择适当的平面波基函数，并将其用作电子波函数的展开系数，从而求解薛定谔方程得到能量和波函数。
这种方法在实际应用中已经取得了很多成功，并对我们理解固体的电子性质提供了重要的帮助。
平面波法是一种常用的计算能带的方法，它基于平面波的特点和布洛赫定理，通过求解薛定谔方程来确定电子的能量和波函数。
二维核磁共振谱全解
O
2 C1
C
O CH2
CH3
H
5
6
5
46
1/2JCH F1 (HZ)
-100
0
100
δ C 150
100
50 F2
反式丙烯酸酯的13C,1H J分解2D谱。 F2轴为13C的化学位移。F1轴为1JCH偶合的多重峰。
21
10
7
9
5,6
16
2
8
O
3
5 4
由于DEPT等测定碳 原子级数的方法能代 替异核J 谱，且检测 速度快，操作方便， 因此异核J 谱较少应 用。
15
在同核1H,1H －2D J分解谱中，被测定的核为 1H核。
16
3
2
5
4
6
H
3
C CH3
4
O
C1
5
6
C
2
O CH2 CH3
H
F1
-10
J0 (HZ)
10
δH 6
4
2 F2
反式丙烯酸乙酯的1H－1H同核二维J分解谱
17
18
化学等价磁 不等价
化学等价磁 不等价
2位的a-H和e-H相互偶合， 裂分为双峰后进一步被3位 氢裂分而显示六重峰。同 样，3-H的a、e质子首先裂 分为双峰，再被邻位2-H分 裂为六重峰，和4-H质子进 一步裂分为十八重峰，但 由于有些峰相互重叠，在 2DJ分解谱中只显示9个点。 5位甲基没有受到偶合，因 此只在F1=0轴上显示单峰。
C6’,H6’ C6,H6
C1’,H1’ 36
F2
7,4
CH3CO
37
（二）1H检测的异核多量子相关谱（HMQC）
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
! TM波FDTD讨论平面波源的加入module data_moduleimplicit noneinteger,parameter::nx0=0,nx1=360,ny0=0,ny1=360,nz0=-100,nz1=1200integer,parameter::nxl1=nx0+80,nxl2=nx1-80,nyl1=ny0+80,nyl2=ny1-80 !连接边界real,parameter::f=2.0e8,c=3.0e8,delt=0.0177,deltt=delt/6.0e8,eps0=8.85e-12,miu0=1.2566e-6,pi= 3.14159real,parameter::w=2*pi*f,s=-0.477369real,parameter::p=-1.0/3.0,q=-miu0*c/6,r=-miu0*c/2,p1=1/(2*miu0*c),p2=1/(2*eps0*c)real,parameter::tal=2e-9,t0=0.8*tal,fai=pi/3.0real cez,chx,chyinteger,parameter::nt=2000,m0=200integer ncomplex Ez3(nx0:nx1,ny0:ny1)real Ez4(nx0:nx1,ny0:ny1),Ez2(nx0:nx1,ny0:ny1) !记录幅值提取时的实部和虚部real sita(nx0:nx1,ny0:ny1),Ez0(nx0:nx1,ny0:ny1) !记录幅值和相位real Ez(nx0:nx1,ny0:ny1),Hx(nx0:nx1,ny0:ny1),Hy(nx0:nx1,ny0:ny1),Ez1(nx0:nx1,ny0:ny1) real Ei0(nz0:nz1),Hi0(nz0:nz1),Ei1(nz0:nz1)real Ezi(nx0:nx1,ny0:ny1),Hxi(nx0:nx1,ny0:ny1),Hyi(nx0:nx1,ny0:ny1)end module data_module!/////////////////////////////////////////////////////////////////////////////////////////////////subroutine inc()use data_moduleimplicit noneinteger i,j,kreal t,dt=n*delttEi1=Ei0do k=nz0,nz1-1Hi0(k)=Hi0(k)-p1*(Ei1(k+1)-Ei1(k))end do!Ezido i=nxl1,nxl2do j=nyl1,nyl2d=real(i-nxl1)*cos(fai)+real(j-nyl1)*sin(fai)Ezi(i,j)=(d-int(d))*Ei0(m0+int(d)+1)+(1-(d-int(d)))*Ei0(m0+int(d))end doend dodo k=nz0+1,nz1-1Ei0(k)=Ei1(k)-p2*(Hi0(k)-Hi0(k-1)) ! 入射波的场量end doEi0(m0-2)=sin(w*t) !exp(-(4*pi*(t-t0)**2/Tal/tal))k=nz1Ei0(k)=Ei1(k-1)-1.0/3.0*(Ei0(k-1)-Ei1(k)) !入射波右吸收边界K=nz0Ei0(k)=Ei1(k+1)-1.0/3.0*(Ei0(k+1)-Ei1(k)) !入射波左吸收边界!Hxij=nyl1do i=nxl1,nxl2d=real(i-nxl1)*cos(fai)+(j-0.5-nyl1)*sin(fai)+0.5Hxi(i,j-1)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*sin(fai) end doj=nyl2do i=nxl1,nxl2d=real(i-nxl1)*cos(fai)+(j+0.5-nyl1)*sin(fai)+0.5Hxi(i,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*sin(fai) end do!Hyii=nxl1do j=nyl1,nyl2d=(i-0.5-nyl1)*cos(fai)+real(j-nyl1)*sin(fai)+0.5Hyi(i-1,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*(-cos(fai)) enddoi=nxl2do j=nyl1,nyl2d=(i+0.5-nyl1)*cos(fai)+real(j-nyl1)*sin(fai)+0.5Hyi(i,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*(-cos(fai)) enddowrite(13,*)n,Ei0(0),Ei0(20),Ei0(30),Ei0(100),sin(w*t)end subroutine inc!/////////////////////////////////////////////////////////////////////////////////////////////////program mainuse data_moduleimplicit noneinteger i,j,kreal topen(11,file='A.txt')open(12,file='fai.txt')open(13,file='Ez.txt',recl=300)Ez=0;Hx=0;Hy=0;Ez1=0Ei0=0;Hi0=0Ez4=0;Ez2=0;Ez3=cmplx(0.0,0.0)cez=1.0;chx=sin(fai);chy=-cos(fai)do n=0,ntcall inc()call FDTD()print *,'step',n,Ez(nx1/2,ny1/2)!write(13,*)n,Ei0(0),Ei0(20),Ei0(30),Ei0(100),sin(w*t)Ez1=Ezend doEz4=Ezdo n=nt+1,nt+42call inc()call FDTD()Ez1=Ezend doEz2=Ezcall AP()do j=ny0,ny1do i=nx0,nx1write(11,*)i,j,Ez0(i,j)write(12,*)i,j,sita(i,j)end doend doclose(11)close(12)close(13)end program main!///////////////////////////////////////////////////////////////////////////////////////////////// subroutine FDTD()use data_moduleimplicit noneinteger i,jcall cal_H()call cal_E()end subroutine FDTD!///////////////////////////////////////////////////////////////////////////////////////////////// subroutine cal_H()use data_moduleimplicit noneinteger i,jreal,external::Eido j=ny0,ny1-1Hx(i,j)=Hx(i,j)-p1*(Ez(i,j+1)-Ez(i,j))end doend dodo j=ny0,ny1do i=nx0,nx1-1Hy(i,j)=Hy(i,j)+p1*(Ez(i+1,j)-Ez(i,j))end doend do!连接边界i=nxl1-1do j=nyl1,nyl2Hy(i,j)=Hy(i,j)-p1*Ezi(i+1,j)end doi=nxl2do j=nyl1,nyl2Hy(i,j)=Hy(i,j)+p1*Ezi(i,j)end doj=nyl1-1do i=nxl1,nxl2Hx(i,j)=Hx(i,j)+p1*Ezi(i,j+1)end doj=nyl2do i=nxl1,nxl2Hx(i,j)=Hx(i,j)-p1*Ezi(i,j)end doend subroutine cal_H!///////////////////////////////////////////////////////////////////////////////////////////////// subroutine cal_E()use data_moduleimplicit noneinteger i,jreal,external::Hido j=ny0+1,ny1-1do i=nx0+1,nx1-1Ez(i,j)=Ez(i,j)+p2*(Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1))end doend do!Ez(nx1/2,ny1/2)=Ez(nx1/2,ny1/2)-p2/delt*sin(w*(n+0.5)*deltt) call mur()! 连接边界!四条边do j=nyl1+1,nyl2-1Ez(i,j)=Ez(i,j)-p2*Hyi(i-1,j)end doi=nxl2do j=nyl1+1,nyl2-1Ez(i,j)=Ez(i,j)+p2*Hyi(i,j)end doj=nyl1do i=nxl1+1,nxl2-1Ez(i,j)=Ez(i,j)+p2*Hxi(i,j-1)end doj=nyl2do i=nxl1+1,nxl2-1Ez(i,j)=Ez(i,j)-p2*Hxi(i,j)end do!四个角i=nxl1j=nyl1Ez(i,j)=Ez(i,j)+p2*(Hxi(i,j-1)-Hyi(i-1,j))i=nxl2Ez(i,j)=Ez(i,j)+p2*(Hxi(i,j-1)+Hyi(i,j))j=nyl2Ez(i,j)=Ez(i,j)+p2*(-Hxi(i,j)+Hyi(i,j))i=nxl1Ez(i,j)=Ez(i,j)+p2*(-Hxi(i,j)-Hyi(i-1,j))end subroutine cal_E!/////////////////////////////////////////////////////////////////////////////////////////////////subroutine mur()use data_moduleimplicit noneinteger i,j!四条边i=nx0do j=ny0+1,ny1-1Ez(i,j)=Ez1(i+1,j)+p*(Ez(i+1,j)-Ez1(i,j))+q*(Hx(i+1,j)+Hx(i,j)-Hx(i+1,j-1)-Hx(i,j-1)) end doi=nx1-1do j=ny0+1,ny1-1Ez(i+1,j)=Ez1(i,j)+p*(Ez(i,j)-Ez1(i+1,j))+q*(Hx(i+1,j)+Hx(i,j)-Hx(i+1,j-1)-Hx(i,j-1)) end doj=ny0do i=nx0+1,nx1-1Ez(i,j)=Ez1(i,j+1)+p*(Ez(i,j+1)-Ez1(i,j))-q*(Hy(i,j+1)+Hy(i,j)-Hy(i-1,j+1)-Hy(i-1,j)) end doj=ny1-1do i=nx0+1,nx1-1Ez(i,j+1)=Ez1(i,j)+p*(Ez(i,j)-Ez1(i,j+1))-q*(Hy(i,j+1)+Hy(i,j)-Hy(i-1,j+1)-Hy(i-1,j)) end do!角点i=nx0j=ny0Ez(i,j)=Ez1(i+1,j+1)+s*(Ez(i+1,j+1)-Ez1(i,j))i=nx1Ez(i,j)=Ez1(i-1,j+1)+s*(Ez(i-1,j+1)-Ez1(i,j))j=ny1Ez(i,j)=Ez1(i-1,j-1)+s*(Ez(i-1,j-1)-Ez1(i,j))i=nx0Ez(i,j)=Ez1(i+1,j-1)+s*(Ez(i+1,j-1)-Ez1(i,j))end subroutine mur!/////////////////////////////////////////////////////////////////////////////////////////////////////// subroutine AP()use data_moduleimplicit noneinteger i,jreal mreal tmp1,tmp2,tmp3complex tmptmp=cmplx(Ez2(nx1/2,ny1/2),Ez4(nx1/2,ny1/2))do j=ny0,ny1do i=nx0,nx1Ez3(i,j)=cmplx(Ez2(i,j),Ez4(i,j))Ez0(i,j)=sqrt(Ez2(i,j)**2+Ez4(i,j)**2)m=Imag(Ez3(i,j)/tmp)/real(Ez3(i,j)/tmp)if(real(Ez3(i,j)/tmp)==0 .and. Imag(Ez3(i,j)/tmp)>0) thensita(i,j)=pi/2elseif(real(Ez3(i,j)/tmp)==0 .and. Imag(Ez3(i,j)/tmp)<0) thensita(i,j)=-pi/2elsesita(I,j)=atan(m)end ifif (abs(sita(i,j))<=1e-6) thensita(i,j)=0end ifend doend doend subroutine AP!///////////////////////////////////////////////////////////////////////////////////////////////////////。