C
C
C     DEFSQ.FORT(A1) Programmed by K.Sakoda (1996) 
C      
C     Defect Mode in 2D Square Lattice of Circular Cylinders  
C         Totally Symmetric (A1) Mode Calculated on 1/8 Lattice        
C       
C       
      IMPLICIT REAL*8 (A-H,O-Z)    
C       
      COMMON /C1/ QEZ1(241,241), QEZ2(241,241)       
      COMMON /C4/ QEPS(241,241)    
C       
C       
      DATA OM1,OM2,NPOINT / 0.45D0, 0.48D0, 60 /  
      DATA NX, LX, NT, NC / 40, 4,  80,  50 /        
      DATA A, R   / 1.27D0, 0.48D0 /        
      DATA QEPS1  / 9.0D0 /        
      DATA QEPS2  / 1.0D0 /        
      DATA QEPSD  / 1.0D0 /        
      DATA X0, Y0 / 0.00D0, 0.0D0 /         
      DATA EMU    / 1.0D0   /      
C     OM1: Omega Minimum  
C     OM2: Omega Maximum  
C     NPOINT : Number of Omega     
C     NX : Number of Points in X direction  
C     NT : Number of Points in a Period     
C     NC : Number of Cycles        
C     LX : Number of Unit Cell in X direction        
C     A  : Lattice Constant        
C     R  : Radius of a Rod         
C     QEPS1 : Dielectric Constant of the Rod         
C     QEPS2 : Dielectric Constant of Background      
C     QEPSD : Dielectric Constant of Defect 
C     X0,Y0 : Position of Dipole Moment     
C     EMU   : Magnitude of Dipole  
C       
      DATA PI / 3.141592653589D0 / 
C
      OPEN(4,FILE='sa1.dat')
C       
      CALL REC1(A,R,QEPS1,QEPS2,QEPSD,NX,LX,X0,Y0,NT,NC)      
      CALL CALEPS(A, R, QEPS1, QEPS2, QEPSD, NX, LX) 
C       
C       
      L0 = X0/A*NX        
      M0 = Y0/A*NX        
      DX = 2.0*PI/DFLOAT(NX)       
C       
      DO 500 I1=0, NPOINT 
      OMEGA = OM1 + I1*(OM2-OM1)/NPOINT     
      DT = PI/OMEGA/DFLOAT(NT)     
      CALL INIRAD(NX, LX) 
C       
      DO 400 I2=1, NC     
C       
      CALL CALRAD(OMEGA,NT,L0,M0,DX,DT,EMU,NX,LX)    
C     CALL CALU(U,NX,LX,DX,OMEGA)  
C     CALL REC2(OMEGA, U)        
C       
  400 CONTINUE   
C       
      CALL CALU(U,NX,LX,DX,OMEGA)  
      CALL REC2(OMEGA, U)        
C     CALL REC3(NX,LX)    
C     CALL REC4(NX,LX)    
C       
C     REC1 : Input Data   
C     REC2 : Omega vs Energy       
C     REC3 : Field Distribution    
C       
  500 CONTINUE   
C       
C       
C     GO TO 10   
C       
C       
   70 STOP       
      END        
C       
C       
C  *****************************************************************   
C       
C       
      SUBROUTINE REC4(NX,LX)       
C       
      IMPLICIT REAL*8 (B-H,O-Z)    
C       
      COMMON /C1/ QEZ1(241,241), QEZ2(241,241)       
      DIMENSION AE(241,241)        
C       
      DATA N1 /  4 /      
C       
C       
      DO 20 J=2, NX*LX+1  
      DO 10 I=1, J-1      
C     QEZ1(I,J) = QEZ1(J,I)        
      AE(I,J) = SNGL(QEZ1(J,I))    
      AE(J,I) = AE(I,J)   
   10 CONTINUE   
   20 CONTINUE   
C       
C       
      DO 30 I=1, NX*LX+1  
      AE(I,I) = SNGL(QEZ1(I,I))    
   30 CONTINUE   
C       
C       
      A1 = 0.0   
      DO 50 J=1, NX*LX+1  
      DO 40 I=1, NX*LX+1  
      IF(ABS(AE(I,J)).GT.ABS(A1)) A1 = AE(I,J)       
   40 CONTINUE   
   50 CONTINUE   
C       
C       
      DO 110 J=-NX*LX, NX*LX, N1   
      AJ = FLOAT(J)/NX    
      JP = IABS(J)        
C       
      DO 100 I=-NX*LX, NX*LX, N1   
      AI = FLOAT(I)/NX    
      IP = IABS(I)        
C       
      WRITE(4,600) AI, AJ, AE(IP+1,JP+1)/A1 
  100 CONTINUE   
  110 CONTINUE   
C       
  600 FORMAT(1H ,F8.4,1H,,F8.4,1H,,F8.4)    
C       
C       
      RETURN     
      END        
C       
C       
C  *****************************************************************   
C       
C       
      SUBROUTINE REC3(NX,LX)       
C       
      IMPLICIT REAL*8 (A-H,O-Z)    
C       
      COMMON /C1/ QEZ1(241,241), QEZ2(241,241)       
C       
C       
      J = 0      
      DO 100 I=0, NX*LX   
      WRITE(4,600) I, J, QEZ1(I+1,J+1)      
  100 CONTINUE   
C       
  600 FORMAT(1H ,I5,5X,I5,5X,E15.7)         
C       
C       
      RETURN     
      END        
C       
C       
C       
C  *****************************************************************   
C       
C       
      SUBROUTINE REC1(A,R,QEPS1,QEPS2,QEPSD,NX,LX,X0,Y0,NT,NC)         
C       
      IMPLICIT REAL*8 (A-H,O-Z)    
C       
      WRITE(4,600) A, R   
  600 FORMAT(1H ,E15.7,5X,E15.7)   
      WRITE(4,600) QEPS1  
      WRITE(4,600) QEPS2  
      WRITE(4,600) QEPSD  
      WRITE(4,601) NX, LX 
  601 FORMAT(1H ,I5,5X,I5)         
      WRITE(4,600) X0, Y0 
      WRITE(4,602) NT, NC 
  602 FORMAT(1H ,I5,5X,I5)         
C       
C       
      RETURN     
      END        
C       
C       
C       
C  *****************************************************************   
C       
C       
      SUBROUTINE REC2(OMEGA, U)  
C       
      IMPLICIT REAL*8 (A-H,O-Z)    
C       
C       
      WRITE(4,600) OMEGA, U    
  600 FORMAT(1H ,E15.7,5X,E15.7)      
C       
      RETURN     
      END        
C       
C       
C       
C  *****************************************************************   
C       
C       
      SUBROUTINE CALRAD(OMEGA,NT,L0,M0,DX,DT,EMU,NX,LX)       
C       
      IMPLICIT REAL*8 (A-H,O-Z)    
C       
      COMMON /C1/ QEZ1(241,241), QEZ2(241,241)       
      COMMON /C4/ QEPS(241,241)    
      DIMENSION EPS(241,241)       
      DATA PI / 3.141592653589D0 / 
C       
C       
      Q1 = 4.0*PI*EMU*OMEGA*OMEGA*DT*DT/QEPS(L0+1,M0+1)       
      A = DT*DT/DX/DX     
      N1 = NX*LX 
C       
C       
      DO 410 M=1, N1+1    
      DO 400 L=M, N1+1    
      EPS(L,M) = QEPS(L,M)/A       
  400 CONTINUE   
  410 CONTINUE   
C       
C       
      DO 300 I=1, NT      
C       
C     L=0, M=0   
C       
      QEZ2(1,1) = 2.0*QEZ1(1,1) - QEZ2(1,1) 
     1 + 4.0*( QEZ1(2,1) - QEZ1(1,1) )/ EPS(1,1)     
C       
C     M=0        
C       
      DO 10 L=1, N1-1     
      QEZ2(L+1,1) = 2.0*QEZ1(L+1,1) - QEZ2(L+1,1)    
     1 +( QEZ1(L+2,1)+QEZ1(L,1)-4.0*QEZ1(L+1,1)      
     1   +2.0*QEZ1(L+1,2) )/ EPS(L+1,1)     
   10 CONTINUE   
C       
C     L=N1, M=0  
C       
      QEZ2(N1+1,1) = 2.0*QEZ1(N1+1,1) - QEZ2(N1+1,1) 
     1 +( 2.0*QEZ1(N1,1)-4.0*QEZ1(N1+1,1)   
     1   +2.0*QEZ1(N1+1,2) )/ EPS(N1+1,1)   
C       
C     L=N1       
C       
      DO 15 M=1, N1-1     
      QEZ2(N1+1,M+1) = 2.0*QEZ1(N1+1,M+1) - QEZ2(N1+1,M+1)    
     1 +( 2.0*QEZ1(N1,M+1)-4.0*QEZ1(N1+1,M+1)        
     1   +QEZ1(N1+1,M+2)+QEZ1(N1+1,M) )/ EPS(N1+1,M+1)        
   15 CONTINUE   
C       
C     L=N1, M=N1 
C       
      QEZ2(N1+1,N1+1) = 2.0*QEZ1(N1+1,N1+1) - QEZ2(N1+1,N1+1) 
     1 +4.0*( QEZ1(N1+1,N1) - QEZ1(N1+1,N1+1) )/ EPS(N1+1,N1+1)        
C       
C     L=M        
C       
      DO 20 L=1, N1-1     
      QEZ2(L+1,L+1) = 2.0*QEZ1(L+1,L+1) - QEZ2(L+1,L+1)       
     1 +2.0*( QEZ1(L+2,L+1)+QEZ1(L+1,L)-2.0*QEZ1(L+1,L+1) )   
     1   / EPS(L+1,L+1)   
   20 CONTINUE   
C       
C       
      DO 40 M=1, N1-2     
      DO 30  L=M+1, N1-1  
      QEZ2(L+1,M+1) = 2.0*QEZ1(L+1,M+1) - QEZ2(L+1,M+1)       
     1 +( QEZ1(L+2,M+1)+QEZ1(L,M+1)-4.0*QEZ1(L+1,M+1)         
     1   +QEZ1(L+1,M+2)+QEZ1(L+1,M) )/ EPS(L+1,M+1)  
   30 CONTINUE   
   40 CONTINUE   
C       
C       
      T = DT*(I-1)*2.0    
      QEZ2(L0+1,M0+1)=QEZ2(L0+1,M0+1) + Q1*DSIN(OMEGA*T)      
C       
C       
C       
C     L=0, M=0   
C       
      QEZ1(1,1) = 2.0*QEZ2(1,1) - QEZ1(1,1) 
     1 + 4.0*( QEZ2(2,1) - QEZ2(1,1) )/ EPS(1,1)     
C       
C     M=0        
C       
      DO 110 L=1, N1-1    
      QEZ1(L+1,1) = 2.0*QEZ2(L+1,1) - QEZ1(L+1,1)    
     1 +( QEZ2(L+2,1)+QEZ2(L,1)-4.0*QEZ2(L+1,1)      
     1   +2.0*QEZ2(L+1,2) )/ EPS(L+1,1)     
  110 CONTINUE   
C       
C     L=N1, M=0  
C       
      QEZ1(N1+1,1) = 2.0*QEZ2(N1+1,1) - QEZ1(N1+1,1) 
     1 +( 2.0*QEZ2(N1,1)-4.0*QEZ2(N1+1,1)   
     1   +2.0*QEZ2(N1+1,2) )/ EPS(N1+1,1)   
C       
C     L=N1       
C       
      DO 115 M=1, N1-1    
      QEZ1(N1+1,M+1) = 2.0*QEZ2(N1+1,M+1) - QEZ1(N1+1,M+1)    
     1 +( 2.0*QEZ2(N1,M+1)-4.0*QEZ2(N1+1,M+1)        
     1   +QEZ2(N1+1,M+2)+QEZ2(N1+1,M) )/ EPS(N1+1,M+1)        
  115 CONTINUE   
C       
C     L=N1, M=N1 
C       
      QEZ1(N1+1,N1+1) = 2.0*QEZ2(N1+1,N1+1) - QEZ1(N1+1,N1+1) 
     1 +4.0*( QEZ2(N1+1,N1) - QEZ2(N1+1,N1+1) )/ EPS(N1+1,N1+1)        
C       
C     L=M        
C       
      DO 120 L=1, N1-1    
      QEZ1(L+1,L+1) = 2.0*QEZ2(L+1,L+1) - QEZ1(L+1,L+1)       
     1 +2.0*( QEZ2(L+2,L+1)+QEZ2(L+1,L)-2.0*QEZ2(L+1,L+1) )   
     1   / EPS(L+1,L+1)   
  120 CONTINUE   
C       
C       
      DO 140 M=1, N1-2    
      DO 130  L=M+1, N1-1 
      QEZ1(L+1,M+1) = 2.0*QEZ2(L+1,M+1) - QEZ1(L+1,M+1)       
     1 +( QEZ2(L+2,M+1)+QEZ2(L,M+1)-4.0*QEZ2(L+1,M+1)         
     1   +QEZ2(L+1,M+2)+QEZ2(L+1,M) )/ EPS(L+1,M+1)  
  130 CONTINUE   
  140 CONTINUE   
C       
C       
      T = DT*(I*2-1)      
      QEZ1(L0+1,M0+1)=QEZ1(L0+1,M0+1) + Q1*DSIN(OMEGA*T)      
  300 CONTINUE   
C       
C       
      RETURN     
      END        
C       
C       
C       
C  *****************************************************************   
C       
C       
      SUBROUTINE CALU(U,NX,LX,DX,OMEGA)     
C       
      IMPLICIT REAL*8 (A-H,O-Z)    
C       
      COMMON /C1/ QEZ1(241,241), QEZ2(241,241)       
      COMMON /C4/ QEPS(241,241)    
      DIMENSION HX(241,241), HY(241,241)    
      DATA PI / 3.141592653589D0 / 
C       
      N1 = NX*LX 
C       
      DO 10 J=1, N1       
      HX(N1+1,J) = QEZ1(N1+1,J)-QEZ1(N1+1,J+1)       
      HY(N1+1,J) = QEZ1(N1,J)-QEZ1(N1+1,J)  
   10 CONTINUE   
C       
      DO 50 I=1, N1       
      HX(I,I) = QEZ1(I,I)-QEZ1(I+1,I)       
      HY(I,I) = QEZ1(I+1,I)-QEZ1(I,I)       
   50 CONTINUE   
C       
C     I=N1, J=N1 
C       
      HX(N1+1,N1+1) = QEZ1(N1+1,N1+1)-QEZ1(N1+1,N1)  
      HY(N1+1,N1+1) = QEZ1(N1+1,N1)-QEZ1(N1+1,N1+1)  
C       
C       
      DO 100 J=1, N1-1    
      DO 90 I=J+1, N1     
      HX(I,J) = QEZ1(I,J)-QEZ1(I,J+1)       
      HY(I,J) = QEZ1(I+1,J)-QEZ1(I,J)       
   90 CONTINUE   
  100 CONTINUE   
C       
C       
      S4 = QEPS(1,1)*QEZ1(1,1)**2  
     1        + 2.0*QEPS(N1+1,1)*QEZ1(N1+1,1)**2     
     2        + QEPS(N1+1,N1+1)*QEZ1(N1+1,N1+1)**2   
C       
      SX4 = HX(1,1)**2 + 2.0*HX(N1+1,1)**2 + HX(N1+1,N1+1)**2 
      SY4 = HY(1,1)**2 + 2.0*HY(N1+1,1)**2 + HY(N1+1,N1+1)**2 
C       
      S2 = 0.0D0 
      SX2 = 0.0D0         
      SY2 = 0.0D0         
C       
      DO 320 I=2, N1      
      S2 = S2 + QEPS(I,1)*QEZ1(I,1)**2      
     1        + QEPS(N1+1,I)*QEZ1(N1+1,I)**2         
      SX2 = SX2 + HX(I,1)**2 + HX(N1+1,I)**2         
      SY2 = SY2 + HY(I,1)**2 + HY(N1+1,I)**2         
  320 CONTINUE   
C       
      S2 = S2*2.0         
      SX2 = SX2*2.0       
      SY2 = SY2*2.0       
C       
C       
      S1 = 0.0D0 
      SX1 = 0.0D0         
      SY1 = 0.0D0         
C       
      DO 360 J=2, N1-1    
      DO 350 I=J+1, N1    
      S1 = S1 + QEPS(I,J)*QEZ1(I,J)**2      
      SX1 = SX1 + HX(I,J)**2       
      SY1 = SY1 + HY(I,J)**2       
  350 CONTINUE   
  360 CONTINUE   
C       
      S1 = S1*2.0         
      SX1 = SX1*2.0       
      SY1 = SY1*2.0       
C       
      DO 370 I=2, N1      
      S1 = S1 + QEPS(I,I)*QEZ1(I,I)**2      
      SX1 = SX1 + HX(I,I)**2       
      SY1 = SY1 + HY(I,I)**2       
  370 CONTINUE   
C       
C       
      U = (S1*4.0 + S2*2.0 + S4)/NX/NX/8.0/PI        
      UX = (SX1*4.0 + SX2*2.0 + SX4)/NX/NX/8.0/PI/OMEGA/OMEGA/DX/DX    
      UY = (SY1*4.0 + SY2*2.0 + SY4)/NX/NX/8.0/PI/OMEGA/OMEGA/DX/DX    
C       
      U = U + UX + UY     
C       
C       
      RETURN     
      END        
C       
C       
C  *****************************************************************   
C       
C       
      SUBROUTINE INIRAD(NX, LX)    
C       
      IMPLICIT REAL*8 (A-H,O-Z)    
C       
      COMMON /C1/ QEZ1(241,241), QEZ2(241,241)       
C       
C       
      DO 100 J=1, NX*LX+1 
      DO  90 I=J, NX*LX+1 
      QEZ1(I,J) = 0.0D0   
      QEZ2(I,J) = 0.0D0   
   90 CONTINUE   
  100 CONTINUE   
C       
C       
      RETURN     
      END        
C       
C       
C  *****************************************************************   
C       
C       
      SUBROUTINE CALEPS(A, R, QEPS1, QEPS2, QEPSD, NX, LX)    
C       
      IMPLICIT REAL*8 (A-H,O-Z)    
C       
      COMMON /C4/ QEPS(241,241)    
C       
C       
      R2 = R*R/A/A        
C       
C       
      DO 100 J=0, NX      
      Y = DFLOAT(J)/DFLOAT(NX)     
C       
      DO 90  I=0, NX      
      X = DFLOAT(I)/DFLOAT(NX)     
C       
      IF(X*X+Y*Y.LT.R2) GO TO 50   
      IF((X-1.0)*(X-1.0)+Y*Y.LT.R2) GO TO 50         
      IF(X*X+(Y-1.0)*(Y-1.0).LT.R2) GO TO 50         
      IF((X-1.0)*(X-1.0)+(Y-1.0)*(Y-1.0).LT.R2) GO TO 50      
      QEPS(I+1,J+1) = QEPS2        
      GO TO 90   
C       
   50 QEPS(I+1,J+1) = QEPS1        
C       
   90 CONTINUE   
  100 CONTINUE   
C       
C       
      DO 200 L1=0, LX-1   
      NL1 = L1*NX         
C       
      DO 190 L2=0, LX-1   
      NL2 = L2*NX         
C       
      DO 180 J=1, NX+1    
      DO 170 I=1, NX+1    
      QEPS(NL1+I,NL2+J) = QEPS(I,J)         
  170 CONTINUE   
  180 CONTINUE   
C       
  190 CONTINUE   
  200 CONTINUE   
C       
C       
      DO 300 J=0, NX      
      Y = DFLOAT(J)/NX    
C       
      DO 290 I=0, NX      
      X = DFLOAT(I)/NX    
C       
      IF(X*X+Y*Y.LT.R2) QEPS(I+1,J+1) = QEPSD        
  290 CONTINUE   
  300 CONTINUE   
C       
C       
      RETURN     
      END        

