square_1.f90

计算二维光子晶体正方空气柱的源代码

! TM,TE both
PROGRAM  MAIN


external SGEEV,SGETRF,FUNC

!  integer::num=11           !  num is odd  and num**2=the number of plane wave 

integer::LDVL,LDVR,LWORK, IPIV, INFO 

CHARACTER:: JOBVL,JOBVR

REAL::WR,WI,WORK,VL,VR

DIMENSION::WR_TE(1:441),WR_TM(1:441),WI(1:441),vl(1:441),vr(1:441)

DIMENSION::WORK(1:1764)

REAL::Kz,ky,k 

REAL::la,L1,L2     !la for lattice constant

REAL::EPa,EPb,PI

REAL::TE(441,441),TM(441,441),EPC ,eig_TE,EIG_TM

DIMENSION EPC(0:20,0:20),eig_TE(1:10),EIG_TM(1:10)

INTEGER I,J,M1,M2,N1,N2,Q,P1
real::KG1,KG2

REAL::F,P,x ,c   !F=P**2,P=RA/LA

   PI=4.0*ATAN(1.0)
   EPa=1.0            
   EPb=12.96            
   L1=SQRT(50.0)
   L2=L1
   la=10.0
   
   LDVL=484
   LDVR=484
   LWORK=1674

   JOBVL="N"
   JOBVR="N"


   open(unit=11,file='s_TE_1.txt',status='replace')
   open(unit=12,file='s_TM_1.txt',status='replace')


!******求介电函数的傅立叶展开******************************

P=L1/LA
F=P*P
DO J=0,20
  DO I=0,20    
	if(i==0.AND.J==0)then 
	  EPC(i,j)=F/EPA+(1-F)/EPB	  
	   else 	     
          EPC(I,J)=(1.0/EPa-1.0/EPB)*FUNc(I,L1)*FUNC(J,L2)      
      end if 
   END DO 
END DO

!********ky,kz****

DO Q=0,150,2
  if(Q<=50) THEN 
    KY=REAL(Q)/100
      Kz=0.0
       ELSE IF(Q>50.AND.Q<=100) THEN     
	  KY=0.5 
     KZ=REAL(Q-50)/100
    ELSE 
   KY=REAL(150-Q)/100
  KZ=KY
END IF 
!********calculate the matrix element**********************
DO J=1,441
   DO I=1,441
     M1=(I-1)/21-10
       N1=i-11-21*((I-1)/21)
         M2=(J-1)/21-10
           N2=J-11-21*((J-1)/21)
          KG1=(KY+M1)*(KY+M1)+(KZ+N1)*(KZ+N1)
   	    KG2=(KY+M2)*(KY+M2)+(KZ+N2)*(KZ+N2)
      TE(I,J)=EPC(ABS(M1-M2),ABS(N1-N2))*SQRT(KG1*KG2)
	 TM(I,J)=EPC(ABS(M1-M2),ABS(N1-N2))*((Ky+M1)*(Ky+M2)+(Kz+N1)*(Kz+N2))
  end do
end do

 !求M的本征值

!   call SGETRF(121,121,A,121,IPIV,INFO)
     call SGEEV(JOBVL,JOBVR,441,TE,441,WR_TE,WI,VL,LDVL,VR,LDVR,WORK,1764,INFO)
	 call SGEEV(JOBVL,JOBVR,441,TM,441,WR_TM,WI,VL,LDVL,VR,LDVR,WORK,1764,INFO)


!**********对本征值重新排序***********************
!TE
   
   do i=1,441
    do j=i,441
	 	  if(WR_TE(J)>=0.0.AND.(wr_TE(i)-wr_TE(j))>=1.0e-6) then 
	    c = wr_TE(j)
		wr_TE(j)=wr_TE(i)
		wr_TE(i)=c
	 end if 
	end do
 end do 
! TM
   do i=1,441
    do j=i,441
	 	  if(WR_TM(J)>=0.0.AND.(wr_TM(i)-wr_TM(j))>=1.0e-6) then 
	    c = wr_TM(j)
		wr_TM(j)=wr_TM(i)
		wr_TM(i)=c
	 end if 
	end do
 end do 
!***************************************
do i=1,10
    eig_TE(i)=sqrt(wr_TE(i))
	eig_TM(i)=sqrt(wr_TM(i))
   write(*,*) eig_TE(i),EIG_TM(I)                
end do 
!*******************
 if(q<=50) then 
  write(11,"(f7.3,10f8.5)") ky,eig_TE
   write(12,"(f7.3,10f8.5)") ky,eig_TM 
!************************************************
 else if(q>50.and.q<=100) then 
      k=0.5+kz
      write(11,"(f7.3,10f8.5)") k,eig_TE
	   write(12,"(f7.3,10f8.5)") k,eig_TM 
!*************************************************
  else if (q>100)then  
  k=sqrt((ky-0.5)**2+(kz-0.5)**2)+1.0
    write(11,"(f7.3,10f8.5)") k,eig_TE
	write(12,"(f7.3,10f8.5)") k,eig_TM 
 !***********************************************
 end if 
 
END DO 

 stop

END 




FUNCTION FUNC(I,L1)
REAL::FUNC
REAL::L1,P
INTEGER::I
REAL::LA=10.0
REAL::PI

PI=4.0*ATAN(1.0)
P=L1/LA

IF(I==0) THEN 
   FUNC=p
    ELSE 
   FUNC=SIN(I*PI*P)/(I*PI)
 END IF 
 RETURN 
 END FUNCTION FUNC