!****************************************************************************
!
!  PROGRAM: Transfer_Matrix_Method
!
!****************************************************************************

	program TMM
	implicit none
	real(8),parameter::pi=3.141592653589793D0,error=0.01D0
	integer,parameter::m=15                  !slice number
	real(8)::a,f,die1,die2,theta_k,phi_k,theta_E,phi_E    !input parameter  
	real(8)::a1,a2,b1,b2,c1,c2,h,sh_i_1,y01,y02,theta_i,ko,kox,koy,koz,Eox,Eoy,Eoz,RT_result(1500,3),&
    RT_fun,temp,sum,temp_result(1500,2)                    !interior parameter
	integer::k0,n,nn,i,m_unit,RT                  
	real(8),dimension(:,:),allocatable::t1,p
	complex(8),dimension(:,:),allocatable::sa,ta,beta,S0,T0,s,layer1,layer2,slab
	complex(8),dimension(:),allocatable::EF0,EFn,Ez
	intrinsic sqrt,dcos,dsin,dcmplx,dabs
	external overlap,construc,STa,ST0,s_matrix,s_layer1,s_layer2,RT_fun,Ez_array
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    write(*,*)'The program can be used to calculate the band structure,&
             &reflectance and tansmission spectra of photonic crystal,&
             &by the plane wave based tansfer matrix method!'
    write(*,'(a)',advance='no')'Input the lattice constant a (mum)='
	read(*,*)a
	write(*,'(a)',advance='no')'Input the fraction filling f (r/a)='
	read(*,*)f
	write(*,'(a)',advance='no')'Input the dielectrc constant of periodic lattice='
	read(*,*)die1
    write(*,'(a)',advance='no')'Input the dielectrc constant of slab matrix='
	read(*,*)die2
	write(*,*)'The following is the incident angle of wave vector:'	
	write(*,'(a)',advance='no')'Polar angle theta(degree)='
	read(*,*)theta_k
	write(*,'(a)',advance='no')'Azimuthal angle phi(degree)='
	read(*,*)phi_k
    write(*,*)'The following is the polarization direcation of electric field.& 
    Amplitudes E0 is normalized. The direction should be vertical to wave vector.'
    write(*,'(a)',advance='no')'Polar angle theta(degree)='
	read(*,*)theta_E
   write(*,'(a)',advance='no')'Azimuthal angle phi(degree)='
	read(*,*)phi_E
    write(*,'(a)',advance='no')'the number of unit cell='
	read(*,*)m_unit
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	!Initialization 
	a=a*1.0D-4
	theta_E=theta_E*pi/180.0D0
	phi_E=phi_E*pi/180.0D0
    theta_k=theta_k*pi/180.0D0
	phi_k=phi_k*pi/180.0D0
    Eox=dsin(theta_E)*dcos(phi_E)
	Eoy=dsin(theta_E)*dsin(phi_E)		               
	Eoz=dcos(theta_E)        
    kox=dsin(theta_k)*dcos(phi_k)
    koy=dsin(theta_k)*dsin(phi_k)
    koz=dcos(theta_k)
	if((kox*Eox+koy*Eoy+koz*Eoz)>1.0D-6)then
	write(*,*)'The electric field isn''t vertical to wave vector, and input argument is wrong!'
	stop
	end if  
	if((f<0.0D0).OR.(f>0.5D0))then
	write(*,*)'the fraction filling ''f'' must be at the range of 0.0 to 0.5!'
	stop
	end if   
	
!**************************************************************************************  
!TM direction!!
!**************************************************************************************  
    n=10
	nn=n/2		
	temp_result=0.0D0 
	temp=1.0D0
	a1=a
	b1=2.0D0*pi/a1
TM:    if(f<=0.433012702D0)then  !*******************************************************
TMsw:do while(temp>error)
	  allocate(p(n,n),sa(n,n),ta(n,n),S0(n,n),T0(n,n),t1(n,n),beta(n,n))
	  allocate(s(2*n,2*n),layer1(2*n,2*n),layer2(2*n,2*n),slab(2*n,2*n))
	  allocate(EF0(n),EFn(n),Ez(nn))
      EF0=(0.0D0,0.0D0)
	  EF0(nn)=dcmplx(Eox)
	  EF0(nn+1)=dcmplx(Eoy)
TMsn:do k0=1,1500             !**************************** the range of wavenumber!!!!
        ko=20.0D0*pi*k0	         
     	RT_result(k0,1)=ko              
		kox=ko*dsin(theta_k)*dcos(phi_k) 
       	koy=ko*dsin(theta_k)*dsin(phi_k) 
	    koz=ko*dcos(theta_k) 	
	  layer1=(0.0D0,0.0D0)         !initialize the layer matrix 
	do i=1,2*n
	  layer1(i,i)=(1.0D0,0.0D0)
	end do 
    call ST0(b1,kox,koy,ko,n,S0,T0)    !S0,T0. air layer matrix

TMsl:do i=1,m+2                           !m+2 slices every layer
      select case(i) 
	  case(1,m+2)
	    h=a*(0.433012702D0-f)
	    c1=0.0D0
		c2=0.0D0
        call construc(a1,c1,c2,b1,kox,koy,ko,die1,die2,n,t1,p)    !t1,p matrix

		call STa(ko,t1,p,beta,sa,ta,n)   !Sa,Ta,Beta       

		call s_matrix(h,n,S0,T0,sa,ta,beta,s)          !Si

		call s_layer1(n,s,layer1)             !the whole Sn for layer1

	  case(2:m+1)
	    theta_i=(i-1.5D0)*pi/m
		h=f*a*(dcos(theta_i-0.5D0*pi/m)-dcos(theta_i+0.5D0*pi/m))
		c1=2.0D0*f*a*dsin(theta_i)
        c2=0.0D0
		call construc(a1,c1,c2,b1,kox,koy,ko,die1,die2,n,t1,p)    !p matrix
        
		call STa(ko,t1,p,beta,sa,ta,n)   !Sa,Ta,Beta
	    
		call s_matrix(h,n,S0,T0,sa,ta,beta,s)          !s-matrix
        
		call s_layer1(n,s,layer1)             !the whole s_matrix for layer1
	  case default
	  write(*,*)'the process of discretizing layer faults!'  
      stop
	  end select
	end do TMsl
     !layer2 matrix through the symmetry transformation
    y01=0.5D0*a1
	call s_layer2(n,layer1,b1,y01,layer2)
    
	call s_layer1(n,layer2,layer1)              !now,layer1 is unit cell layer. 
    !slab scattering matrix
	slab=layer1
	do i=2,m_unit
	  call s_layer1(n,layer1,slab)
	end do 
	!reflective and transmission spectra
	RT=1      !reflective spectra
	call Ez_array(ko,kox,koy,b1,n,nn,EF0,slab,T0,EFn,Ez,RT)
    RT_result(k0,2)=RT_fun(ko,kox,koy,koz,b1,n,nn,EFn,Ez)
	RT=2
	call Ez_array(ko,kox,koy,b1,n,nn,EF0,slab,T0,EFn,Ez,RT)
	RT_result(k0,3)=RT_fun(ko,kox,koy,koz,b1,n,nn,EFn,Ez)
	end do TMsn
	!the precision control
    temp=0.0D0
	sum=0.0D0
	do k0=1,2
	  do i=1,1500
	    temp=temp+dabs(RT_result(i,k0+1)-temp_result(i,k0))
		sum=sum+dabs(RT_result(i,k0+1))
	  end do
	end do
	temp_result(1:1500,1:2)=RT_result(1:1500,2:3)
    deallocate(p,sa,ta,S0,T0,t1,beta)
	deallocate(s,layer1,layer2,slab)
	deallocate(EF0,EFn,Ez)
	  temp=temp/sum
	  write(*,*)'TM:the orders of diffration wave=',(nn-1)/2
	  write(*,*)'TM:the relative precision=',temp
	  n=n+4
	  nn=n/2
	  end do TMsw

    open(unit=1,file='TM_R&T_spectra.dat',status='replace')
    write(1,*)'The TM direction!'
	write(1,*)'The correlative paremeter:'
	write(1,*)'the lattice constant(cm)=',a
	write(1,*)'the fraction filling(r/a)=',f
	write(1,*)'the number of unit cell=',m_unit 
	write(1,*)'the dielectric constants of periodic structure/matrix(slab)=',die1,'  /',die2
    write(1,*)'the incident wave vector angle:Polar/Azimuthal angle(radian)=',theta_k,'  /',phi_k
    write(1,*)'the  direction of electics field:Polar/Azimuthal angle(radian)=',theta_E,'  /',phi_E
	write(1,*)'the used plane wave numbers (orders of diffraction waves)=',(nn-3)/2
	write(1,*)'the range of wave number(cm-1):*************************************'
	do k0=1,1500
	  write(1,*)RT_result(k0,1)
	end do
	
	write(1,*)'the reflective spectra**********************************************'
	do k0=1,1500
 	  write(1,*)RT_result(k0,2)
	end do
	
	write(1,*)'the transmission spectra********************************************'
	do k0=1,1500
	write(1,*)RT_result(k0,3)
	end do
	close(1)
   
   else  !*****************************************************************************
TMbw:do while(temp>error)
	  allocate(p(n,n),sa(n,n),ta(n,n),S0(n,n),T0(n,n),t1(n,n),beta(n,n))
	  allocate(s(2*n,2*n),layer1(2*n,2*n),layer2(2*n,2*n),slab(2*n,2*n))
	  allocate(EF0(n),EFn(n),Ez(nn))
      EF0=(0.0D0,0.0D0)
	  EF0(nn)=dcmplx(Eox)
	  EF0(nn+1)=dcmplx(Eoy)
TMbn:do k0=1,1500                 !******************** the range of wavenumber!!!!!!!!!
        ko=20.0D0*pi*k0	         
     	RT_result(k0,1)=ko              
		kox=ko*dsin(theta_k)*dcos(phi_k) 
       	koy=ko*dsin(theta_k)*dsin(phi_k) 
	    koz=ko*dcos(theta_k) 	
	  layer1=(0.0D0,0.0D0)         !initialize the layer matrix 
	do i=1,2*n
	  layer1(i,i)=(1.0D0,0.0D0)
	end do 
    call ST0(b1,kox,koy,ko,n,S0,T0)    !S0,T0. air layer matrix

       sh_i_1=0.0D0
TMbl:do i=1,m                           !m+2 slices every layer
	    call overlap(i,a,f,m,h,sh_i_1,c1,c2,1)
	    
        call construc(a1,c1,c2,b1,kox,koy,ko,die1,die2,n,t1,p)    !p matrix
        
		call STa(ko,t1,p,beta,sa,ta,n)   !Sa,Ta,Beta
	    
		call s_matrix(h,n,S0,T0,sa,ta,beta,s)          !s-matrix
        
		call s_layer1(n,s,layer1)             !the whole s_matrix for layer1
	end do TMbl
     !layer2 matrix through the symmetry transformation
    y01=0.5D0*a1
	call s_layer2(n,layer1,b1,y01,layer2)
    
	call s_layer1(n,layer2,layer1)              !now,layer1 is unit cell layer. 
    !slab scattering matrix
	slab=layer1
	do i=2,m_unit
	  call s_layer1(n,layer1,slab)
	end do 
	!reflective and transmission spectra
	RT=1      !reflective spectra
	call Ez_array(ko,kox,koy,b1,n,nn,EF0,slab,T0,EFn,Ez,RT)
    RT_result(k0,2)=RT_fun(ko,kox,koy,koz,b1,n,nn,EFn,Ez)
	RT=2
	call Ez_array(ko,kox,koy,b1,n,nn,EF0,slab,T0,EFn,Ez,RT)
	RT_result(k0,3)=RT_fun(ko,kox,koy,koz,b1,n,nn,EFn,Ez)
	end do TMbn
	!the precision control
    temp=0.0D0
	sum=0.0D0
	do k0=1,2
	  do i=1,1500
	    temp=temp+dabs(RT_result(i,k0+1)-temp_result(i,k0))
		sum=sum+dabs(RT_result(i,k0+1))
	  end do
	end do
	temp_result(1:1500,1:2)=RT_result(1:1500,2:3)
    deallocate(p,sa,ta,S0,T0,t1,beta)
	deallocate(s,layer1,layer2,slab)
	deallocate(EF0,EFn,Ez)
	  temp=temp/sum
	 write(*,*)'TM:the order numbers of diffrction waves=',(nn-1)/2
 	 write(*,*)'TM:the relative precision=',temp  
	  n=n+4
	  nn=n/2 
	  end do TMbw

    open(unit=1,file='TM_R&T_spectra.dat',status='replace')
    write(1,*)'The TM direction!'
	write(1,*)'The correlative paremeter:'
	write(1,*)'the lattice constant(cm)=',a 
	write(1,*)'the fraction filling(r/a)=',f
	write(1,*)'the number of unit cell=',m_unit 
	write(1,*)'the dielectric constants of periodic structures/matrix(slab)=',die1,'  /',die2
    write(1,*)'the incident wave vector angle:Polar/Azimuthal angle(radian)=',theta_k,'  /',phi_k
    write(1,*)'the direction of electrics field:Polar/Azimuthal angle(radian)=',theta_E,'  /',phi_E
	write(1,*)'the used plane wave numbers (orders of diffraction waves)=',(nn-3)/2
	write(1,*)'the range of wave number(cm-1):*************************************'
	do k0=1,1500
	  write(1,*)RT_result(k0,1)
	end do
	
	write(1,*)'the reflective spectra**********************************************'
	do k0=1,1500
 	  write(1,*)RT_result(k0,2)
	end do
	
	write(1,*)'the transmission spectra********************************************'
	do k0=1,1500
	write(1,*)RT_result(k0,3)
	end do
	close(1)
   end if TM    !**********************************************************************
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!***********************************************************************************
! the TK direction!
!***********************************************************************************
! the initialization!

    n=n-8
	nn=n/2
	temp_result=0.0D0 
	temp=1.0D0
	a2=sqrt(3.0D0)*a
	b2=2.0D0*pi/(a2)
TK:   if(f<=0.25D0)then  !*******************************************************
TKsw:do while(temp>error)
	  allocate(p(n,n),sa(n,n),ta(n,n),S0(n,n),T0(n,n),t1(n,n),beta(n,n))
	  allocate(s(2*n,2*n),layer1(2*n,2*n),layer2(2*n,2*n),slab(2*n,2*n))
	  allocate(EF0(n),EFn(n),Ez(nn))
      EF0=(0.0D0,0.0D0)
	  EF0(nn)=dcmplx(Eox)
	  EF0(nn+1)=dcmplx(Eoy)
TKsn:do k0=1,1500             !**************************** the range of wavenumber!!!!
        ko=20.0D0*pi*k0	         
     	RT_result(k0,1)=ko              
		kox=ko*dsin(theta_k)*dcos(phi_k) 
       	koy=ko*dsin(theta_k)*dsin(phi_k) 
	    koz=ko*dcos(theta_k) 	
	  layer1=(0.0D0,0.0D0)         !initialize the layer matrix 
	do i=1,2*n
	  layer1(i,i)=(1.0D0,0.0D0)
	end do 
    call ST0(b2,kox,koy,ko,n,S0,T0)    !S0,T0. air layer matrix

TKsl:do i=1,m+2                           !m+2 slices every layer
      select case(i) 
	  case(1,m+2)
	   h=a*(0.25D0-f)
	    c1=0.0D0
		c2=0.0D0
        call construc(a2,c1,c2,b2,kox,koy,ko,die1,die2,n,t1,p)    !t1,p matrix

		call STa(ko,t1,p,beta,sa,ta,n)   !Sa,Ta,Beta       

		call s_matrix(h,n,S0,T0,sa,ta,beta,s)          !Si

		call s_layer1(n,s,layer1)             !the whole Sn for layer1

	  case(2:m+1)
	    theta_i=(i-1.5D0)*pi/m
		h=f*a*(dcos(theta_i-0.5D0*pi/m)-dcos(theta_i+0.5D0*pi/m))
		c1=2.0D0*f*a*dsin(theta_i)
        c2=0.0D0
		call construc(a2,c1,c2,b2,kox,koy,ko,die1,die2,n,t1,p)    !p matrix
        
		call STa(ko,t1,p,beta,sa,ta,n)   !Sa,Ta,Beta
	    
		call s_matrix(h,n,S0,T0,sa,ta,beta,s)          !s-matrix
        
		call s_layer1(n,s,layer1)             !the whole s_matrix for layer1
	  case default
	  write(*,*)'the process of discretizing layer faults!'  
      stop
	  end select
	end do TKsl
     !layer2 matrix through the symmetry transformation
    y02=a2/2.0D0
	call s_layer2(n,layer1,b2,y02,layer2)
    
	call s_layer1(n,layer2,layer1)              !now,layer1 is unit cell layer. 
    !slab scattering matrix
	slab=layer1
	do i=2,m_unit
	  call s_layer1(n,layer1,slab)
	end do 
	!reflective and transmission spectra
	RT=1      !reflective spectra
	call Ez_array(ko,kox,koy,b2,n,nn,EF0,slab,T0,EFn,Ez,RT)
    RT_result(k0,2)=RT_fun(ko,kox,koy,koz,b2,n,nn,EFn,Ez)
	RT=2
	call Ez_array(ko,kox,koy,b2,n,nn,EF0,slab,T0,EFn,Ez,RT)
	RT_result(k0,3)=RT_fun(ko,kox,koy,koz,b2,n,nn,EFn,Ez)
	end do TKsn
	!the precision control
    temp=0.0D0
	sum=0.0D0
	do k0=1,2
	  do i=1,1500
	    temp=temp+dabs(RT_result(i,k0+1)-temp_result(i,k0))
		sum=sum+dabs(RT_result(i,k0+1))
	  end do
	end do
	temp_result(1:1500,1:2)=RT_result(1:1500,2:3)
    deallocate(p,sa,ta,S0,T0,t1,beta)
	deallocate(s,layer1,layer2,slab)
	deallocate(EF0,EFn,Ez)
	  temp=temp/sum
	  write(*,*)'TK:the orders of diffration wave=',(nn-1)/2
	  write(*,*)'TK:the relative precision=',temp
	  n=n+4
	  nn=n/2
	  end do TKsw

    open(unit=1,file='TK_R&T_spectra.dat',status='replace')
    write(1,*)'The TK direction!'
    write(1,*)'The correlative paremeter:'
	write(1,*)'the lattice constant(cm)=',a 
	write(1,*)'the fraction filling(r/a)=',f
	write(1,*)'the number of unit cell=',m_unit 
	write(1,*)'the dielectric constants of periodic structure/matrix(slab)=',die1,'  /',die2
    write(1,*)'the incident wave vector angle:Polar/Azimuthal angle(radian)=',theta_k,'  /',phi_k
    write(1,*)'the direction of electrics field:Polar/Azimuthal angle(radian)=',theta_E,'  /',phi_E
	write(1,*)'the used plane wave numbers (orders of diffraction waves)=',(nn-3)/2
	write(1,*)'the range of wave number(cm-1):*************************************'
	do k0=1,1500
	  write(1,*)RT_result(k0,1)
	end do
	
	write(1,*)'the reflective spectra**********************************************'
	do k0=1,1500
 	  write(1,*)RT_result(k0,2)
	end do
	
	write(1,*)'the transmission spectra********************************************'
	do k0=1,1500
	write(1,*)RT_result(k0,3)
	end do
	close(1)
   
   else  !*****************************************************************************
TKbw:do while(temp>error)