test2.f90

在光子晶体中

!plane wave method for 2d photonic crystals
!正方形排列、正方形截面
use msimsl
implicit none
real*8, parameter:: pb=0.68d0,pi=3.14159265d0,angle=30.0d0     !lb=0.45d0,fb is the volume factor ratio 
integer,parameter:: nn=625 
integer,parameter:: nx=12 , ny=12, in=40 !the total number of plane wave is nx*ny
real*8, parameter:: eb=1.0D+0,em=12.96d0,ub=1.0D+0,um=12.96d0 ! m-matrix, e-element, ro: mass density 
real*8:: kgpx,kgpy,kgx,kgy,lb ! kgp: k+Gp; kg: k+G 
real*8, dimension(nn,nn):: Txy,TZ,eg,egi,txyg,tzg,ug,ugi,KTxyg,KTzg
 	 ! Txy=|k+G|*|k+G'| use to caculate TE mode
	 ! TZ =(k+G).(k+G')	use to caculate TM mode
integer:: n1x,n2x,n1y,n2y,n1z,n2z,i,j,k,kk ! G: (n1x,n1y)a/2pi; G prime: (n2x,n2y)a/2pi
integer:: Gx,Gy,Gz  ! Gx,Gy,Gz: G minus G prime
integer:: ii,ik ! ii: the variable describing wave vector k
real*8:: kx,ky ! kx,ky: complements of k value along certain direction
doublecomplex:: evalr(nn), evalzr(nn) ! eigen value: eval/evalr
open(1,file='625air_die_square4_f68_r30tm.DAT',STATUS='UNKNOWN')
open(2,file='625air_die_square4_f68_r30te.dat',status='unknown')
open(3,file='625air_die_square4_f68_r30tmxm.DAT',STATUS='UNKNOWN')
open(4,file='625air_die_square4_f68_r30texm.dat',status='unknown')

!!!!!!!!!compute elastic contants!!!!!!!!!!
lb=dsqrt(pb)/2

!!!!!!!!!!!!!!!!! give the value to the eigenmatrix !!!!!!!!!!!!!!!!!!!!!!!!!
!************* define the direction of the wave vector K*********************
do ik=1,3		   !1:T->M; 2:T->X; 3: X->M
do ii=0,in 
  
if(ik==1) then   ! ik=1 means the wave vector along the TM direction
 kx=ii/(2.0d0*in)-0.5d0       ! k is defined as k*a/2pi
 ky=ii/(2.0d0*in)-0.5d0
 
elseif(ik==2) then   ! ik=2 means the wave vector along the TX direction
 
 kx=ii/(2.0d0*in)
 ky=0.0d0
else
 kx=0.5d0
 ky=ii/(2.0d0*in)

endif

 
!*********************** end of the wave vector K direction difinition ***************
i=1
j=1
k=1																				
kk=1
Txy=0.0d0
Txyg=0.0d0
Tz=0.0d0
Tzg=0.0d0
eg=0.0d0
egi=0.0d0
do n1x=-nx, nx
   kgx=kx+n1x			   ! (k+G) at x-direction

  do n1y=-ny,ny

	 kgy=ky+n1y
    ! summation of G; G is defined as (n1x,n1y,n1z)*a/2pi
	 
	  do n2x=-nx,nx
	     Gx=n1x-n2x  !Gx mean G minus G prime  G-G' 
	     kgpx=kx+n2x         ! (k+G') at x-direction
		do n2y=-ny,ny
		   ! summation of G prime

		   	      kgpy=ky+n2y
				  Gy=n1y-n2y  !Gy mean G minus G prime 
				   
			  				  
!!!!!!!!!!!!!!!!!!!!! calculate the Furior coefficient of ro,cij!!!!!!!				  
				  eG(kk,k)=Fur(em,eb,gx,gy,pb,pi)
                  uG(kk,k)=Fur(um,ub,gx,gy,pb,pi)		!(edit)		  
                  !!!!!!!!!!!!!!!!!!!  equation u1 !!!!!!!!!!!!!!!!!!!!!!!	  
				  Tz(kk,k)=kgpx*kgx+kgpy*kgy		! (k+G).(k+G')  use to caculate TM mode
                  !Txy(kk,k)=dsqrt(kgx*kgx+kgy*kgy)*dsqrt(kgpx*kgpx+kgpy*kgpy)  ! use to caculate TE mode	!(edit)
                  !!!!!!!!!!!!!!!!!!!  equation u1 !!!!!!!!!!!!!!!!!!!!!!!	  
				  
				  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
			k=k+1	
            j=j+1	     ! the column number of Txy should be changed when G prime changed
          
		  enddo
		enddo
	           ! end of the do-loop of G prime
      
	  j=1
	  k=1
	  kk=kk+1
	  i=i+1          ! the row number of T should be changed when G changed
 
	
  enddo
enddo                ! end of the do-loop of G
CALL DLINRG (nn,eg,nn,egi,nn)		 !矩阵求逆
CALL DLINRG (nn,ug,nn,ugi,nn)		 !(edit)
Txyg=tz*ugi						 !(edit)
tzg=tz*egi
KTzg=MATMUL(ugi,tzg)
KTxyg=MATMUL(egi,txyg)

call DEvlrg(nn,KTxyg,nn,evalr)		 ! TE mode
call DEvlrg(nn,KTzg,nn,evalzr)		 ! TM mode

!if(ik==3) then
!write(1,*) 0.5+ky,(dble(cdsqrt(evalzr(i))),i=nn,1)
!write(2,*) 0.5+ky,(dble(cdsqrt(evalr(i))),i=nn,1)

!else
!write(1,*) kx,(dble(cdsqrt(evalzr(i))),i=nn,1)
!write(2,*) kx,(dble(cdsqrt(evalr(i))),i=nn,1)
!endif





do i=nn,nn-8,-1
 if(ik==3) then
 write(3,*) 0.5+ky, dble(cdsqrt(evalzr(i)))
 write(4,*) 0.5+ky, dble(cdsqrt(evalr(i)))
 else
 write(1,*) kx,dble(cdsqrt(evalzr(i)))
 write(2,*) kx,dble(cdsqrt(evalr(i))) ! output the eigenvalue
 endif
enddo

!do i=nn,nn-8,-1
 !if(ik==3) then
 !write(4,*) 0.5+ky, dble(cdsqrt(evalr(i)))
 !else
 !write(2,*) kx,dble(cdsqrt(evalr(i))) ! output the eigenvalue
 !endif
!enddo

enddo  ! end the do-loop of wave vector K
enddo

contains
function Fur(cm,cb,ginx,giny,pb,pi)
double precision:: Fur,cb,cm,pb,pbg,pi,gnx,gny,gnxp,gnyp,seta
integer::          ginx,giny,ij,ijk,il
  
  seta=pi*angle/180.0
  gnxp=2*pi*ginx
  gnyp=2*pi*giny
  gnx=gnxp*dcos(seta)+gnyp*dsin(seta)
  gny=-gnxp*dsin(seta)+gnyp*dcos(seta)
  if(abs(gnx)<=0.00000010d0.and.abs(gny)<=0.00000010d0) then
    Fur=cb*pb+cm*(1-pb)
  elseif(abs(gny)>0.00000010d0.and.abs(gnx)<=0.00000010d0) then
    pbg=pb*dsin(gny*lb)/(gny*lb)
	Fur=(cb-cm)*pbg
	  elseif(abs(gnx)>0.00000010d0.and.abs(gny)<=0.00000010d0) then
        pbg=pb*dsin(gnx*lb)/(gnx*lb)
	    Fur=(cb-cm)*pbg
         else
           pbg=pb*dsin(gnx*lb)*dsin(gny*lb)/(gnx*lb*gny*lb)
	       Fur=(cb-cm)*pbg
       
  endif

  return
end function Fur

end