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