📄 erwei.f90
字号:
program main
implicit none
!************************************************************************************************************************************************************************
!* 参变量设定 *
!************************************************************************************************************************************************************************
integer n,k,zongdian !后两变量为一维波源点，及水平位置中点
integer nlayer,net(0:20),nnn,nnn1 !路面结构层层数（包括空气层），各层网格数，总网格数,路面结构截取水平网格数
real Thickness(0:50),weigh,d,c !路面结构层各层的实际厚度及截取计算的路面水平长度，d为常变量=power(2,0.5),c为光速
real ez(0:500,0:500),hx(0:500,0:500),hy(0:500,0:500) !一维及二维电场与磁场变量
real er1(0:50),sig1(0:50),vp(0:500,0:500) !路面结构层各层的介电常数及电导率，以及一维二维在各网格点上的传播速度
real er_mod1(0:500),sig_mod1(0:500),er_mod(0:500,0:500),sig_mod(0:500,0:500) !一维及二维介电常数及电导率变量 !总场－散射场边界条件的范围划定值
real pulse,pi
real dt,dx,F,e0,u0,eta !一维及二维时间及空间步长及基本介电常数及电导率
real er,sig,ur
real e1(0:500,0:500),e2(0:500,0:500),e3(0:500,0:500),ez1(0:500,0:500) !一维及二维吸收边界条件差分公式参数
real d1(0:500,0:500),d2(0:500,0:500),c1(0:500,0:500),c2(0:500,0:500) !麦克斯韦方程二维维差分方程的系数
integer i,j,l,nsteps,mmm !nsteps为程序运行的总步长数
!***********************************************************************************************************************************************************************
! 读取各类已知参数及打开各类需要文件 *
!***********************************************************************************************************************************************************************
i=0
j=0
l=0
nsteps=0
!--------
open(3,file='input-parameter.txt',status="old") !该文件为读取相关结构层的厚度层数，各层及目标体的基本介电常数，电导率
read(3,*) weigh !读取总的计算区域的长度(cm)
read(3,*) nlayer !读取包括空气层在内的总层数
read(3,*) (Thickness(i),i=1,nlayer) !读取总共的层数（包括空气层）
read(3,*) (er1(i),i=1,nlayer) !读取各层的介电常数（包括空气层）
read(3,*) (sig1(j),j=1,nlayer) !读取各层的电导率（包括空气层）
read(3,*) dt,dx,nsteps,F,e0,u0,eta
open(11,file="zongchang.txt") !该文件为记录总场文件(不带顺序数n)
open(20,file="jiedianfenbu.txt")
!***********************************************************************************************************************************************************************
! 划分模型网格 *
!***********************************************************************************************************************************************************************
do i=1,20
net(i)=0
enddo
do i=1,200
er_mod1(i)=1.0
sig_mod1(i)=0.0
do j=1,200
er_mod(j,i)=1.0
sig_mod(j,i)=0.0
enddo
enddo
nnn1=0
zongdian=0.0
nnn=0
mmm=0
er1(0)=1.0
sig1(0)=0.0
ur=1.00
er=1.00
sig=0.00
do k=1,200
do i=1,200 !设置一维及二维的各变量初值
ez(i,k)=0.0
hx(i,k)=0.0
hy(i,k)=0.0
enddo
enddo
do i=1,200
do j=1,200
e1(i,j)=0.0
e2(i,j)=0.0
e3(i,j)=0.0
ez1(i,j)=0.0
d1(i,j)=0.0
d2(i,j)=0.0
c1(i,j)=0.0
c2(i,j)=0.0
enddo
enddo
!**********************************************************以上为部分变量付初值
!**********************************************************
nnn1=int(weigh/dx/100) !该部分为计算水平网格总数
zongdian=(nnn1/2)+1
nnn1=nnn1+1
do i=1,nlayer
net(i)=int(Thickness(i)/dx/100) !该部分为计算纵向网格总数
nnn=nnn+net(i)
write(30,*) net(i),"每层总共网格数截至处",nnn+1
enddo
write(30,*) nnn,"总网格数"
net(0)=0
do i=1,nlayer
nnn=0
mmm=1 !给各层一维的介电常数及电导率赋值
do j=1,i
nnn=nnn+net(j)
enddo
do j=1,i
mmm=mmm+net(j-1)
enddo
do l=mmm,nnn+1
er_mod1(l)=er1(i)
sig_mod1(l)=sig1(i)
enddo
enddo
nnn=0
do i=1,nlayer
nnn=nnn+net(i)
enddo
nnn=nnn+1
!--------------------------------------------------
write(20,*) "最终纵向每单元网格层的介电常数及电导率分布"
write(20,*) "***************************************************************************"
write(20,*) " 序列号 介电常数 电导率"
do i=1,nnn
write(20,*) i,er_mod1(i),sig_mod1(i) !输出总的介电常数及电导率
enddo
do i=1,nnn1
do j=1,nnn
er_mod(i,j)=er_mod1(j)
sig_mod(i,j)= sig_mod1(j)
enddo
enddo
c = 2.996e8 ! m/s: velocity of propagation in air
do j=1,nnn1
do i=1,nnn
vp(j,i)=c/sqrt(er_mod(j,i))
enddo
enddo
!下段为一维及二维麦克斯韦方程差分公式系数
!****************************************************************************************************
!------二维
do j=1,nnn1
do i=1,nnn
d1(j,i) = dt/(u0*ur*dx)
d2(j,i) = 1
c1(j,i) = (1.0-sig_mod(j,i)*dt/(2.0*er_mod(j,i)*e0))/(1.0+sig_mod(j,i)*dt/(2.0*er_mod(j,i)*e0)) !二维差分计算公式系数
c2(j,i) = dt/(dx*er_mod(j,i)*e0*(1.0+sig_mod(j,i)*dt/(er_mod(j,i)*e0*2.0)))
enddo
enddo
!***********************************************************************************************************************************************************************
! 开始计算引入波源，电场，磁场，及总场散射场边界，吸收边界条件 * *
!***********************************************************************************************************************************************************************
pi=3.1415926
write(20,*)
do n=1,nsteps
!write(*,*) zongdian
pulse=(dt*n)**2*exp(-0.93*2*pi*dt*n*8.0e+8)*sin(2*pi*dt*n*8.0e+8) !加入一维入射源
write(20,*) n,pulse
!do i=1,50
ez(54,2)=pulse+ez(54,2)
!enddo
!********************************************************************!Mur一维边界条件计算
do j=1,nnn
do i=1,nnn1
ez(i,j)=c1(i,j)*ez(i,j)+c2(i,j)*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1)) !计算二维电场
enddo
enddo
!**********************************************************************!以下为Mur二阶吸收边界条件
!mur-2 abc
d=1.414214
do j=1,nnn1
do i=1,nnn
e1(j,i)=(vp(j,i)*dt-dx)/(vp(j,i)*dt+dx)
e2(j,i)=(vp(j,i)*vp(j,i)*dt*u0)/(2.0*(vp(j,i)*dt+dx)) !二维四边的各边边界条件差分计算公式系数计算
e3(j,i)=(vp(j,i)*dt-d*dx)/(vp(j,i)*dt+d*dx)
enddo
enddo
!-----------------------------------------------------左边边界计算
!left左边
do j=2,nnn-1
ez(1,j)=ez1(2,j)+e1(1,j)*(ez(2,j)-ez1(1,j))-e2(1,j)*(hx(1,j)-hx(1,j-1)+hx(2,j)-hx(2,j-1))
enddo
!-----------------------------------------------------右边边界计算
!right右边
do j=2,nnn-1
ez(nnn1,j)=ez1(nnn1-1,j)+e1(nnn1,j)*(ez(nnn1-1,j)-ez1(nnn1,j))-e2(nnn1,j)*(hx(nnn1,j)-hx(nnn1,j-1)+hx(nnn1-1,j)-hx(nnn1-1,j-1))
enddo
!-----------------------------------------------------底边边界计算
!bottom底边
do i=2,nnn1-1
ez(i,1)=ez1(i,2)+e1(i,1)*(ez(i,2)-ez1(i,1))+e2(i,1)*(hy(i,2)-hy(i-1,2)+hy(i,1)-hy(i-1,1))
enddo
!------------------------------------------------------顶边边界计算
!top顶边
do i=2,nnn1-1
ez(i,nnn)=ez1(i,nnn-1)+e1(i,nnn)*(ez(i,nnn-1)-ez1(i,nnn))+e2(i,nnn)*(hy(i,nnn)-hy(i-1,nnn)+hy(i,nnn-1)-hy(i-1,nnn-1))
enddo
!******************************************************四个角点
!corner四个角点
ez(1,1)=ez1(2,2)+e3(1,1)*(ez(2,2)-ez1(1,1))
ez(1,nnn)=ez1(2,nnn-1)+e3(1,nnn)*(ez(2,nnn-1)-ez1(1,nnn))
ez(nnn1,1)=ez1(nnn1-1,2)+e3(nnn1,1)*(ez(nnn1-1,2)-ez1(nnn1,1))
ez(nnn1,nnn)=ez1(nnn1-1,nnn-1)+e3(nnn1,nnn)*(ez(nnn1-1,nnn-1)-ez1(nnn1,nnn))
do j=1,nnn
do i=1,nnn1
ez1(i,j)=ez(i,j)
enddo
enddo
!******************************************************************************************************************************************
do j=1,nnn-1
do i=1,nnn1
hx(i,j)=hx(i,j)-(dt/(u0*dx))*(ez(i,j+1)-ez(i,j)) !计算二维磁场hx
enddo
enddo
do j=1,nnn
do i=1,nnn1-1
hy(i,j)=hy(i,j)+(dt/(u0*dx))*(ez(i+1,j)-ez(i,j)) !计算二维磁场hy
enddo
enddo
write(11,*) ez(60,2) !输出总场结果
enddo
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -