nearfield.f

FDTD近场分析

	    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                                     ******
**************************************************************************

****************************