nearfield.f

FDTD近场分析

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *                                                                               * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
* *                                                                         * *
* *     2D Finite-Difference Time-Domain Algorithm Program                  * * 
* *                                                                         * *
* *   TYPE        : MAIN PROGRAM                                            * *
* *   SOURCE FILE : FDTD.F                                                  * *
* *   LANGUAGE    : Fortran77                                               * *
* *   CPU         : Intel                                                   * *
* *   CALLS       : NOTE                                                    * *
* *   ADDRESS     : XiDian University,Xi`an,710071                          * *
* *   Input File  : ??????(.SCT)                                            * *
* *   Output File : ?????? AZ.PIC   ?????? AX.PIC   ?????? AY.PIC           * *
* *               : ?????? PZ.PIC   ?????? PX.PIC   ?????? PY.PIC           * *
* *               : ??????.BDN                                              * *
* *   Note        : Output Data Normalized by z-Component of Incident Wave  * *
* *                                                                         * *
* *                                                                         * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
* * * * * * * * * * * * * * * 变量定义* * * * * * * * * * * * * * * * * * * * *
      real Ez(-160:160,-160:160)          ! FDTD采样点Ez(TM波，TE波
	                                    ! 为Hz,以下类似)
	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)                 !

	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 Phi                             !入射角
	integer TimeStop                     !总时间步

	real MaxX
      real MinX
	real MaxY
	real MinY
	real Step
	
	integer XNum      
      integer YNum
	
	integer Imax                         !FDTD区域总边界 
      integer Imin 
      integer Jmax
	integer Jmin
	
	integer Itmax                        !连接边界
	integer Itmin
	integer Jtmax
	integer Jtmin
	
	integer Iomax                       !输出边界
	integer Iomin
	integer Jomax
	integer Jomin
	
	integer Ismax                       !目标区边界     
	integer Ismin 
      integer Jsmax
	integer Jsmin

	character*30 FileName            !目标信息文件名
      character*1 Blank                !空格字符
      integer FileNameLength
	logical FileExist

      character*1 OUTflag   !输出标志：“R”输出场量的实部，“I”为虚部
      character*1 Respond   !响应回车键字符
      character*2 WaveMode
      character*30 rem

	real k,Z,WaveLength      !自由空间波数、波阻抗、波长

	real T1,T2,T3,T4,T5,T6   !临时变量
	complex CT1,CT2,CT3,CT4
	integer I,J,N,II,III,IFlag,NN

	real Media(1:4,0:10)      !介质参数数组
	integer MediaNumber       !目标所包含介质总数目
	
* * * * * * * * * * * * * * * 初始化 * * * * * * * * * * * * * * * * * *
	
10000 write(*,*) char(7)
      write(*,'(/////////////////////////)') 
	write(*,'(11X,48h2D-FDTD Program Ver. 6.0 for Intel Compatibility)')
      write(*,'(11X,"32h DarkStairs Ltd.")')
	write(*,'(11X,"46h Maximum FDTD zone is:-160:160,-160:160")')          
      write(*,'(11X,46h Maximum total field zone is:-150:150,-150:150)')     
	write(*,'(11X,"32h Maximum kinds of media is: 10")')                   
	write(*,*)
      write(*,'(11X,"34hThe scatterer data file name is:",$)')
	 
	read(*,'(a30)')FileName

	Blank=''
      FileNameLength=index(FileName,Blank)-1               
      FileName=FileName(1:FileNameLength)//'.SCT'

	inquire(file=FileName,exist=FileExist)

	if(FileExist)then
	  OPEN(1,file=FileName,access='direct',recl=4)
	else
	  write(*,'(1x,15x,15hFile not found.)')
        write(*,'(1x,15x,33hPress [Enter] to redo from start.)')
	  read(*,'(a10)') Respond
	  goto 10000
	end if

	read(1,rec=1) WaveLength
      read(1,rec=2) Temp
	WL=Temp                                 

      write(*,'(11X,"34hChose wave mode please(TM or TE):,$")')
	read(*,'(a2)') WaveMode
	if(wavemode.eq.'TM'.or.wavemode.eq.'tm')then
	  TEMFlag=1
	else if(wavemode.eq.'TE'.or.wavemode.eq.'te')then
	  TEMFlag=-1
	else
	  close(1)
        write(*,'(1x,15x,16hWrong wavemode !)')
        write(*,'(1x,15x,33hPress [Enter] to redo from start.)')
        read(*,'(a1)') Respond
        goto 10000
      end if

      read(1,rec=3) Temp
	Ismin=Temp
      read(1,rec=4) Temp
      Ismax=Temp
      read(1,rec=5) Temp
      Jsmin=Temp
      read(1,rec=6) Temp
      Jsmax=Temp

	print *,''
      write(*,'(11X,"15hKnown message:")')
      write(*,'(11X,11hWaveLength=,e12.4,4H m=, i4,5H Grid)')
     +     WaveLength,WL
      write(*,'(11X,38hMinimum connective boundary is (Grid):)')
      write(*,'(15X,5hXmin=,I4,5x,5hXmax=,I4)')Ismin-5,Ismax+5
      write(*,'(15x,5hYmin=,I4,5x,5hYmax=,I4)')Jsmin-5,Jsmax+5
      print *,''

      write(*,'(11X,"19hMake choice below:")')
      write(*,'(11X,25hConnective boundary is:,$)')
	read(*,*)Itmin,Itmax,Jtmin,Jtmax
	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(*,'(11X,"25hAbsorbing boundary is:",$)')
      read(*,*)Imin,Imax,Jmin,Jmax
      write(*,'(11X,25hThe output boundary is:,$)')
      read(*,*)Iomin,Iomax,Jomin,Jomax

      write(*,'(11X,27hIncident angle(degree) is:,$)')
      read(*,*) Phi

	Temp=(Itmax-Itmin)*(Itmax-Itmin)+(Jtmax-Jtmin)*(Jtmax-Jtmin)
	Temp=6*sqrt(Temp)                                                  
      write(*,'(11X,12hTimestep ( >,i5,8h ) is : , $ )')int(Temp)
      read(*,*)TimeStop

      write(*,*)char(7)
      write(*,'(11X,"23hOK ! Now initializing...")')

	read(1,rec=7)Temp
	MediaNumber=int(Temp)
	nn=7

	do j=1,MediaNumber
	   do i=1,4
	     nn=nn+1
	     read(1,rec=nn)Media(i,j)
	   end do
	end do

	do i=Imin,Imax                         
	   do j=Jmin,Jmax
	     ob(i,j)=0
	   end do
      end do
       
      do i=Ismin,Ismax                       
	  do j=Jsmin,Jsmax
	    nn=nn+1
	    read(1,rec=nn)temp
          ob(i,j)=temp
	  end do
      end do
	close(1)

	IncidentStart=-710               
      IncidentEnd=710

	Pi=3.14159265
	MU0=Pi*.0000004
	EPS0=8.85E-12
	Z=sqrt(Mu0/Eps0)
	k=2*Pi/WaveLength

c sigma=3.72E+07 sigmam=0 sigr=-sigma*Z/k sigmr=-sigmam/(Z*k)
      
	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

	Phi=Pi*Phi/180.
	if(TEMflag.EQ.-1)then
	  EPS0=Pi*.0000004
	  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(Phi)
	sp=SIN(Phi)

	if(TEMFlag.eq.1)then     !计算FDTD迭代公式系数
	do i=0,MediaNumber
	  do j=0,MediaNumber
	    FE1(i,j)=((Media(1,i)+Media(1,j))+                  
     +           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))+
     +           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 if 
	
***********************************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)
*******************************连接边界处入射场各分量****************************************
      do I=Itmin,Itmax


	  T1=float(I)*cp+float(Jtmin)*sp
	  II=int(T1)
	  if(T1.GE.0.) then
	    III=II+1
	    T1=float(III)-T1
	  else
	    III=II-1
	    T1=T1-float(III)
	  end if 
	  ECB(I,1)=T1*(Ein(II)-Ein(III))+Ein(III)    !底部连接边界处入射电场分量


	  T1=float(I)*cp+float(Jtmax)*sp
	  II=int(T1)
	  if(T1.GE.0.) then
	    III=II+1
	    T1=float(III)-T1
	  else
	    III=II-1
	    T1=T1-float(III)
	  end if 
	  ECB(I,3)=T1*(Ein(II)-Ein(III))+Ein(III)    !顶部连接边界处入射电场分量

      end do

	do J=Jtmin,Jtmax

	  T1=float(J)*sp+float(Itmin)*cp
	  II=int(T1)
	  if(T1.GE.0.) then
	    III=II+1
	    T1=float(III)-T1
	  else
	    III=II-1
	    T1=T1-float(III)
	  end if 
	  ECB(J,4)=T1*(Ein(II)-Ein(III))+Ein(III)    !左边连接边界处入射电场分量
      

	  T1=float(J)*sp+float(Itmax)*cp
	  II=int(T1)
	  if(T1.GE.0.) then
	    III=II+1
	    T1=float(III)-T1
	  else
	    III=II-1
	    T1=T1-float(III)
	  end if 
	  ECB(J,2)=T1*(Ein(II)-Ein(III))+Ein(III)    !右边连接边界处入射电场分量

      end do

      DO I=Itmin,Itmax
	 
	  T1=float(I)*cp+(float(Jtmin)-.5)*sp
        if(T1.GT..5) then
	    II=int(T1-.5)
		III=II+1
		T1=float(III)+.5-T1
	  else
	    II=int(T1+.5)
		III=II-1
		T1=T1-.5-float(III)
	  end if
	  HCB(I,1)=T1*(Hin(II)-Hin(III))+Hin(III)
	  HCB(I,1)=HCB(I,1)*sp            !底部连接边界入射波磁场分量(Hx)
	  
      
	  T1=float(I)*cp+(float(Jtmax)+.5)*sp
        if(T1.GT..5) then
	    II=int(T1-.5)
		III=II+1
		T1=float(III)+.5-T1
	  else
	    II=int(T1+.5)
		III=II-1
		T1=T1-.5-float(III)
	  end if
	  HCB(I,3)=T1*(Hin(II)-Hin(III))+Hin(III)
	  HCB(I,3)=HCB(I,3)*sp            !顶部连接边界入射波磁场分量(Hx)
	  

	end do


      DO J=Jtmin,Jtmax
	   T1=float(J)*sp+(float(Itmin)-.5)*cp
         if(T1.GT..5) then
	     II=int(T1-.5)
		 III=II+1
		 T1=float(III)+.5-T1
	   else
	     II=int(T1+.5)
		 III=II-1
		 T1=T1-.5-float(III)
	   end if
	   HCB(J,4)=T1*(Hin(II)-Hin(III))+Hin(III)
	   HCB(J,4)=-HCB(J,4)*cp            !左边连接边界入射波磁场分量(Hx)
	  
      
	   T1=float(J)*sp+(float(Itmax)+.5)*cp
         if(T1.GT..5) then
	     II=int(T1-.5)
		 III=II+1
		 T1=float(III)+.5-T1
	   else
	     II=int(T1+.5)
		 III=II-1
		 T1=T1-.5-float(III)
	   end if
	   HCB(J,2)=T1*(Hin(II)-Hin(III))+Hin(III)
	   HCB(J,2)=-HCB(J,2)*cp            !右边连接边界入射波磁场分量(Hx)
	end do

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