C    
C   
C     LINESQ.FORT(EVEN) Programmed by K.Sakoda (1997)    
C   
C     Line Defect Mode in 2D Square Lattice of Circular Cylinders     
C        E Polarization     
C           Even Parity   
C    
C   
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q) 
C     
      COMMON /C1/ QEZ1(41,321), QEZ2(41,321)    
      COMMON /C3/ QHX(41,321), QHY(41,321)      
      COMMON /C4/ QEPS(41,321)      
C     
C     
      DATA AK1,AK2,KPOINT / 0.01D0, 0.01D0, 1 /
      DATA OM1,OM2,NPOINT / 0.34D0, 0.38D0,  80/
      DATA NX, NY, LY, NT, NC / 40, 152, 4,  320,  50 /      
      DATA A, B, R   / 1.27, 1.016, 0.2550/     
      DATA QEPS1  / ( 10.0D0, 0.0D0) /     
      DATA QEPS2  / ( 1.0D0, 0.0D0)  /    
      DATA X0, Y0 / 0.0D0, 0.2D0 /  
      DATA EMU    / 1.0D0   / 
C     K1 : K Minimum in unit of 2pi/A
C     K2 : K Maximum in unit of 2pi/A
C     KPOINT : Number of K
C     OM1: Omega Minimum
C     OM2: Omega Maximum
C     NPOINT : Number of Omega
C     NX : Number of Points in X direction
C     NY : Number of Points in Y direction
C     NT : Number of Points in a Period
C     NC : Number of Cycles
C     LY : Number of Unit Cell in Y direction
C     A  : Lattice Constant
C     B  : Center of the 1st Lattice
C     R  : Radius of Rods
C     QEPS1 : Dielectric Constant of the Rod
C     QEPS2 : Dielectric Constant of Background
C     X0,Y0 : Position of Dipole Moment
C     EMU   : Magnitude of Dipole
C     
      DATA PI / 3.14159265D0 /
C
      OPEN(4,FILE='line.dat')
C     
      CALL REC1(A,B,R,QEPS1,QEPS2,NX,NY,LY,X0,Y0,EMU,NT,NC) 
C     
      CALL CALEPS(A, B, R, QEPS1, QEPS2, NX, NY, LY)  
C     
C     
      W = B + A*(DFLOAT(LY)-1.0D0)  
      L0 = (X0/A+0.5D0)*DFLOAT(NX)  
      M0 = Y0/W*DFLOAT(NY)    
      DX = 2.0*PI/DFLOAT(NX)  
      DY = 2.0*PI/DFLOAT(NY)*W/A    
C    
C     
      DO 600 J1=1, KPOINT     
      AK = AK1 + J1*(AK2-AK1)/KPOINT      
C     
      DO 500 I1=0, NPOINT     
      OMEGA = OM1 + I1*(OM2-OM1)/NPOINT   
      DT  = 2.0*PI/OMEGA/DFLOAT(NT)/2.0   
      CALL INIRAD(NX, NY)     
C     
      DO 400 I2=1, NC   
C     
      CALL CALRAD(AK,OMEGA,NT,L0,M0,DX,DY,DT,EMU,NX,NY)     
C     CALL CALU(U,NX,NY,DX,DY,AK,OMEGA)  
C     CALL REC2(AK,OMEGA,U)  
C     
  400 CONTINUE    
C     
      CALL CALU(U,NX,NY,DX,DY,AK,OMEGA)   
      CALL REC2(AK,OMEGA,U)   
C     CALL REC3(NX,NY)  
C     CALL REC4(NX,NY,A,W)    
C     
  500 CONTINUE    
C     
  600 CONTINUE    
C     
C    
   70 STOP  
      END   
C     
C     
C     
C  ****************************************************************
C    
C     
      SUBROUTINE REC4(NX,NY,A,W)   
C     
      IMPLICIT REAL*8(A,W), COMPLEX*16 (Q)     
C   
      COMMON /C1/ QEZ1(41,321), QEZ2(41,321)    
      COMMON /C3/ QHX(41,321), QHY(41,321)     
      COMMON /C4/ QEPS(41,321)     
      DIMENSION EZ(41,321)  
C   
C   
      S = 0.0    
C    
      DO 20 J=1, NY+1  
      DO 10 I=1, NX+1  
      EZ(I,J) = DREAL(QEZ1(I,J))   
      IF(ABS(EZ(I,J)).GT.S) S = EZ(I,J)   
   10 CONTINUE  
   20 CONTINUE   
C   
C     
      DO 100 J=-NY, 0, 2     
      Y = DFLOAT(J)/DFLOAT(NY)*W/A      
      DO 90 I=0, NX, 2  
      X = FLOAT(I)/FLOAT(NX) - 0.5 
      B = EZ(I+1,-J+1)/S     
      WRITE(4,600) X, Y, B   
   90 CONTINUE  
  100 CONTINUE   
C    
      DO 200 J=2, NY, 2      
      Y = DFLOAT(J)/DFLOAT(NY)*W/A      
      DO 190 I=0, NX, 2      
      X = FLOAT(I)/FLOAT(NX) - 0.5 
      B = EZ(I+1,J+1)/S      
      WRITE(4,600) X, Y, B   
  190 CONTINUE   
  200 CONTINUE    
C    
  600 FORMAT(1H ,F8.4,5H,    ,F8.4,5H,    ,F8.4)      
C   
C    
      RETURN     
      END  
C    
C    
C    
C  ****************************************************************
C     
C     
      SUBROUTINE REC3(NX,NY)  
C     
      IMPLICIT COMPLEX*16 (Q) 
C     
      COMMON /C1/ QEZ1(41,321), QEZ2(41,321)    
      COMMON /C3/ QHX(41,321), QHY(41,321)      
      COMMON /C4/ QEPS(41,321)      
C     
C     
      I = NX/2    
C     
C     
      DO 100 J=0, NY    
      A = DREAL(QEZ1(I+1,J+1))      
      B = DIMAG(QEZ1(I+1,J+1))      
      WRITE(4,600) I, J, A, B 
  100 CONTINUE    
C     
  600 FORMAT(1H ,I5,5H,    ,I5,5H,    ,E15.7,5H,    ,E15.7) 
C     
C     
      RETURN      
      END   
C     
C     
C     
C  *****************************************************************
C     
C     
      SUBROUTINE REC1(A,B,R,QEPS1,QEPS2,NX,NY,LY,X0,Y0,EMU,NT,NC) 
C     
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)   
C     
      WRITE(4,603) A, B, R    
  603 FORMAT(1H ,E15.7,5X,E15.7,5X,E15.7) 
C     
      WRITE(4,600) QEPS1      
      WRITE(4,600) QEPS2      
  600 FORMAT(1H ,E15.7,5X,E15.7)    
C     
      WRITE(4,601) NX, NY, LY 
  601 FORMAT(1H ,2(I5,5X),I5) 
C     
      WRITE(4,603) X0, Y0, EMU      
C     
      WRITE(4,602) NT, NC     
  602 FORMAT(1H ,I5,5X,I5)    
C     
C     
      RETURN      
      END   
C     
C     
C     
C  *****************************************************************
C     
C     
      SUBROUTINE REC2(AK,OMEGA,U)   
C     
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)   
C     
C     
      WRITE(4,600) AK, OMEGA, U     
  600 FORMAT(1H ,E15.7,5X,E15.7,5X,E15.7) 
C     
      RETURN      
      END   
C     
C     
C     
C  ****************************************************************
C     
C     
      SUBROUTINE CALRAD(AK,OMEGA,NT,L0,M0,DX,DY,DT,EMU,NX,NY)     
C     
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)   
C     
      COMMON /C1/ QEZ1(41,321), QEZ2(41,321)    
      COMMON /C4/ QEPS(41,321)      
      DATA QI, PI / (0.0D0, 1.0D0), 3.14159265D0 /    
      DATA Q2     / (2.0D0, 0.0D0) /      
C     
C     
      QA = DCMPLX( DT*DT/DX/DX, 0.0D0 )   
      QB = DCMPLX( DT*DT/DY/DY, 0.0D0 )   
      QP = CDEXP(2.0D0*PI*AK*QI)    
C     
C     
      DO 300 I=1, NT    
C     
C     L=0, M=0    
C     
      QEZ2(1,1) = Q2*QEZ1(1,1) - QEZ2(1,1)      
     1 +( QA*(QEZ1(2,1)+QEZ1(NX,1)/QP-Q2*QEZ1(1,1))   
     1   +QB*(QEZ1(1,2)+QEZ1(1,2)-Q2*QEZ1(1,1)))/QEPS(1,1)  
C     
C     L=NX, M=0   
C     
      QEZ2(NX+1,1) = QEZ2(1,1)*QP   
C     
C     M=0   
C     
      DO 10 L=1, NX-1   
      QEZ2(L+1,1) = Q2*QEZ1(L+1,1) - QEZ2(L+1,1)      
     1 +( QA*(QEZ1(L+2,1)+QEZ1(L,1)-Q2*QEZ1(L+1,1))   
     1   +QB*(QEZ1(L+1,2)+QEZ1(L+1,2)-Q2*QEZ1(L+1,1)))/QEPS(L+1,1)      
   10 CONTINUE    
C     
C     L=0, M=NY   
C     
      QEZ2(1,NY+1) = Q2*QEZ1(1,NY+1) - QEZ2(1,NY+1)   
     1 +( QA*(QEZ1(2,NY+1)+QEZ1(NX,NY+1)/QP-Q2*QEZ1(1,NY+1))      
     1   +QB*(QEZ1(1,NY)+QEZ1(1,NY)-Q2*QEZ1(1,NY+1)))/QEPS(1,NY+1) 
C     
C     L=NX, M=NY  
C     
      QEZ2(NX+1,NY+1) = QEZ2(1,NY+1)*QP   
C     
C     M=NY  
C     
      DO 20 L=1, NX-1   
      QEZ2(L+1,NY+1) = Q2*QEZ1(L+1,NY+1) - QEZ2(L+1,NY+1)   
     1 +( QA*(QEZ1(L+2,NY+1)+QEZ1(L,NY+1)-Q2*QEZ1(L+1,NY+1))      
     1   +QB*(QEZ1(L+1,NY)+QEZ1(L+1,NY)-Q2*QEZ1(L+1,NY+1))) 
     1   /QEPS(L+1,NY+1)      
   20 CONTINUE    
C     
C     L=0, L=NX   
C     
      DO 25 M=1, NY-1   
      QEZ2(1,M+1) = Q2*QEZ1(1,M+1) - QEZ2(1,M+1)      
     1 +( QA*(QEZ1(2,M+1)+QEZ1(NX,M+1)/QP-Q2*QEZ1(1,M+1))   
     1   +QB*(QEZ1(1,M+2)+QEZ1(1,M)-Q2*QEZ1(1,M+1)))/QEPS(1,M+1)  
C     
      QEZ2(NX+1,M+1) = QEZ2(1,M+1)*QP     
   25 CONTINUE    
C     
C     
      DO 40 M=1, NY-1   
      DO 30 L=1, NX-1   
      QEZ2(L+1,M+1) = Q2*QEZ1(L+1,M+1) - QEZ2(L+1,M+1)      
     1 +( QA*(QEZ1(L+2,M+1)+QEZ1(L,M+1)-Q2*QEZ1(L+1,M+1))   
     1   +QB*(QEZ1(L+1,M+2)+QEZ1(L+1,M)-Q2*QEZ1(L+1,M+1)))/QEPS(L+1,M+1)
   30 CONTINUE    
   40 CONTINUE    
C     
C     
      T = DT*(I-1)*2.0  
      QQ = 4.0*PI*EMU*OMEGA*OMEGA*DT*DT*CDEXP(-QI*(OMEGA*T+PI/2.0))
      QEZ2(L0+1,M0+1)=QEZ2(L0+1,M0+1) + QQ/QEPS(L0+1,M0+1)  
C     
C     
C     
C     L=0, M=0    
C     
      QEZ1(1,1) = Q2*QEZ2(1,1) - QEZ1(1,1)      
     1 +( QA*(QEZ2(2,1)+QEZ2(NX,1)/QP-Q2*QEZ2(1,1))   
     1   +QB*(QEZ2(1,2)+QEZ2(1,2)-Q2*QEZ2(1,1)))/QEPS(1,1)  
C     
C     L=NX, M=0   
C     
      QEZ1(NX+1,1) = QEZ1(1,1)*QP   
C     
C     M=0   
C     
      DO 110 L=1, NX-1  
      QEZ1(L+1,1) = Q2*QEZ2(L+1,1) - QEZ1(L+1,1)      
     1 +( QA*(QEZ2(L+2,1)+QEZ2(L,1)-Q2*QEZ2(L+1,1))   
     1   +QB*(QEZ2(L+1,2)+QEZ2(L+1,2)-Q2*QEZ2(L+1,1)))/QEPS(L+1,1)      
  110 CONTINUE    
C     
C     L=0, M=NY   
C     
      QEZ1(1,NY+1) = Q2*QEZ2(1,NY+1) - QEZ1(1,NY+1)   
     1 +( QA*(QEZ2(2,NY+1)+QEZ2(NX,NY+1)/QP-Q2*QEZ2(1,NY+1))      
     1   +QB*(QEZ2(1,NY)+QEZ2(1,NY)-Q2*QEZ2(1,NY+1)))/QEPS(1,NY+1)      
C     
C     L=NX, M=NY  
C     
      QEZ1(NX+1,NY+1) = QEZ1(1,NY+1)*QP   
C     
C     M=NY  
C     
      DO 120 L=1, NX-1  
      QEZ1(L+1,NY+1) = Q2*QEZ2(L+1,NY+1) - QEZ1(L+1,NY+1)   
     1 +( QA*(QEZ2(L+2,NY+1)+QEZ2(L,NY+1)-Q2*QEZ2(L+1,NY+1))      
     1   +QB*(QEZ2(L+1,NY)+QEZ2(L+1,NY)-Q2*QEZ2(L+1,NY+1))) 
     1   /QEPS(L+1,NY+1)      
  120 CONTINUE    
C     
C     L=0, L=NX   
C     
      DO 125 M=1, NY-1  
      QEZ1(1,M+1) = Q2*QEZ2(1,M+1) - QEZ1(1,M+1)      
     1 +( QA*(QEZ2(2,M+1)+QEZ2(NX,M+1)/QP-Q2*QEZ2(1,M+1))   
     1   +QB*(QEZ2(1,M+2)+QEZ2(1,M)-Q2*QEZ2(1,M+1)))/QEPS(1,M+1)  
C     
      QEZ1(NX+1,M+1) = QEZ1(1,M+1)*QP     
  125 CONTINUE    
C     
C     
      DO 140 M=1, NY-1  
      DO 130 L=1, NX-1  
      QEZ1(L+1,M+1) = Q2*QEZ2(L+1,M+1) - QEZ1(L+1,M+1)      
     1 +( QA*(QEZ2(L+2,M+1)+QEZ2(L,M+1)-Q2*QEZ2(L+1,M+1))   
     1   +QB*(QEZ2(L+1,M+2)+QEZ2(L+1,M)-Q2*QEZ2(L+1,M+1)))/QEPS(L+1,M+1)
  130 CONTINUE    
  140 CONTINUE    
C     
C     
      T = DT*(I*2-1)    
      QQ = 4.0*PI*EMU*OMEGA*OMEGA*DT*DT*CDEXP(-QI*(OMEGA*T+PI/2.0))
      QEZ1(L0+1,M0+1)=QEZ1(L0+1,M0+1) + QQ/QEPS(L0+1,M0+1)  
C     
C     
  300 CONTINUE    
C     
C     
      RETURN      
      END   
C     
C     
C     
C  ****************************************************************
C     
C     
      SUBROUTINE CALU(U,NX,NY,DX,DY,AK,OMEGA)   
C     
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)   
C     
      COMMON /C1/ QEZ1(41,321), QEZ2(41,321)    
      COMMON /C3/ QHX(41,321), QHY(41,321)      
      COMMON /C4/ QEPS(41,321)      
      DATA QI, PI / (0.0D0, 1.0D0), 3.14159265D0 /    
C     
C     
      QX = QI*OMEGA*DY  
      QY = QI*OMEGA*DX  
      QP = CDEXP(2.0D0*PI*AK*QI)    
C     
C     I=NX+1      
C     
      DO 10 J=1, NY     
      QHX(NX+1,J) = (QEZ1(NX+1,J+1)-QEZ1(NX+1,J))/QX  
      QHY(NX+1,J) = (QEZ1(NX+1,J)-QEZ1(2,J)*QP)/QY    
   10 CONTINUE    
C     
C     
C     J=NY+1      
C     
      DO 20 I=1, NX     
      QHX(I,NY+1) = (QEZ1(I,NY)-QEZ1(I,NY+1))/QX      
      QHY(I,NY+1) = (QEZ1(I,NY+1)-QEZ1(I+1,NY+1))/QY  
   20 CONTINUE    
C     
C     
C     I=NX+1, J=NY+1    
C     
      QHX(NX+1,NY+1) = (QEZ1(NX+1,NY)-QEZ1(NX+1,NY+1))/QX   
      QHY(NX+1,NY+1) = (QEZ1(NX+1,NY+1)-QEZ1(2,NY+1)*QP)/QY 
C     
C     
      DO 100 J=1, NY    
      DO 90 I=1, NX     
      QHX(I,J) = (QEZ1(I,J+1)-QEZ1(I,J))/QX     
      QHY(I,J) = (QEZ1(I,J)-QEZ1(I+1,J))/QY     
   90 CONTINUE    
  100 CONTINUE    
C     
C     
C     
      S4 = DREAL(QEPS(1,1))*(DREAL(QEZ1(1,1))**2+DIMAG(QEZ1(1,1))**2)   
      S4 = S4 + DREAL(QEPS(1,NY+1)) 
     1     *(DREAL(QEZ1(1,NY+1))**2+DIMAG(QEZ1(1,NY+1))**2) 
      S4 = S4 + DREAL(QEPS(NX+1,1)) 
     1     *(DREAL(QEZ1(NX+1,1))**2+DIMAG(QEZ1(NX+1,1))**2) 
      S4 = S4 + DREAL(QEPS(NX+1,NY+1))    
     1     *(DREAL(QEZ1(NX+1,NY+1))**2+DIMAG(QEZ1(NX+1,NY+1))**2) 
C     
      SX4 = DREAL(QHX(1,1))**2+DIMAG(QHX(1,1))**2     
      SX4 = SX4 + DREAL(QHX(1,NY+1))**2+DIMAG(QHX(1,NY+1))**2     
      SX4 = SX4 + DREAL(QHX(NX+1,1))**2+DIMAG(QHX(NX+1,1))**2     
      SX4 = SX4 + DREAL(QHX(NX+1,NY+1))**2+DIMAG(QHX(NX+1,NY+1))**2
C     
      SY4 = DREAL(QHY(1,1))**2+DIMAG(QHY(1,1))**2     
      SY4 = SY4 + DREAL(QHY(1,NY+1))**2+DIMAG(QHY(1,NY+1))**2     
      SY4 = SY4 + DREAL(QHY(NX+1,1))**2+DIMAG(QHY(NX+1,1))**2     
      SY4 = SY4 + DREAL(QHY(NX+1,NY+1))**2+DIMAG(QHY(NX+1,NY+1))**2
C     
C     
      S2 = 0.0D0  
      SX2 = 0.0D0 
      SY2 = 0.0D0 
C     
      DO 320 I=2, NY    
      S2 = S2 + DREAL(QEPS(1,I))*(DREAL(QEZ1(1,I))**2 
     1      +DIMAG(QEZ1(1,I))**2)   
     1  + DREAL(QEPS(NX+1,I))*(DREAL(QEZ1(NX+1,I))**2 
     1      +DIMAG(QEZ1(NX+1,I))**2)      
      SX2 = SX2 + DREAL(QHX(1,I))**2+DIMAG(QHX(1,I))**2     
     1    + DREAL(QHX(NX+1,I))**2+DIMAG(QHX(NX+1,I))**2     
      SY2 = SY2 + DREAL(QHY(1,I))**2+DIMAG(QHY(1,I))**2     
     1    + DREAL(QHY(NX+1,I))**2+DIMAG(QHY(NX+1,I))**2     
  320 CONTINUE    
C     
C     
      DO 330 I=2, NX    
      S2 = S2 + DREAL(QEPS(I,1))*(DREAL(QEZ1(I,1))**2 
     1      +DIMAG(QEZ1(I,1))**2)   
     1  + DREAL(QEPS(I,NY+1))*(DREAL(QEZ1(I,NY+1))**2 
     1      +DIMAG(QEZ1(I,NY+1))**2)      
      SX2 = SX2 + DREAL(QHX(I,1))**2+DIMAG(QHX(I,1))**2     
     1    + DREAL(QHX(I,NY+1))**2+DIMAG(QHX(I,NY+1))**2     
      SY2 = SY2 + DREAL(QHY(I,1))**2+DIMAG(QHY(I,1))**2     
     1    + DREAL(QHY(I,NY+1))**2+DIMAG(QHY(I,NY+1))**2     
  330 CONTINUE    
C     
C     
      S1 = 0.0D0  
      SX1 = 0.0D0 
      SY1 = 0.0D0 
C     
      DO 360 J=2, NY    
      DO 350 I=2, NX    
      S1 = S1 + DREAL(QEPS(I,J))*(DREAL(QEZ1(I,J))**2 
     1      +DIMAG(QEZ1(I,J))**2)   
      SX1 = SX1 + DREAL(QHX(I,J))**2+DIMAG(QHX(I,J))**2     
      SY1 = SY1 + DREAL(QHY(I,J))**2+DIMAG(QHY(I,J))**2     
  350 CONTINUE    
  360 CONTINUE    
C     
C     
      U = (S1*4.0 + S2*2.0 + S4)/NX/NY/32.0/PI  
      UX = (SX1*4.0 + SX2*2.0 + SX4)/NX/NY/32.0/PI    
      UY = (SY1*4.0 + SY2*2.0 + SY4)/NX/NY/32.0/PI    
C     
      U = U + UX + UY   
C     
C     
      RETURN      
      END   
C     
C     
C  ***************************************************************
C     
C     
      SUBROUTINE INIRAD(NX, NY)     
C     
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)   
C     
      COMMON /C1/ QEZ1(41,321), QEZ2(41,321)    
C     
C     
      DO 100 J=1, NY+1  
      DO  90 I=1, NX+1  
      QEZ1(I,J) = DCMPLX(0.0D0, 0.0D0)    
      QEZ2(I,J) = DCMPLX(0.0D0, 0.0D0)    
   90 CONTINUE    
  100 CONTINUE    
C     
C     
      RETURN      
      END   
C     
C     
C  ****************************************************************
C     
C     
      SUBROUTINE CALEPS(A, B, R, QEPS1, QEPS2, NX, NY, LY)  
C     
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)   
C     
      COMMON /C4/ QEPS(41,321)      
C     
C     
      R2 = R*R    
      W = B + A*(DFLOAT(LY)-1.0D0)  
C     
C     
      DO 100 J=0, NY    
      Y = DFLOAT(J)/DFLOAT(NY) * W  
C     
      DO 90  I=0, NX    
      X = (DFLOAT(I)/DFLOAT(NX) - 0.5D0) * A    
C     
      DO 80 L=1, LY     
      YL = B + A * (DFLOAT(L)-1.0D0)      
      IF(X*X+(Y-YL)*(Y-YL).LT.R2) GO TO 50      
   80 CONTINUE    
C     
      QEPS(I+1,J+1) = QEPS2   
      GO TO 90    
C     
   50 QEPS(I+1,J+1) = QEPS1   
C     
   90 CONTINUE    
  100 CONTINUE    
C     
C     
      RETURN      
      END   
