?? nearfield.f
字號:
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*T1
*************************吸收邊界條件Ez****************************************
do I=Imin+1,Imax-1
T1=Hy(I,Jmin)-Hy(I-1,Jmin)
T1=T1+Hy(I,Jmin+1)-Hy(I-1,Jmin+1)
T1=T1*TEMFlag
T2=Ez(I,Jmin+1)-Ez(I,Jmin)
Ez(I,Jmin)=EAB(I,1)-.33333333*T2
Ez(I,Jmin)=Ez(I,Jmin)+.16666667*T1
EAB(I,1)=Ez(I,Jmin+1)
T1=Hy(I,Jmax)-Hy(I-1,Jmax)
T1=T1+Hy(I,Jmax-1)-Hy(I-1,Jmax-1)
T1=T1*TEMFlag
T2=Ez(I,Jmax-1)-Ez(I,Jmax)
Ez(I,Jmax)=EAB(I,3)-.33333333*T2
Ez(I,Jmax)=Ez(I,Jmax)+.16666667*T1
EAB(I,3)=Ez(I,Jmax-1)
end do
do J=Jmin+1,Jmax-1
T1=Hx(Imin,J)-Hx(Imin,J-1)
T1=T1+Hx(Imin+1,J)-Hx(Imin+1,J-1)
T1=T1*TEMFlag
T2=Ez(Imin+1,J)-Ez(Imin,J)
Ez(Imin,J)=EAB(J,4)-.33333333*T2
Ez(Imin,J)=Ez(Imin,J)-.16666667*T1
EAB(J,4)=Ez(Imin+1,J)
T1=Hx(Imax,J)-Hx(Imax,J-1)
T1=T1+Hx(Imax-1,J)-Hx(Imax-1,J-1)
T1=T1*TEMFlag
T2=Ez(Imax-1,J)-Ez(Imax,J)
Ez(Imax,J)=EAB(J,2)-.33333333*T2
Ez(Imax,J)=Ez(Imax,J)-.16666667*T1
EAB(J,2)=Ez(Imax-1,J)
end do
****************************吸收邊界角點Ez*************************************
Iflag=1-Iflag ! 交替存儲標志
Ez(Imin,Jmin)=.2928932*EAC(1,1,Iflag)+.7071068*EAC(1,2,Iflag)
EAC(1,1,Iflag)=Ez(Imin,Jmin)
EAC(1,2,Iflag)=Ez(Imin+1,Jmin+1)
Ez(Imax,Jmin)=.2928932*EAC(2,1,Iflag)+.7071068*EAC(2,2,Iflag)
EAC(2,1,Iflag)=Ez(Imax,Jmin)
EAC(2,2,Iflag)=Ez(Imax-1,Jmin+1)
Ez(Imin,Jmax)=.2928932*EAC(3,1,Iflag)+.7071068*EAC(3,2,Iflag)
EAC(3,1,Iflag)=Ez(Imin,Jmax)
EAC(3,2,Iflag)=Ez(Imin+1,Jmax-1)
Ez(Imax,Jmax)=.2928932*EAC(4,1,Iflag)+.7071068*EAC(4,2,Iflag)
EAC(4,1,Iflag)=Ez(Imax,Jmax)
EAC(4,2,Iflag)=Ez(Imax-1,Jmax-1)
********************************Hx,Hy的FDTD迭代******************************
do I=Imin,Imax
do J=Jmin,Jmax-1
T1=Ez(I,J+1)-Ez(I,J)
Hx(I,J)=FH1(ob(i,j),ob(i,j+1))*Hx(I,J)-
+ FH2(ob(i,j),ob(i,j+1))*T1
end do
end do
do I=Imin,Imax-1
do J=Jmin,Jmax
T1=Ez(I+1,J)-Ez(I,J)
Hy(I,J)=FH1(ob(i,j),ob(i+1,j))*Hy(I,J)+
+ FH2(ob(i,j),ob(i+1,j))*T1
end do
end do
************************ 連接邊界引入入射波Hx,Hy分量*********************
do I=Itmin,Itmax
Hx(I,Jtmin-1)=Hx(I,Jtmin-1)+FH*ECB(I,1)
Hx(I,Jtmax)=Hx(I,Jtmax)-FH*ECB(I,3)
end do
do J=Jtmin,Jtmax
Hy(Itmin-1,J)=Hy(Itmin-1,J)-FH*ECB(J,4)
Hy(Itmax,J)=Hy(Itmax,J)+FH*ECB(J,2)
end do
end do
************************輸出結果*****************************
if(OUTflag.EQ.'I') then ! 如果輸出標志為'I',則輸出場量為虛部
write(*,'(1H+,10X,"40h Now outputing IMAGE part......")')
OPEN(1,FILE='FDTD.EI',ACCESS='DIRECT',RECL=4)
OPEN(4,FILE='FDTD.HxI',ACCESS='DIRECT',RECL=4)
OPEN(5,FILE='FDTD.HyI',ACCESS='DIRECT',RECL=4)
OPEN(2,FILE='FDTD.BI')
EI0=Ein(0)
else
write(*,'(1H+,10X,"40H Now outputing REAL part......")')
OPEN(1,FILE='FDTD.ER',ACCESS='DIRECT',RECL=4)
OPEN(4,FILE='FDTD.HxR',ACCESS='DIRECT',RECL=4)
OPEN(5,FILE='FDTD.HyR',ACCESS='DIRECT',RECL=4)
OPEN(2,FILE='FDTD.BR')
ER0=Ein(0)
end if
NN=0
do I=Imin,Imax
do J=Jmin,Jmax
NN=NN+1
write(1,REC=NN) Ez(I,J)
write(4,REC=NN) Hx(I,J)/Z
write(5,REC=NN) Hy(I,J)/Z
end do
end do
CLOSE(1)
CLOSE(4)
CLOSE(5)
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
GOTO 999
end if
****************************將輸出結果組成復數******************************
write(*,'(1H+,10X,"40h Make.PIC file...... ")')
CT1=(1.,0.)/CMPLX(ER0,EI0)
CT2=CEXP(CMPLX(0,-.5*OMIGA))*CT1 !CT2比CT1倒數滯后半個時間步長
!CT1!CT2將用來分別對電場和磁場歸一
OPEN(1,FILE='FDTD.EI',ACCESS='DIRECT',RECL=4)
OPEN(2,FILE='FDTD.ER',ACCESS='DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'AZ.PIC'
OPEN(4,FILE=FileName,ACCESS='DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'PZ.PIC'
OPEN(5,FILE=FileName,ACCESS='DIRECT',RECL=4)
NN=0
do I=Imin,Imax
do J=Jmin,Jmax
NN=NN+1
read(1,REC=NN) T1
read(2,REC=NN) T2
CT3=CMPLX(T2,T1)*CT1
write(4,REC=NN) ABS(CT3)
T1=ATAN2(IMAG(CT3),real(CT3))
write(5,REC=NN)T1
end do
end do
CLOSE(1)
CLOSE(2)
CLOSE(4)
CLOSE(5)
OPEN(6,FILE='FDTD.HxI',ACCESS='DIRECT',RECL=4)
OPEN(7,FILE='FDTD.HxR',ACCESS='DIRECT',RECL=4)
OPEN(8,FILE='FDTD.HyI',ACCESS='DIRECT',RECL=4)
OPEN(9,FILE='FDTD.HyR',ACCESS='DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'AX.dat'
OPEN(10,FILE=FileName,ACCESS='DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'PX.PIC'
OPEN(11,FILE=FileName,ACCESS='DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'AY.PIC'
OPEN(12,FILE=FileName,ACCESS='DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'PY.PIC'
OPEN(13,FILE=FileName,ACCESS='DIRECT',RECL=4)
NN=0
do I=Imin,Imax
do J=Jmin,Jmax
NN=NN+1
read(6,REC=NN)T1
read(7,REC=NN)T2
CT3=CMPLX(T2,T1)*CT1
write(10,REC=NN) ABS(CT3)
T1=ATAN2(IMAG(CT3),real(CT3))
write(11,REC=NN) T1
read(8,REC=NN) T1
read(9,REC=NN)T2
CT3=CMPLX(T2,T1)*CT1
write(12,REC=NN) ABS(CT3)
T1=ATAN2(IMAG(CT3),real(CT3))
write(13,REC=NN) T1
end do
end do
CLOSE(6)
CLOSE(7)
CLOSE(8)
CLOSE(9)
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(13)
OPEN(1,FILE='FDTD.BI')
OPEN(2,FILE='FDTD.BR')
FileName=FileName(1:FileNameLength)//'.BND'
OPEN(3,FILE=FileName,access='direct',recl=8)
write(3,rec=1)WaveLength
!輸出輸出邊界的場值(幅值和相位,用CT1或歸CT2一)
write(3,rec=2)WL
write(3,rec=3)TEMFlag
write(3,rec=4)Phi*180/Pi
write(3,rec=5)Imin
write(3,rec=6)Imax
write(3,rec=7)Jmin
write(3,rec=8)Jmax
write(3,rec=9)Itmin
write(3,rec=10)Itmax
write(3,rec=11)Jtmin
write(3,rec=12)Jtmax
write(3,rec=13)Iomin
write(3,rec=14)Iomax
write(3,rec=15)Jomin
write(3,rec=16)Jomax
write(3,rec=17)TimeStop-WL/2
NN=17
do iii=1,10000
read(1,*,end=1111)T1,T2
read(2,*)T3,T4
NN=NN+1
CT3=CMPLX(T3,T1)*CT1
write(3,rec=NN)CT3
NN=NN+1
CT4=CMPLX(T4,T2)*CT2
write(3,rec=NN)CT4
end do
1111 continue
CLOSE(1)
CLOSE(2)
CLOSE(3)
!輸出.REM(與.PIC對應)文件以便用Photo.exe顯示計算區域場值
Step=WaveLength/WL
MaxX=Imax*Step
MinX=Imin*Step
MaxY=Jmax*Step
MinY=Jmin*Step
XNum=Imax-Imin+1
YNum=Jmax-Jmin+1
FileName=FileName(1:FileNameLength)//'AZ.REM'
open(1,file=FileName)
write(1,*)'"Title"""""""'
write(1,*)'" "'
write(rem,101)MaxX
101 format(1x,'"MaxX:',f9.4,'m"')
write(1,*)rem
write(rem,102)MinX
102 format(1x,'"MinX:',f9.4,'m"')
write(1,*)rem
write(rem,103)MaxY
103 format(1x,'"MaxY:',f9.4,'m"')
write(1,*)rem
write(rem,104)MinY
104 format(1x,'"MinY:',f9.4,'m"')
write(1,*)rem
write(rem,105)Step
105 format(1x,'"Step:',f9.4,'m"')
write(1,*)rem
write(rem,106)XNum
106 format(1x,'"XNum:',I4,'(grid)"')
write(1,*)rem
write(rem,107)YNum
107 format(1x,'"YNum:',I4,'(grid)"')
write(1,*)rem
write(1,*)'" "'
write(1,*)'"WaveMode:'//WaveMode,' "'
write(rem,108)WaveLength
108 format(1x,'"WL=',f9.4,'m"')
write(1,*)rem
write(rem,109)WL
109 format(1x,'"=',I4,'grid"')
write(1,*)rem
write(1,*)'"Incident Angle"'
write(rem,110)Phi*180/Pi
110 format(1x,'" ',f7.2,'degree"')
write(1,*)rem
write(rem,111)TimeStop-WL/2
111 format(1x,'"TimeStep=',I5,' "')
write(1,*)rem
write(1,*)'" "'
write(1,*)'" "'
write(1,*)MaxX,MinX,MaxY,MinY
write(1,*)XNum,YNum,Step
close(1)
FileName=FileName(1:FileNameLength)//'PZ.REM'
open(1,file=FileName)
write(1,*)'"Title"""""""'
write(1,*)'" "'
write(rem,101)MaxX
write(1,*)rem
write(rem,102)MinX
write(1,*)rem
write(rem,103)MaxY
write(1,*)rem
write(rem,104)MinY
write(1,*)rem
write(rem,105)Step
write(1,*)rem
write(rem,106)XNum
write(1,*)rem
write(rem,107)YNum
write(1,*)rem
write(1,*)'" "'
write(1,*)'"WaveMode:'//WaveMode,'"'
write(rem,108)WaveLength
write(1,*)rem
write(rem,109)WL
write(1,*)rem
write(1,*)'"Incident Angle"'
write(rem,110)Phi*180/Pi
write(1,*)rem
write(rem,111)TimeStop-WL/2
write(1,*)rem
write(1,*)'" "'
write(1,*)'" "'
write(1,*)MaxX,MinX,MaxY,MinY
write(1,*)XNum,YNum,Step
close(1)
open(1,file='FDTD.EI') ! 刪除臨時文件
close(1,status='delete')
open(2,file='FDTD.ER')
close(2,status='delete')
open(3,file='FDTD.HXI')
close(3,status='delete')
open(4,file='FDTD.HXR')
close(4,status='delete')
open(5,file='FDTD.HYI')
close(5,status='delete')
open(6,file='FDTD.HYR')
close(6,status='delete')
open(7,file='FDTD.BI')
close(7,status='delete')
open(8,file='FDTD.BR')
close(8,status='delete')
end
**************************************************************************
**** END OF PROGRAM ******
**************************************************************************
****************************
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -