!*********************************************************************************************************    
!****************************************END OF PROGRAM****************************************    
!***********************************************************************************************************

!外推遠區場（或者RCS）

!***************************************************************
!**  Extrapolation of Radar Cross Section of Far Field       **
!**      Based on FDTD Calculations                          **
!**                                                          **
!**      TYPE               :main program                    **
!**      CALLS              :Bessel Functions                **
!**     LANGUAGE            :Fortran90                       **
!**     INTPUT FILE NAME    :????????.bnd                    **
!**     OUTPUT FILE NAME    :????????.dat                    **
!***********************************************************************
!***********************變量定義****************************************
     complex Jz      !輸出面上電流密度（TM極化，TE極化為磁流密度，以下類似）
     complex Mx      !磁流密度
     complex My      !磁流密度
	 complex Ctemp,ctemp1   !臨時復數變量
	 complex E(-1250:1250,4) !輸出邊界切向電場
	 complex H(-1250:1250,4) !輸出邊界切向磁場
	 real WaveLength         !波長
	 integer WL              !數值波長
	 integer TEMFlag         ! TM=1 TE=-1
	 real IncidentAngle      !FDTD計算中的入射角（度）
	 real Phi                !觀察角（度）
	 real StartAngle,EndAngle,DegreeStep !輸入角起始值、終止值、布進
	 character *30 FileName  !散射數據（由FDTD近場計算程序得到）
     character *1 Blank,Respond
	 logical FileExist
	 integer FileNameLength, Record
 !*************初始化**************

 !write(*,*) char(7)
 10002       write(*,'(///////////////////)')
	   write(*,*)''
	   write(*,*)'Extrapolate Base on FD-TD Boundary to Far Zone'
	   write(*,*)'to Caculate Radar Cross Section'
	   write(*,*)'(Two-Dimension&TM mode)'
	   write(*,*)''
       write(*,*)''
	   Blank=''
	   WRITE(*,'(11X,34hThe scatterer data file name is :,$)')!
	   READ(*,'(a30)')FileName
	   FileNameLength=index(FileName,Blank)-1
	   inquire(file=FileName(1:FileNameLength)//'.BND',exist=FileExist)
	   if(FileExist)then
	     OPEN(1,file=FileName(1:FileNameLength)//'.BND',access='direct',recl=8)
	   else
	      write(*,'(1x,5x,15hFile not found.)')
		  write(*,'(1x,15x,33hPress [Enter] to redo from start.)')
		  read(*,'(a1)')Respond
		  goto 10002
		  endif
	   read(1,rec=1)WaveLength
       read(1,rec=2)WL
	   read(1,rec=3)TEMFlag
	   read(1,rec=4)IncidentAngle
	   read(1,rec=5)Imin
	   read(1,rec=6)Imax
	   read(1,rec=7)Jmin
	   read(1,rec=8)Jmax
	   read(1,rec=9)Itmin
	   read(1,rec=10)Itmax
	   read(1,rec=11)Jtmin
	   read(1,rec=12)Jtmax
	   read(1,rec=13)Iomin
	   read(1,rec=14)Iomax
	   read(1,rec=15)Jomin
	   read(1,rec=16)Jomax
	   read(1,rec=17)TimeStep
       write(*,*)''
	   WRITE(*,'(11X,15hKnown message :)')
	   WRITE(*,'(13x,11hWaveLength=,e12.4,4H m = ,i4,5H Grid)')WaveLength,WL
	     if(TEMFlag.eq.1) then
		   write(*,'(13x,12hWaveMode:TM,5x,17hIncident angle =,f7.2,8h degrees)')IncidentAngle
	   else
           write(*,'(13x,12hWaveMode:TE,5x,17hIncident angle =,f7.2,8h degrees)')IncidentAngle
	   end if

	   write(*,'(13X,25hAbsorbing boundary is : ,4i6)')Imin,Imax,Jmin,Jmax
	   write(*,'(13X,25hConnective boundary is : ,4i6)')Itmin,Itmax,Jtmin,Jtmax
	   write(*,'(13X,25hThe output boundary is :,4i6)')Iomin,Iomax,Jomin,Jomax
	   write(*,'(13X,14hTimeStop is : ,i6)')TimeStep
	   write(*,*)''
	   write(*,'(11X,19hMake choice below :)')
	   write(*,'(13x,34hStartDegree,EndDegree,DegreeStep: ,$)')
	   read(*,*)StartAngle,EndAngle,DegreeStep
!begin
	   pi=3.1415927
	   rk=2*pi/float(WL)           !contant of propogation K!
	   Z=sqrt(pi*4.0e-7/8.854e-12) !wave inpedence of free space!
	   !********************讀入散射場數據Es,Hs***********************
	   !********************并計算n x Es, n x Hs**********************
	   Record=17
	   do i=Iomin,Iomax
	      Record=Record+1
		  read(1,rec=Record)E(i,1)     !Mz0=Esz,Jx0=Hsx
		  Record=Record+1
          read(1,rec=Record)H(i,1) 
          Record=Record+1
		  read(1,rec=Record)Ctemp      
		  E(i,3)=-Ctemp               !Mz0=-Esz
		  Record=Record+1
		  read(1,rec=record)Ctemp1
		  H(i,3)=-Ctemp1               !Jx0=-Hsx
		  end do
		  
		  
		  do i=Jomin,Jomax
		      Record=Record+1
		      read(1,rec=Record)E(i,2)    !Mz0=Esz,Jx0=Hsx
			  Record=Record+1
			  read(1,rec=Record)H(i,2)
			  Record=Record+1
			  read(1,rec=Record)Ctemp
			  E(i,4)=-Ctemp             !Mz0=-Esz
			  Record=Record+1
			  read(1,rec=Record)Ctemp1
			  H(i,4)=-Ctemp1            !Jy0=-Hsy
			 end do
			 close(1)
	   
	   if(TEMFlag.eq.-1)then    !若為TE模，則做以下轉換
	   do i=iomin,iomax
	     E(i,1)=E(i,1)            !Ez-->Hz
	     E(i,3)=E(i,3)
	     H(i,1)=-H(i,1)           !Hx-->-Ex
	     H(i,3)=-H(i,3)
	   end do
	   do i=jomin,jomax
	     E(i,2)=E(i,2)            !Ez-->Hz
	     E(i,4)=E(i,4)
	     H(i,2)=-H(i,2)           !Hy-->-Ey
	     H(i,4)=-H(i,4)
	   end do
	   Z=1./Z
	   end if                    !reciprocal
	   !******************沿著四條輸出邊界積分Jz,Mx,My**********************
	   OPEN(2,file=FileName(1:FileNameLength)//'.dat') !存儲RCS的文件
	   

	   do dg=StartAngle,EndAngle,DegreeStep             !散射角循環
	     Jz=(0.,0.)
	     Mx=(0.,0.)
	     My=(0.,0.)
	     write(*,1)dg
1        format(1H+,10x,'Now caculate degree',f6.1)
	     Phi=dg*pi/180.


	     do I=Iomin,Iomax
	       Ctemp=cmplx(0.,(float(I)*cos(Phi)+(jomin-0.5)*sin(Phi))*rk)
	       Ctemp=cexp(Ctemp)                                 !exp(j*k*r);x=i,j=jomin-.5
	       Mx=Mx+E(i,1)*Ctemp                                 !下邊界
	       Jz=Jz+H(i,1)*Ctemp
           Ctemp=cmplx(0.,(float(I)*cos(Phi)+(jomax+.5)*sin(Phi))*rk)  
	       Ctemp=cexp(Ctemp)                                !exp(j*k*r);x=i,j=jomin-.5
	       Mx=Mx+E(i,3)*Ctemp                                !上邊界
		   Jz=Jz+H(i,3)*Ctemp
		end do

		do I=Jomin,Jomax
        Ctemp=cmplx(0.,(float(I)*sin(Phi)+(iomax+.5)*cos(Phi))*rk)
		Ctemp=cexp(Ctemp)                                   !exp(j*k*r);x=iomax+.5,y=i
        My=My+E(i,2)*Ctemp                                  !右邊界
		Jz=Jz+H(i,2)*Ctemp
        Ctemp=cmplx(0.,(float(I)*sin(Phi)+(iomin-.5)*cos(Phi))*rk) 
		Ctemp=cexp(Ctemp)                                   !exp(j*k*r);x=iomin-.5,y=i
		My=My+E(i,4)*Ctemp                                  !左邊界
         Jz=Jz+H(i,4)*Ctemp
		 end do 
!*************************計算雷達散射截面(RCS)并輸出************************
            
			
			 RCS=abs(Mx*sin(Phi)-My*cos(Phi)+Z*Jz)
             RCS=.25*rk*rcs*rcs*Wavelength/float(WL)/Wavelength
			 write(1,*)dg,rcs
			 RCS=10.*log10(RCS)
			 write(2,*)dg,RCS                      
			 end do                                         !結束散射角循環
            close(2)
            write(*,*) CHAR(7)                
			end
!****************************************************************************
!****************************END OF PROGRAM**********************************
!****************************************************************************