tmm.f90

用于分析2d介质柱阵列模型中电磁场分布 附TXT使用说明

	  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)
TKbn: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

       sh_i_1=0.0D0
TKbl:do i=1,m                           !m+2 slices every layer
	    call overlap(i,a,f,m,h,sh_i_1,c1,c2,2)
	    
        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
	end do TKbl
     !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 TKbn
	!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 order numbers of diffrction waves=',(nn-1)/2
 	 write(*,*)'TK:the relative precision=',temp  
	  n=n+4
	  nn=n/2
	  end do TKbw

    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 structures/matrix(slab)=',die1,'  /',die2
    write(1,*)'the incident wave vector angle:Polar/Azimuthal angle(radian)=',theta_k,'  /',phi_k
    write(1,*)'the direcion 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 TK 
	end program TMM
	

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    subroutine overlap(i,a,f,m,h,sh_i_1,c1,c2,KM)
	implicit none
	real(8),parameter::pi=3.141592653589793D0
	integer,intent(in)::i,m,KM
	real(8),intent(in)::a,f
	real(8),intent(out)::h,c1,c2
	real(8),intent(inout)::sh_i_1
	real(8)::alpha,theta_i,L,psi,area0,area
	intrinsic dacos,sqrt,dcos,dsin
	  select case(KM)
	  case(1)            !TM direction!
	    alpha=dacos(sqrt(3.0D0)/(4.0D0*f))
	  case(2)            !TK direction!      
        alpha=dacos(1.0D0/(4.0D0*f))
	  end select
	  theta_i=alpha+(pi-2.0D0*alpha)/(2.0D0*m)+(i-1)*(pi-2.0D0*alpha)/m
	  h=f*a*(dcos(theta_i-0.5D0*(pi-2.0D0*alpha)/m)-dcos(theta_i+0.5D0*(pi-2.0D0*alpha)/m))
	  c1=2.0D0*f*a*dsin(theta_i)
        !c2
	  select case(KM)
	  case(1)            !TM direction!
	    L=sqrt(3.0D0)*a/4.0D0+sh_i_1+h
	  case(2)
	    L=a/4.0D0+sh_i_1+h
      end select
	    if(L<=f*a)then
		  L=L-0.5D0*h
		  c2=2.0D0*sqrt((f*a)**2-L**2)
		else 
		  if((L-h)<f*a)then
		    L=L-h
		    psi=2.0D0*dacos(L/(f*a))
		    area0=0.5D0*psi*((f*a)**2)
		    area=L*sqrt((f*a)**2-L**2)
		    c2=(area0-area)/h
		  else
            c2=0.0D0
          end if
		end if
	  sh_i_1=sh_i_1+h
	  return
    end subroutine overlap
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    function RT_fun(ko,kox,koy,koz,b,n,nn,EFn,Ez)
	implicit none
    integer,intent(in)::n,nn
	real(8),intent(in)::ko,kox,koy,koz,b
	complex(8),intent(in)::EFn(n),Ez(nn)
	real(8)::beta,temp,tempk
	real(8)::RT_fun
	integer::i,j,k
	intrinsic dsqrt,cdabs,dabs
	j=1
	k=-(nn-1)/2	
	RT_fun=0.0D0
	do i=1,n,2
	  tempk=ko**2-kox**2-(koy+k*b)**2
	  if(tempk>=0.0D0)then
	    beta=dsqrt(tempk)
        temp=(cdabs(EFn(i))**2+cdabs(EFn(i+1))**2+cdabs(Ez(j))**2)*dabs(beta)/dabs(koz)
	    RT_fun=RT_fun+temp
	  end if
		j=j+1
		k=k+1
	end do
	return
	end function RT_fun
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  subroutine Ez_array(ko,kox,koy,b,n,nn,EF0,slab,T0,EFn,Ez,RT)
	implicit none
	integer,intent(in)::n,nn,RT
	real(8),intent(in)::ko,kox,koy,b
	complex(8),intent(in)::EF0(n),slab(2*n,2*n),T0(n,n)
	complex(8),intent(out)::EFn(n),Ez(nn)
	integer::i,j,k
    real(8)::ky
	intrinsic matmul
    j=1
	k=-(nn-1)/2
	select case(RT)
	case(1)          !1 for reflective index
	EFn=matmul(slab(n+1:2*n,1:n),EF0)    !Enx,Eny
    do i=1,n,2
	  ky=koy+k*b
	  Ez(j)=EFn(i)*(kox*T0(i+1,i)-ky*T0(i,i))/ko+EFn(i+1)*(kox*T0(i+1,i+1)-ky*T0(i,i+1))/ko
      j=j+1
	  k=k+1
	end do
	case(2) 
	EFn=matmul(slab(1:n,1:n),EF0)         !2 for transnission index
    do i=1,n,2
	  ky=koy+k*b
	  Ez(j)=EFn(i)*(ky*T0(i,i)-kox*T0(i+1,i))/ko+EFn(i+1)*(ky*T0(i,i+1)-kox*T0(i+1,i+1))/ko
	  j=j+1
      k=k+1
    end do
	end select
	return
	end subroutine Ez_array
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    subroutine s_layer2(n,layer1,b,y0,layer2)   !symmetry transformation 
    implicit none
	integer,intent(in)::n
	real(8),intent(in)::b,y0
	complex(8),intent(in)::layer1(2*n,2*n)
	complex(8),intent(out)::layer2(2*n,2*n)
	complex(8)::temp
	integer::i,j,k,l
    intrinsic cdexp
	k=-(n/2-1)/2
	do i=1,n,2
	  l=-(n/2-1)/2
	  do j=1,n,2
	    temp=cdexp((0.0D0,1.0D0)*(k-l)*b*y0)
	    layer2(i:i+1,j:j+1)=temp*layer1(i:i+1,j:j+1)

		layer2(i:i+1,n+j:n+j+1)=temp*layer1(i:i+1,n+j:n+j+1)
		
		layer2(n+i:n+i+1,j:j+1)=temp*layer1(n+i:n+i+1,j:j+1)
		
		layer2(n+i:n+i+1,n+j:n+j+1)=temp*layer1(n+i:n+i+1,n+j:n+j+1)
		l=l+1
	  end do
	    k=k+1
	end do
	return
	end subroutine s_layer2
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    subroutine s_layer1(n,s,layer)   !iteration matrix 
	implicit none
	integer,intent(in)::n	
	complex(8),intent(in)::s(2*n,2*n)
    complex(8),intent(inout)::layer(2*n,2*n)
	complex(8)::temp(n,n),atemp(n,n),temp11(n,n),temp12(n,n),temp21(n,n),temp22(n,n),unit(n,n)
	integer::i,itemp,jtemp
	intrinsic matmul
	external dlincg
      unit=(0.0D0,0.0D0)
	do i=1,n
	  unit(i,i)=(1.0D0,0.0D0)
	end do
	  temp11=layer(1:n,1:n)         !layer11
	  temp12=layer(1:n,n+1:2*n)     !layer12 
	  temp21=layer(n+1:2*n,1:n)     !layer21
	  temp22=layer(n+1:2*n,n+1:2*n) !layer22
	  temp=unit-matmul(temp12,s(n+1:2*n,1:n))
	call dlincg(n,temp,n,atemp,n)
      layer(1:n,1:n)=matmul(matmul(s(1:n,1:n),atemp),temp11)
	  temp=unit-matmul(s(n+1:2*n,1:n),temp12)
	call dlincg(n,temp,n,atemp,n)
	  layer(1:n,n+1:2*n)=s(1:n,n+1:2*n)+matmul(matmul(s(1:n,1:n),temp12),matmul(atemp,s(n+1:2*n,n+1:2*n)))
      temp=unit-matmul(temp12,s(n+1:2*n,1:n))
	call dlincg(n,temp,n,atemp,n)
	  layer(n+1:2*n,1:n)=temp21+matmul(matmul(temp22,s(n+1:2*n,1:n)),matmul(atemp,temp11))
	  temp=unit-matmul(s(n+1:2*n,1:n),temp12)
	call dlincg(n,temp,n,atemp,n)
	  layer(n+1:2*n,n+1:2*n)=matmul(matmul(temp22,atemp),s(n+1:2*n,n+1:2*n))
	  return
	end subroutine s_layer1
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    subroutine s_matrix(h,n,S0,T0,sa,ta,beta,s)  
	implicit none
	integer,intent(in)::n
	real(8),intent(in)::h
	complex(8),intent(in)::S0(n,n),T0(n,n),sa(n,n),ta(n,n)
	complex(8),intent(out)::s(2*n,2*n) 
	complex(8),intent(inout)::beta(n,n)
	complex(8)::asa(n,n),ata(n,n),a11(n,n),aa11(n,n),a12(n,n),p1(n,n),p2(n,n),t1(n,n),t2(n,n),temp(n,n)
    integer::i
	intrinsic matmul,cdexp
	external dlincg
	  call dlincg(n,sa,n,asa,n) !!!!!!!!!!!!!!!!!!!!!!!!用于确定SA?
	  call dlincg(n,ta,n,ata,n)
	  a11=0.5D0*(matmul(asa,S0)+matmul(ata,T0))
	  a12=0.5D0*(matmul(asa,S0)-matmul(ata,T0))
	  call dlincg(n,a11,n,aa11,n)   !a11-1
      do i=1,n
	    beta(i,i)=cdexp((0.0D0,1.0D0)*beta(i,i)*h)
	  end do
	    temp=a11-matmul(matmul(matmul(beta,a12),matmul(aa11,beta)),a12)
	  call dlincg(n,temp,n,p1,n)    !p1
		p2=matmul(matmul(aa11,beta),matmul(a12,p1))
        t1=matmul(beta,a11)
		t2=-a12	
       s(1:n,1:n)=matmul(p1,t1)+matmul(p2,t2)        !S11
	   s(n+1:2*n,1:n)=matmul(p1,t2)+matmul(p2,t1)    !S21
       s(1:n,n+1:2*n)=s(n+1:2*n,1:n)                 !S12
	   s(n+1:2*n,n+1:2*n)=s(1:n,1:n)                 !S22
	  return
	end subroutine s_matrix
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    subroutine ST0(b,kox,koy,ko,n,S0,T0)   !S0,T0
	implicit none
	integer,intent(in)::n
	real(8),intent(in)::b,ko,kox,koy
	complex(8),intent(out)::S0(n,n),T0(n,n)
	real(8)::tempk,ky
	complex(8)::beta0
	integer::i,k
	intrinsic dsqrt
	S0=(0.0D0,0.0D0)
	do i=1,n
	  S0(i,i)=(1.0D0,0.0D0)
	end do
	T0=(0.0D0,0.0D0)
	k=-(n/2-1)/2
	do i=1,n,2
	  ky=koy+k*b
	  tempk=ko**2-kox**2-ky**2
      if(tempk>=0.0D0)then
	    beta0=(1.0D0,0.0D0)*dsqrt(tempk)
      else 
	    beta0=(0.0D0,1.0D0)*dsqrt(-tempk)
	  end if
	  T0(i,i)=-kox*ky/(ko*beta0)
	  T0(i,i+1)=(kox**2-ko**2)/(ko*beta0)
	  T0(i+1,i)=(ko**2-ky**2)/(ko*beta0)
	  T0(i+1,i+1)=(kox*ky)/(ko*beta0)
	  k=k+1	 
	end do
	return
	end subroutine ST0
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    subroutine STa(ko,t1,p,beta,sa,ta,n)  !make Sa,Ta,beta
    implicit none
   	integer,intent(in)::n
	real(8),intent(in)::ko,t1(n,n),p(n,n)
    complex(8),intent(out)::beta(n,n),sa(n,n),ta(n,n)
	real(8)::at1(n,n)
	complex(8)::lambda(n)
	integer::i
    intrinsic cdsqrt,matmul
	external devcrg,dlinrg
	  call dlinrg(n,t1,n,at1,n)
  	  call devcrg(n,p,n,lambda,sa,n)
 	  beta=(0.0D0,0.0D0)
	  do i=1,n
	    beta(i,i)=cdsqrt(-lambda(i)) 
	  end do
      ta=ko*matmul(matmul(at1,sa),beta)
	  return
	end subroutine STa
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    subroutine construc(a,c1,c2,b,kox,koy,ko,die1,die2,n,t1,p) !T1,P
    implicit none
	real(8),intent(in)::a,c1,c2,b,kox,koy,ko,die1,die2
	integer,intent(in)::n
	real(8),intent(out)::t1(n,n),p(n,n)  
    real(8)::ft,aft,ftrans,delta,deltat,t2(n,n)   !ft fourier transform
	integer::i,j,k,l
	intrinsic matmul
	external ftrans,delta
	k=-(n/2-1)/2
	do i=1,n,2
	  l=-(n/2-1)/2
	  do j=1,n,2
	       !t1 matrix
		aft=ftrans(a,b,c1,c2,die1,die2,k,l,2)
		deltat=delta(k,l)
	    t1(i,j)=kox*(koy+l*b)*aft
        t1(i,j+1)=ko*ko*deltat-kox*kox*aft
		t1(i+1,j)=(koy+k*b)*(koy+l*b)*aft-deltat*ko*ko
        t1(i+1,j+1)=-kox*(koy+k*b)*aft
           !t2 matrix
		ft=ftrans(a,b,c1,c2,die1,die2,k,l,1)
        t2(i,j)=-kox*(koy+l*b)*deltat
		t2(i,j+1)=kox*kox*deltat-ko*ko*ft
		t2(i+1,j)=ko*ko*ft-(koy+k*b)*(koy+l*b)*deltat
		t2(i+1,j+1)=kox*(koy+k*b)*deltat
	      l=l+1
	  end do
          k=k+1
	end do
    p=-(1.0D0/(ko**2))*matmul(t1,t2)
	return
	end subroutine construc
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    function delta(i,j)
	implicit none
	integer,intent(in)::i,j
	real(8)::delta
	if(i==j)then
	delta=1.0D0
	else
	delta=0.0D0
	end if
	return
	end function delta 
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    function ftrans(a,b,c1,c2,die1,die2,k,l,n)
	implicit none
	real(8),parameter::pi=3.141592653589793D0
	real(8),intent(in)::a,b,c1,c2,die1,die2
	integer,intent(in)::k,l,n
	real(8)::ftrans
	integer::m
	real(8)::g
	intrinsic dsin
	m=k-l
	g=m*b
	select case(n)
	case(1)      !dielectric function
	  if(m==0)then
    	ftrans=die2+(die1-die2)*(c1+c2)/a
      else
	    ftrans=(die1-die2)*(dsin(c1*g/2.0D0)+dsin(c2*g/2.0D0)*((-1)**m))/(m*pi)
	  end if
	case(2)      !the inverse of dielectric function
      if(m==0)then
	    ftrans=1.0D0/die2+(1.0D0/die1-1.0D0/die2)*(c1+c2)/a
	  else
	    ftrans=(1.0D0/die1-1.0D0/die2)*(dsin(c1*g/2.0D0)+dsin(c2*g/2.0D0)*((-1)**m))/(m*pi)
	  end if
	case default
	write(*,*)'the Fouriertransformation is error!'
	stop
	end select
    end function ftrans