integer,parameter::Imin=-75,Imax=75, Jmin=-75, Jmax=75 ! 总边界(吸收边界)
integer,parameter::Itmin=-65,Itmax=65,Jtmin=-65,Jtmax=65 ! 连接边界
integer,parameter::Iomin=-70,Iomax=70,Jomin=-70,Jomax=70 ! 输出边界
real Ez(-160:160,-160:160) ! FDTD 采样点Ez
real Hx(-160:160,-160:159) ! FDTD 采样点Hx
real Hy(-160:159,-160:160) ! FDTD 采样点Hy
real Ein(-800:800) ! 入射波E采样点
real Hin(-800:799) ! 入射波H采样点
integer IncidentStart, IncidentEnd,Isource ! 一维入射波数组的端点和原点位
real EBin(4) ! 入射波吸收边界条件临时变量
real ECB(-150:150,4) ! 连接边界处E入射波分量.
real HCB(-150:150,4) ! 连接边界处H入射波分量.
real EAB(-160:160,4) ! 吸收边界条件临时变量
real EAC(4,2,0:1) ! 吸收边界条件角点临时变量
integer ob(-150:150,-150:150) ! 目标数组,存储FDTD网格所属介质编号
real FE1(0:10,0:10) ! FDTD迭代公式系数
real FE2(0:10,0:10) !
real FH1(0:10,0:10) !
real FH2(0:10,0:10) !
integer TEMflag ! 1 代表 TM 极化; -1 代表 TE极化
real MU0 ! 自由空间电磁参数:磁导率
real EPS0 ! 介电常数
real SMU0 ! 磁导率的平方根
real SEPS0 ! 介电常数的平方根
real PI ! 圆周率
real OMIGA ! 数值角频率
integer WL ! 数值波长
real IncidentAngle ! 入射角86
integer TimeStop ! 总时间步
character*1 OUTflag !输出标志: R为输出场量的实部, I为虚部
character*2 WaveMode
real k,Z,WaveLength !自由空间波数、波阻抗、波长
real T1,T2,T3,T4 !临时存储变量.
complex CT1,CT2,CT3,CT4
integer I,J,N,II,III,IFlag
real Media(1:4,0:10) !介质参数数组
integer MediaNumber !目标所包含介质总数目
complex Jz !输出面上的电流密度
complex Mx !磁流密度
complex My !磁流密度
complex Ctemp,ctemp1 !临时复数变量
complex E(-1250:1250,4) !输出边界切向电场
complex H(-1250:1250,4) !输出边界切向磁场
real Phi !观察角度
real StartAngle,EndAngle,DegreeStep !输出角87
C 建模
C Initiate
wavelength=0.12 !入射波长
wl=40 !一个波长采样40个点
MediaNumber=1 !0:真空, 1:代表目标媒质类型
media(1,1)=4. !epsilon_r
media(2,1)=1. !mu_r
media(3,1)=3.72E+7 !sigma
media(4,1)=0. !sigma_m
C 建立散射体
ob=0
do i=Imin,Imax
do j=Jmin,Jmax
if((i.ge.-wl).and.(i.le.wl-1).and. !判断网格是否在目标上
(j.ge.-wl).and.(j.le.wl-1))then
ob(i,j)=1
endif
enddo
enddo88
C FDTD 计算近场
C* 参数初始化
write(*,*)'Chose wave mode please(TM or TE):'
read(*,*) WaveMode
if (wavemode.eq.'TM'.or.wavemode.eq.'tm') then
TEMFlag=1
else if (wavemode.eq.'TE'.or.wavemode.eq.'te') then
TEMFlag=-1
end if
T1=Itmin*Itmin+Jtmin*Jtmin
T2=Itmax*Itmax+Jtmin*Jtmin
T3=Itmin*Itmin+Jtmax*Jtmax
T4=Itmax*Itmax+Jtmax*Jtmax
Isource=-sqrt(max(T1,T2,T3,T4))-3
write(*,*)'Incident angle(degree) is :'
read(*,*) IncidentAngle
Temp=(Itmax-Itmin)*(Itmax-Itmin)+(Jtmax-Jtmin)*(Jtmax-Jtmin)
Temp=6*sqrt(Temp) ! 时间步估算
write(*,'(11X,12hTimestep ( >,i5,8h ) is : ,$)')int(Temp)
read(*,*) TimeStop
IncidentStart = -710
IncidentEnd = 710
Pi = 3.14159265 ! 参数初始化
MU0 = Pi * .0000004
EPS0 = 8.85E-12
Z = sqrt(Mu0 / Eps0)
k=2*Pi/WaveLength89
Media(1,0)=1. ! 背景介质为自由空间
Media(2,0)=1.
Media(3,0)=0.
Media(4,0)=0.
do i=1,MediaNumber
Media(3,i)=-Media(3,i)*Z/k
Media(4,i)=-Media(4,i)/(Z*k)
end do
IncidentAngle = Pi * IncidentAngle / 180.
if (TEMflag.EQ.-1) then
EPS0=Pi*.0000004 ! TE 模式
MU0=8.85E-12
Z = sqrt(Mu0 / Eps0)
end if
FE=.5*TEMflag
FH=.5*TEMflag
SMU0=SQRT(MU0)
SEPS0=SQRT(EPS0)
OMIGA=Pi/WL
cp=COS(IncidentAngle)
sp=SIN(IncidentAngle)90
if (TEMFlag.eq.1) then ! 计算FDTD迭代公式系数 ,TM
do i=0,MediaNumber
do j=0,MediaNumber
FE1(i,j)=((Media(1,i)+Media(1,j))+ !CP,CQ,CA,CB的均值
+ Pi*.5*(Media(3,i)+Media(3,j))/WL)/
+ ((Media(1,i)+Media(1,j))-
+ Pi*.5*(Media(3,i)+Media(3,j))/WL)
FE2(i,j)=1./((Media(1,i)+Media(1,j))-
+ Pi*.5*(Media(3,i)+Media(3,j))/WL)
FH1(i,j)=((Media(2,i)+Media(2,j))+
+ Pi*.5*(Media(4,i)+Media(4,j))/WL)/
+ ((Media(2,i)+Media(2,j))-
+ Pi*.5*(Media(4,i)+Media(4,j))/WL)
FH2(i,j)=1./((Media(2,i)+Media(2,j))-
+ Pi*.5*(Media(4,i)+Media(4,j))/WL)
end do
end do
else
do i=0,MediaNumber
do j=0,MediaNumber
FH1(i,j)=((Media(1,i)+Media(1,j))+ !CP,CQ,CA,CB的均值,TE
+ Pi*.5*(Media(3,i)+Media(3,j))/WL)/
+ ((Media(1,i)+Media(1,j))-
+ Pi*.5*(Media(3,i)+Media(3,j))/WL)
FH2(i,j)=-1./((Media(1,i)+Media(1,j))-
+ Pi*.5*(Media(3,i)+Media(3,j))/WL)
FE1(i,j)=((Media(2,i)+Media(2,j))+
+ Pi*.5*(Media(4,i)+Media(4,j))/WL)/
+ ((Media(2,i)+Media(2,j))-
+ Pi*.5*(Media(4,i)+Media(4,j))/WL)
FE2(i,j)=-1./((Media(2,i)+Media(2,j))-
+ Pi*.5*(Media(4,i)+Media(4,j))/WL)
end do
end do
end if91
C Main Loop
OUTflag='I'
N=0
999 CONTINUE
do while (N.LT.TimeStop) ! 主循环
N=N+1
write(*,1) N
1 FORMAT(1H+,10X,'Time step',I5,' is in processing...... ')
******************* 产生入射波 *************************
do i=IncidentStart,IncidentEnd-1
Hin(i)=Hin(i)-FH*(Ein(i+1)-Ein(i)) !1D FDTD for Hin
end do
do i=IncidentStart+1,IncidentEnd-1
Ein(i)=Ein(i)+FE*(Hin(i-1)-Hin(i)) !1D FDTD for Ein
end do
c Set the unit source
Ein(Isource)=sin(OMIGA*N) !入射波源
if (N.le.WL) Ein(Isource)=Ein(Isource)*.5*(1.-cos(OMIGA*N))
!开关函数(升余弦函数)
Ein(IncidentStart) = EBin(4) !入射波吸收边界
EBin(4) = EBin(3)
EBin(3) = Ein(IncidentStart+1)
Ein(IncidentEnd) = EBin(2)
EBin(2) = EBin(1)
EBin(1) = Ein(IncidentEnd-1)92
*************** Ez 的FDTD迭代 ****************
do I=Imin+1,Imax-1
do J=Jmin+1,Jmax-1
T1=Hy(I,J)-Hy(I-1,J)-Hx(I,J)+Hx(I,J-1)
Ez(I,J)=FE1(ob(i,j),ob(i,j))*Ez(I,J)+
+ FE2(ob(i,j),ob(i,j))*T1
end do
end do
****** 连接边界引入入射波Ez分量 **************
do I=Itmin+1,Itmax-1
Ez(I,Jtmax)=Ez(I,Jtmax)-FE*HCB(I,3)
Ez(I,Jtmin)=Ez(I,Jtmin)+FE*HCB(I,1)
end do
do J=Jtmin+1,Jtmax-1
Ez(Itmax,J)=Ez(Itmax,J)+FE*HCB(J, 2)
Ez(Itmin,J)=Ez(Itmin,J)-FE*HCB(J, 4)
end do
****** 连接边界角点引入入射波 ****************
T1=HCB(Jtmax,2)-HCB(Itmax,3)
Ez(Itmax,Jtmax)=Ez(Itmax,Jtmax)+FE*T1
T1=-HCB(Jtmin,4)+HCB(Itmin,1)
Ez(Itmin,Jtmin)=Ez(Itmin,Jtmin)+FE*T1
T1=HCB(Jtmin,2)+HCB(Itmax,1)
Ez(Itmax,Jtmin)=Ez(Itmax,Jtmin)+FE*T1
T1=-HCB(Itmin,3)-HCB(Jtmax,4)
Ez(Itmin,Jtmax)=Ez(Itmin,Jtmax)+FE*T193
C 输出结果
if (OUTflag.EQ.'I') then !如果输出标志为'I',则输出场量为虚部
write(*,'(1H+,10X,40h Now outputing IMAGE part ...... )')
OPEN(2,FILE='FDTD.BI')
EI0=Ein(0)
else !如果输出标志为'R',则输出场量为实部
write(*,'(1H+,10X,40H Now outputing REAL part ...... )')
OPEN(2,FILE='FDTD.BR')
ER0=Ein(0)
end if
do I=Iomin,Iomax ! 输出输出边界上场量的实部或虚部到临时文件
T1=.5*(Ez(I,Jomin)+Ez(I,Jomin-1)) !下
write(2,*)T1,Hx(I,Jomin-1)/Z
T1=.5*(Ez(I,Jomax)+Ez(I,Jomax+1)) !上
write(2,*)T1,Hx(I,Jomax)/Z
end do
do J=Jomin,Jomax
T1=.5*(Ez(Iomax,J)+Ez(Iomax+1,J)) !右
write(2,*)T1,Hy(Iomax,J)/Z
T1=.5*(Ez(Iomin,J)+Ez(Iomin-1,J)) !左
write(2,*)T1,Hy(Iomin-1,J)/Z
end do
CLOSE(2)
if (OUTflag.EQ.'I') then
OUTflag='R'
TimeStop=TimeStop+WL/2 !延迟WL/2时间步以输出实部
GOTO 999
end if94
C 将输出结果组成复数
write(*,'(1H+,10X,40h Make .PIC files ...... )')
CT1=(1.,0.)/CMPLX(ER0,EI0) ! 入射波电场在原点处的值的倒数
CT2=CEXP(CMPLX(0,-.5*OMIGA))*CT1 !CT2比CT1倒数滞后半个时间步
!CT1,CT2将用来分别对电场和磁场归一
OPEN(1,FILE='FDTD.BI')
OPEN(2,FILE='FDTD.BR')
OPEN(3,FILE='rect.BND')
do iii=1,10000
read(1,*,end=1111)T1,T2
read(2,*)T3,T4
CT3=CMPLX(T3,T1)*CT1
write(3,*)CT3
CT4=CMPLX(T4,T2)*CT2
write(3,*)CT4
end do
1111 continue
CLOSE (1)
CLOSE (2)
CLOSE (3)95
C 外推远场区
C 初始化
open(1,file='rect.BND')
! write(*,*)'StartDegree,EndDegree,DegreeStep:'
! read(*,*)StartAngle,EndAngle,DegreeStep
StartAngle=0
EndAngle=180
DegreeStep=1
pi=3.1415927
rk=2*pi/float(wl) !contant of propogation K!
z=sqrt(pi*4.0e-7/8.854e-12) !wave inpedence of free space!
C 读入散射场数据Es,Hs, 并计算n x Es,n x Hs
do i=Iomin,Iomax
read(1,*)E(i,1) !Mz0=Esz,Jx0=Hsx*
read(1,*)H(i,1)
read(1,*)Ctemp
E(i,3)=-Ctemp !Mz0=-Esz
read(1,*)Ctemp1
H(i,3)=-Ctemp1 !Jx0=-Hsx
enddo
do i=Jomin,Jomax
read(1,*)E(i,2) !Mz0=Esz,Jx0=Hsx
read(1,*)H(i,2)
read(1,*)Ctemp
E(i,4)=-Ctemp !Mz0=-Esz
read(1,*)Ctemp1
H(i,4)=-Ctemp1 !Jx0=-Hsx
enddo
close(1)96
if(TEMflag.eq.-1) then !若为TE模,则作如下转换
do i=Iomin,Iomax
评论2