simple_new.f90

TMM模拟电磁波在堆膜中传播过程 和场分布,可简易改变参数 作不同体系的透射反射谱图和场分布图

program main	 
    implicit none
       real(8),parameter::pi=3.141592653589793D0,c=3.0D8
       real(8)::d0=15.0,d1=10.0,index2=1.0d0,index0=1.0d0,index3=1.5d0,index1=1.5d0
	   real(8)::wavelengh,w=0.0D0 
	   real(8)::R,T        
	   integer::i
	   integer::layernumber=5
	   complex(8)::k0,k1,k2
       complex(8),dimension(2,2)::D01,D12,P0,P1,M,D23,D30,D10,D21

       intrinsic matmul
do i=1,1000,1
       w=w+(10.0*c)/(index1*d1*500)
       k0=index0*w/c
	   k1=index1*w/c
!	   k2=index2*w/c

    	D01(1,1)=0.5*(1+index0/index1)
		D01(1,2)=0.5*(1-index0/index1)
		D01(2,1)=0.5*(1-index0/index1)
		D01(2,2)=0.5*(1+index0/index1)

    	D12(1,1)=0.5*(1+index1/index2)
		D12(1,2)=0.5*(1-index1/index2)
		D12(2,1)=0.5*(1-index1/index2)
		D12(2,2)=0.5*(1+index1/index2)

		D23(1,1)=0.5*(1+index2/index3)
		D23(1,2)=0.5*(1-index2/index3)
		D23(2,1)=0.5*(1-index2/index3)
		D23(2,2)=0.5*(1+index2/index3)

    	D10(1,1)=0.5*(1+index3/index0)
		D10(1,2)=0.5*(1-index3/index0)
		D10(2,1)=0.5*(1-index3/index0)
		D10(2,2)=0.5*(1+index3/index0)

        P1(1,1)=cdexp((0.0d0,1.0d0)*k1*d1)
        P1(1,2)=(0.0,0.0)
        P1(2,1)=(0.0,0.0)
        P1(2,2)=cdexp((0.0d0,-1.0d0)*k1*d1)

        P0(1,1)=cdexp((0.0d0,1.0d0)*k0*d0)
        P0(1,2)=(0.0,0.0)
        P0(2,1)=(0.0,0.0)
        P0(2,2)=cdexp((0.0d0,-1.0d0)*k0*d0)
       
       M=matmul(D10,P1)
	   M=matmul(M,D01)
	   M=matmul(M,P0)
	   M=matmul(M,D10)
	   M=matmul(M,P1)
	   M=matmul(M,D01)
	   M=matmul(M,P0)
	   M=matmul(M,D10)
	   M=matmul(M,P1)
	   M=matmul(M,D01)
	   M=matmul(M,P0)
	   M=matmul(M,D10)
	   M=matmul(M,P1)
	   M=matmul(M,D01)
	   M=matmul(M,P0)
	   M=matmul(M,D10)
	   M=matmul(M,P1)
	   	   M=matmul(M,D01)
	   M=matmul(M,P0)
	   M=matmul(M,D10)
	   M=matmul(M,P1)
	   	   M=matmul(M,D01)
	   M=matmul(M,P0)
	   M=matmul(M,D10)
	   M=matmul(M,P1)
	   	   M=matmul(M,D01)
	   M=matmul(M,P0)
	   M=matmul(M,D10)
	   M=matmul(M,P1)
	   	   M=matmul(M,D01)
	   M=matmul(M,P0)
	   M=matmul(M,D10)
	   M=matmul(M,P1)
	   	   M=matmul(M,D01)
	   M=matmul(M,P0)
	   M=matmul(M,D10)
	   M=matmul(M,P1)
	   M=matmul(M,D01)
	   


	   T=abs(1.0/M(1,1))
	   T=T**2

    open(unit=1,file='w--T.txt',status='replace')
	write(1,10)w,T

10	   format(1x,e18.6,e18.6)

end do
end program