total2.f90

fortran source code for fdtd parall computing

!************************************************************************************
!***************** 2D Finite-Difference Time-Domain Algorithm Program****************
!************************************************************************************


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采样点
real Einc(-800:800)
integer IncidentStart,IncidentEnd,Isource   !一维入射波数组的端点和原始位置
real EBin(4)                                 !入射波吸收边界条件临时变量

real ECB(-150:150,4)                         !连接边界处E入射波分量
real HCB(-150:150,4)                         !连接边界处H入射波分量

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                            ! 目标所包含介质总数目

!*******************************初始化*********************************************
 ! write(*,*)char(7)
10001 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,32hThe scatterer data file name is: ,$)')
	 ! write(*,'(11X,32hThe scatterer data file name is:,$)')
      !read(*,'(a30)') FileName
      FileName='rect'
      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 10001
      end if

      read(1,rec=1) WaveLength
      read(1,rec=2) Temp
      WL=Temp                                                ! 数值波长
      WaveMode='TM'
      write(*,'(11X,33hChose 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 10001
      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,10hTimestep(>,i5,9h )is : ,$)')int(Temp)
 read(*,*)TimeStop

 write(*,*)char(7)
 write(*,'(11X,27hOK    ! 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

! 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             !TE模式
    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

 !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           ! 左边连接别界入射波磁场分量（Hy）

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                      !右边连接边界入射波磁场分量(Hy)
   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)
     Ez(I,J)=FE1(ob(i,j),ob(i,j))*Ez(I,J)+FE2(ob(i,j),ob(i,j))*T1