C
C
C     DEFTRI.FORT(A1) Programmed by K.Sakoda (1997)
C
C     Defect Mode in 2D Triangular Lattice of Circular Cylinders
C        Totally Symmetric (A1) Mode Calculated on 1/12 Lattice
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      COMMON /C1/ QEZ1(361,361), QEZ2(361,361)
      COMMON /C2/ N1, N2
      COMMON /C4/ QEPS(361,361)
      COMMON /C5/ S3
C
C
      DATA OM1,OM2,NPOINT / 0.4680D0, 0.4680D0,  1 / 
      DATA NX, LX, NT, NC / 15, 8, 720, 100 /
      DATA A, R, RD  / 1.27D0, 0.48D0, 0.00D0 /
      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     RD : Radius of Defect 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='ta1p.dat')
C     
      S3 = DSQRT(3.0D0) 
      NX = NX*3  
      N1 = NX*LX  
      N2 = N1/3*2 
C     
C     
C  10 READ(3,*) QEPSD   
C     IF(QEPSD.LT.0.0D0) GO TO 70   
C     
C     
      CALL REC1(A,R,RD,QEPS1,QEPS2,QEPSD,NX,LX,X0,Y0,NT,NC) 
      CALL CALEPS(A, R, RD, 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=1, NPOINT     
      OMEGA = OM1 + I1*(OM2-OM1)/NPOINT   
      DT = PI/OMEGA/DFLOAT(NT)      
      CALL INIRAD 
C     
      DO 400 I2=1, NC   
C     
      CALL CALRAD(OMEGA,NT,L0,M0,DX,DT,EMU)     
C     CALL CALU(U,NX,DX,OMEGA)      
C     CALL REC2(OMEGA, U)   
C    
  400 CONTINUE    
C     
      CALL CALU(U,NX,DX,OMEGA)     
      CALL REC2(OMEGA, U)   
      CALL REC4(NX,LX)  
C     
C     REC1 : Input Data      
C     REC2 : Omega vs Energy 
C     REC4 : 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 (Q)      
C      
      COMMON /C1/ QEZ1(361,361), QEZ2(361,361)   
      DIMENSION X1(1), Y1(1), E1(1)  
      DIMENSION X2(360), Y2(360), E2(360)  
      DIMENSION X3(360), Y3(360), E3(360)  
      DIMENSION X4(5000), Y4(5000), E4(5000)     
      DIMENSION X5(10000), Y5(10000), E5(10000)  
      DIMENSION X6(360), Y6(360), E6(360)  
C      
      DATA ND /  8 /     
C     ND: EVEN     
C      
      N1 = NX*LX   
      N2 = N1/3*2  
C      
C      
      X1(1) = 0.0  
      Y1(1) = 0.0  
      E1(1) = SNGL(QEZ1(1,1))  
      A1 = E1(1)   
C      
C     
      M2 = 0 
      DO 10 I=ND, N1, ND 
      M2 = M2+1    
      X2(M2) = FLOAT(I)/NX     
      Y2(M2) = 0.0 
      E2(M2) = SNGL(QEZ1(I+1,1))     
      IF(ABS(E2(M2)).GT.ABS(A1)) A1 = E2(M2)     
   10 CONTINUE     
C      
C      
      M3 = 0 
      DO 20 I=ND, N2, ND 
      M3 = M3+1    
      X3(M3) = FLOAT(I)/NX*1.5 
      Y3(M3) = FLOAT(I)/NX*0.86603   
      E3(M3) = SNGL(QEZ1(I+1,I+1))   
      IF(ABS(E3(M3)).GT.ABS(A1)) A1 = E3(M3)     
   20 CONTINUE     
C      
C     
      M4 = 0 
      DO 40 J=ND, N2-ND, ND    
      DO 30 I=J+ND, N1-J/2, ND 
      M4 = M4+1    
      X4(M4) = FLOAT(I)/NX + FLOAT(J)/NX/2 
      Y4(M4) = FLOAT(J)/NX*0.86603   
      E4(M4) = SNGL(QEZ1(I+1,J+1))   
      IF(ABS(E4(M4)).GT.ABS(A1)) A1 = E4(M4)     
   30 CONTINUE     
   40 CONTINUE     
C      
C      
      E1(1) = E1(1)/A1   
C      
C      
      DO 50 I=1, M2      
      E2(I) = E2(I)/A1   
   50 CONTINUE     
C      
C      
      DO 60 I=1, M3      
      E3(I) = E3(I)/A1   
   60 CONTINUE     
C      
C      
      DO 70 I=1, M4      
      E4(I) = E4(I)/A1   
   70 CONTINUE     
C      
C      
      DO 80 I=1, M4      
      X5(I) = (X4(I) + 1.73205*Y4(I))/2.0  
      Y5(I) = (1.73205*X4(I) - Y4(I))/2.0  
      E5(I) = E4(I)      
   80 CONTINUE     
C      
C      
      DO 90 I=1, M2      
      X5(M4+I) = (X2(I) + 1.73205*Y2(I))/2.0     
      Y5(M4+I) = (1.73205*X2(I) - Y2(I))/2.0     
      E5(M4+I) = E2(I)   
   90 CONTINUE     
C      
C      
      M5=M4+M2     
      DO 100 I=1, M4     
      X5(M5+I) = (X4(I) - 1.73205*Y4(I))/2.0     
      Y5(M5+I) = (1.73205*X4(I) + Y4(I))/2.0     
      E5(M5+I) = E4(I)   
  100 CONTINUE     
      M5=M5+M4     
C      
C      
      DO 110 I=1, M3     
      X6(I) = (X3(I) - 1.73205*Y3(I))/2.0  
      Y6(I) = (1.73205*X3(I) + Y3(I))/2.0  
      E6(I) = E3(I)      
  110 CONTINUE     
      M6=M3  
C      
C      
      WRITE(4,600) X1(1), Y1(1), E1(1)     
C      
C      
      DO 200 I=1, M2     
      WRITE(4,600) X2(I), Y2(I), E2(I)     
  200 CONTINUE     
C      
C      
      DO 210 I=1, M3     
      WRITE(4,600) X3(I), Y3(I), E3(I)     
  210 CONTINUE     
C      
C      
      DO 220 I=1, M4     
      WRITE(4,600) X4(I), Y4(I), E4(I)     
  220 CONTINUE     
C      
C      
      DO 230 I=1, M5     
      WRITE(4,600) X5(I), Y5(I), E5(I)     
  230 CONTINUE     
C      
C      
      DO 240 I=1, M6     
      WRITE(4,600) X6(I), Y6(I), E6(I)     
  240 CONTINUE     
C      
C      
      DO 300 I=1, M2     
      WRITE(4,600) -X2(I), Y2(I), E2(I)    
  300 CONTINUE     
C      
C      
      DO 310 I=1, M3     
      WRITE(4,600) -X3(I), Y3(I), E3(I)    
  310 CONTINUE     
C      
C      
      DO 320 I=1, M4     
      WRITE(4,600) -X4(I), Y4(I), E4(I)    
  320 CONTINUE     
C      
C      
      DO 330 I=1, M5     
      WRITE(4,600) -X5(I), Y5(I), E5(I)    
  330 CONTINUE     
C      
C      
      DO 410 I=1, M3     
      WRITE(4,600) X3(I), -Y3(I), E3(I)    
  410 CONTINUE     
C      
C      
      DO 420 I=1, M4     
      WRITE(4,600) X4(I), -Y4(I), E4(I)    
  420 CONTINUE     
C      
C      
      DO 430 I=1, M5     
      WRITE(4,600) X5(I), -Y5(I), E5(I)    
  430 CONTINUE     
C      
C      
      DO 440 I=1, M6     
      WRITE(4,600) X6(I), -Y6(I), E6(I)    
  440 CONTINUE     
C      
C      
      DO 510 I=1, M3     
      WRITE(4,600) -X3(I), -Y3(I), E3(I)   
  510 CONTINUE     
C      
C      
      DO 520 I=1, M4     
      WRITE(4,600) -X4(I), -Y4(I), E4(I)   
  520 CONTINUE     
C      
C      
      DO 530 I=1, M5     
      WRITE(4,600) -X5(I), -Y5(I), E5(I)   
  530 CONTINUE     
C      
C      
  600 FORMAT(1H ,F8.4,1H,,F8.4,1H,,F8.4)   
C      
C      
      RETURN 
      END    
C      
C      
C  *****************************************************************     
C      
C      
      SUBROUTINE REC1(A,R,RD,QEPS1,QEPS2,QEPSD,NX,LX,X0,Y0,NT,NC) 
C      
      IMPLICIT REAL*8 (A-H,O-Z)      
C      
      WRITE(4,603) A, R, RD    
  603 FORMAT(1H ,E15.7,5X,E15.7,5X,E15.7)  
      WRITE(4,600) QEPS1 
      WRITE(4,600) QEPS2 
      WRITE(4,600) QEPSD 
  600 FORMAT(1H ,E15.7,5X,E15.7)     
      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)      
C      
      IMPLICIT REAL*8 (A-H,O-Z)      
C      
      COMMON /C1/ QEZ1(361,361), QEZ2(361,361)   
      COMMON /C2/ N1, N2 
      COMMON /C4/ QEPS(361,361)      
      DIMENSION EPS(361,361)   
      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/3.0D0*2.0D0    
C      
C      
      DO 410 M=1, N2+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.0D0*QEZ1(1,1) - QEZ2(1,1)    
     1 +( QEZ1(2,1) + QEZ1(2,1) + QEZ1(2,1)      
     2  + QEZ1(2,1) + QEZ1(2,1) + QEZ1(2,1)      
     3  - 6.0D0*QEZ1(1,1) )    
     4  / EPS(1,1) 
C      
C     L=1, M=0     
C      
      QEZ2(2,1) = 2.0D0*QEZ1(2,1) - QEZ2(2,1)    
     1 +( QEZ1(3,1) + QEZ1(2,2) + QEZ1(2,1)      
     2  + QEZ1(1,1) + QEZ1(2,1) + QEZ1(2,2)      
     3  - 6.0D0*QEZ1(2,1) )    
     4  / EPS(2,1) 
C      
C     M=0    
C      
      DO 10 L=2, N1-1    
      QEZ2(L+1,1) = 2.0D0*QEZ1(L+1,1) - QEZ2(L+1,1)    
     1 +( QEZ1(L+2,1) + QEZ1(L+1,2) + QEZ1(L  ,2)      
     2  + QEZ1(L  ,1) + QEZ1(L,2) + QEZ1(L+1,2)  
     3  - 6.0D0*QEZ1(L+1,1) )  
     4  / EPS(L+1,1)     
   10 CONTINUE     
C      
C     L=N1, M=0    
C      
      QEZ2(N1+1,1) = 2.0D0*QEZ1(N1+1,1) - QEZ2(N1+1,1) 
     1 +( QEZ1(N1  ,1) + QEZ1(N1  ,2) + QEZ1(N1  ,2)   
     2  + QEZ1(N1  ,1) + QEZ1(N1  ,2) + QEZ1(N1  ,2)   
     3  - 6.0D0*QEZ1(N1+1,1) ) 
     4  / EPS(N1+1,1)    
C      
C     L=MAX  
C      
      DO 15 M=2, N2-2, 2 
      L = N1-M/2   
      QEZ2(L+1,M+1) = 2.0D0*QEZ1(L+1,M+1) - QEZ2(L+1,M+1)    
     1 +( QEZ1(L  ,M+1) + QEZ1(L  ,M+2) + QEZ1(L  ,M+2)      
     2  + QEZ1(L  ,M+1) + QEZ1(L+1,M  ) + QEZ1(L+1,M  )      
     3  - 6.0D0*QEZ1(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
   15 CONTINUE     
C      
C     L=MAX-1      
C      
      DO 17 M=1, N2-3, 2 
      L = N1-1-(M-1)/2   
      QEZ2(L+1,M+1) = 2.0D0*QEZ1(L+1,M+1) - QEZ2(L+1,M+1)    
     1 +( QEZ1(L+1,M+1) + QEZ1(L+1,M+2) + QEZ1(L  ,M+2)      
     2  + QEZ1(L  ,M+1) + QEZ1(L+1,M  ) + QEZ1(L+2,M  )      
     3  - 6.0D0*QEZ1(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
   17 CONTINUE     
C      
C     L=N2, M=N2-1 
C      
      QEZ2(N2+1,N2) = 2.0D0*QEZ1(N2+1,N2) - QEZ2(N2+1,N2)    
     1 +( QEZ1(N2+1,N2  ) + QEZ1(N2+1,N2+1) + QEZ1(N2+1,N2  )      
     2  + QEZ1(N2  ,N2  ) + QEZ1(N2+1,N2-1) + QEZ1(N2+2,N2-1)      
     3  - 6.0D0*QEZ1(N2+1,N2) )      
     4  / EPS(N2+1,N2)   
C      
C     L=N2, M=N2   
C      
      QEZ2(N2+1,N2+1) = 2.0D0*QEZ1(N2+1,N2+1) - QEZ2(N2+1,N2+1)    
     1 +( QEZ1(N2+1,N2  ) + QEZ1(N2+1,N2  ) + QEZ1(N2+1,N2  )      
     2  + QEZ1(N2+1,N2  ) + QEZ1(N2+1,N2  ) + QEZ1(N2+1,N2  )      
     3  - 6.0D0*QEZ1(N2+1,N2+1) )    
     4  / EPS(N2+1,N2+1) 
C      
C     L=MIN  
C      
      DO 19 M=1, N2-1    
      L = M  
      QEZ2(L+1,M+1) = 2.0D0*QEZ1(L+1,M+1) - QEZ2(L+1,M+1)    
     1 +( QEZ1(L+2,M+1) + QEZ1(L+2,M+1) + QEZ1(L+2,M  )      
     2  + QEZ1(L+1,M  ) + QEZ1(L+1,M  ) + QEZ1(L+2,M  )      
     3  - 6.0D0*QEZ1(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
   19 CONTINUE     
C      
C     L=MIN+1      
C      
      DO 20 M=1, N2-2    
      L = M+1      
      QEZ2(L+1,M+1) = 2.0D0*QEZ1(L+1,M+1) - QEZ2(L+1,M+1)    
     1 +( QEZ1(L+2,M+1) + QEZ1(L+1,M+2) + QEZ1(L+1,M+1)      
     2  + QEZ1(L  ,M+1) + QEZ1(L+1,M  ) + QEZ1(L+2,M  )      
     3  - 6.0D0*QEZ1(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
   20 CONTINUE     
C      
C      
      DO 31 M=1, N2-3, 2 
      M1 = (M-1)/2 
      DO 30 L=M+2, N1-2-M1     
      QEZ2(L+1,M+1) = 2.0D0*QEZ1(L+1,M+1) - QEZ2(L+1,M+1)    
     1 +( QEZ1(L+2,M+1) + QEZ1(L+1,M+2) + QEZ1(L  ,M+2)      
     2  + QEZ1(L  ,M+1) + QEZ1(L+1,M  ) + QEZ1(L+2,M  )      
     3  - 6.0D0*QEZ1(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
   30 CONTINUE     
   31 CONTINUE     
C      
C      
      DO 41 M=2, N2-2, 2 
      M1 = (M-2)/2 
      DO 40 L=M+2, N1-2-M1     
      QEZ2(L+1,M+1) = 2.0D0*QEZ1(L+1,M+1) - QEZ2(L+1,M+1)    
     1 +( QEZ1(L+2,M+1) + QEZ1(L+1,M+2) + QEZ1(L  ,M+2)      
     2  + QEZ1(L  ,M+1) + QEZ1(L+1,M  ) + QEZ1(L+2,M  )      
     3  - 6.0D0*QEZ1(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
   40 CONTINUE     
   41 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.0D0*QEZ2(1,1) - QEZ1(1,1)    
     1 +( QEZ2(2,1) + QEZ2(2,1) + QEZ2(2,1)      
     2  + QEZ2(2,1) + QEZ2(2,1) + QEZ2(2,1)      
     3  - 6.0D0*QEZ2(1,1) )    
     4  / EPS(1,1) 
C      
C     L=1, M=0     
C      
      QEZ1(2,1) = 2.0D0*QEZ2(2,1) - QEZ1(2,1)    
     1 +( QEZ2(3,1) + QEZ2(2,2) + QEZ2(2,1)      
     2  + QEZ2(1,1) + QEZ2(2,1) + QEZ2(2,2)      
     3  - 6.0D0*QEZ2(2,1) )    
     4  / EPS(2,1) 
C      
C     M=0    
C      
      DO 110 L=2, N1-1   
      QEZ1(L+1,1) = 2.0D0*QEZ2(L+1,1) - QEZ1(L+1,1)    
     1 +( QEZ2(L+2,1) + QEZ2(L+1,2) + QEZ2(L  ,2)      
     2  + QEZ2(L  ,1) + QEZ2(L,2) + QEZ2(L+1,2)  
     3  - 6.0D0*QEZ2(L+1,1) )  
     4  / EPS(L+1,1)     
  110 CONTINUE     
C      
C     L=N1, M=0    
C      
      QEZ1(N1+1,1) = 2.0D0*QEZ2(N1+1,1) - QEZ1(N1+1,1) 
     1 +( QEZ2(N1  ,1) + QEZ2(N1  ,2) + QEZ2(N1  ,2)   
     2  + QEZ2(N1  ,1) + QEZ2(N1  ,2) + QEZ2(N1  ,2)   
     3  - 6.0D0*QEZ2(N1+1,1) ) 
     4  / EPS(N1+1,1)    
C      
C     L=MAX  
C      
      DO 115 M=2, N2-2, 2      
      L = N1-M/2   
      QEZ1(L+1,M+1) = 2.0D0*QEZ2(L+1,M+1) - QEZ1(L+1,M+1)    
     1 +( QEZ2(L  ,M+1) + QEZ2(L  ,M+2) + QEZ2(L  ,M+2)      
     2  + QEZ2(L  ,M+1) + QEZ2(L+1,M  ) + QEZ2(L+1,M  )      
     3  - 6.0D0*QEZ2(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
  115 CONTINUE     
C      
C     L=MAX-1      
C      
      DO 117 M=1, N2-3, 2      
      L = N1-1-(M-1)/2   
      QEZ1(L+1,M+1) = 2.0D0*QEZ2(L+1,M+1) - QEZ1(L+1,M+1)    
     1 +( QEZ2(L+1,M+1) + QEZ2(L+1,M+2) + QEZ2(L  ,M+2)      
     2  + QEZ2(L  ,M+1) + QEZ2(L+1,M  ) + QEZ2(L+2,M  )      
     3  - 6.0D0*QEZ2(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
  117 CONTINUE     
C      
C     L=N2, M=N2-1 
C      
      QEZ1(N2+1,N2) = 2.0D0*QEZ2(N2+1,N2) - QEZ1(N2+1,N2)    
     1 +( QEZ2(N2+1,N2  ) + QEZ2(N2+1,N2+1) + QEZ2(N2+1,N2  )      
     2  + QEZ2(N2  ,N2  ) + QEZ2(N2+1,N2-1) + QEZ2(N2+2,N2-1)      
     3  - 6.0D0*QEZ2(N2+1,N2) )      
     4  / EPS(N2+1,N2)   
C      
C     L=N2, M=N2   
C      
      QEZ1(N2+1,N2+1) = 2.0D0*QEZ2(N2+1,N2+1) - QEZ1(N2+1,N2+1)    
     1 +( QEZ2(N2+1,N2  ) + QEZ2(N2+1,N2  ) + QEZ2(N2+1,N2  )      
     2  + QEZ2(N2+1,N2  ) + QEZ2(N2+1,N2  ) + QEZ2(N2+1,N2  )      
     3  - 6.0D0*QEZ2(N2+1,N2+1) )    
     4  / EPS(N2+1,N2+1) 
C      
C     L=MIN  
C      
      DO 119 M=1, N2-1   
      L = M  
      QEZ1(L+1,M+1) = 2.0D0*QEZ2(L+1,M+1) - QEZ1(L+1,M+1)    
     1 +( QEZ2(L+2,M+1) + QEZ2(L+2,M+1) + QEZ2(L+2,M  )      
     2  + QEZ2(L+1,M  ) + QEZ2(L+1,M  ) + QEZ2(L+2,M  )      
     3  - 6.0D0*QEZ2(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
  119 CONTINUE     
C      
C     L=MIN+1      
C      
      DO 120 M=1, N2-2   
      L = M+1      
      QEZ1(L+1,M+1) = 2.0D0*QEZ2(L+1,M+1) - QEZ1(L+1,M+1)    
     1 +( QEZ2(L+2,M+1) + QEZ2(L+1,M+2) + QEZ2(L+1,M+1)      
     2  + QEZ2(L  ,M+1) + QEZ2(L+1,M  ) + QEZ2(L+2,M  )      
     3  - 6.0D0*QEZ2(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
  120 CONTINUE     
C      
C      
      DO 131 M=1, N2-3, 2      
      M1 = (M-1)/2 
      DO 130 L=M+2, N1-2-M1    
      QEZ1(L+1,M+1) = 2.0D0*QEZ2(L+1,M+1) - QEZ1(L+1,M+1)    
     1 +( QEZ2(L+2,M+1) + QEZ2(L+1,M+2) + QEZ2(L  ,M+2)      
     2  + QEZ2(L  ,M+1) + QEZ2(L+1,M  ) + QEZ2(L+2,M  )      
     3  - 6.0D0*QEZ2(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
  130 CONTINUE     
  131 CONTINUE     
C      
C      
      DO 141 M=2, N2-2, 2      
      M1 = (M-2)/2 
      DO 140 L=M+2, N1-2-M1    
      QEZ1(L+1,M+1) = 2.0D0*QEZ2(L+1,M+1) - QEZ1(L+1,M+1)    
     1 +( QEZ2(L+2,M+1) + QEZ2(L+1,M+2) + QEZ2(L  ,M+2)      
     2  + QEZ2(L  ,M+1) + QEZ2(L+1,M  ) + QEZ2(L+2,M  )      
     3  - 6.0D0*QEZ2(L+1,M+1) )      
     4  / EPS(L+1,M+1)   
  140 CONTINUE     
  141 CONTINUE     
C      
C      
      T = DT*(I*2-1)     
      QEZ1(L0+1,M0+1)=QEZ1(L0+1,M0+1) + Q1*DSIN(OMEGA*T)     
C      
C      
C     CALL REC3    
C      
C      
  300 CONTINUE     
C      
C      
      RETURN 
      END    
C      
C      
C      
C  *****************************************************************     
C      
C      
      SUBROUTINE CALU(U,NX,DX,OMEGA) 
C      
      IMPLICIT REAL*8 (A-H,O-Z)      
C      
      COMMON /C1/ QEZ1(361,361), QEZ2(361,361)   
      COMMON /C2/ N1, N2 
      COMMON /C4/ QEPS(361,361)      
      COMMON /C5/ S3     
      DIMENSION HX(361,361), HY(361,361)   
      DATA PI / 3.141592653589D0 /   
C      
C     I=MAX  
C      
      DO 5 J=1, N2-1, 2  
      I = N1+1-(J-1)/2   
      HX(I,J) = QEZ1(I-1,J)+QEZ1(I,J)-2.0D0*QEZ1(I-1,J+1)    
      HY(I,J) = QEZ1(I-1,J)-QEZ1(I,J)      
    5 CONTINUE     
C      
C     I=MAX-1      
C      
      DO 10 J=2, N2, 2   
      I = N1-(J-2)/2     
      HX(I,J) = QEZ1(I,J)+QEZ1(I,J)-2.0D0*QEZ1(I,J+1)  
      HY(I,J) = QEZ1(I,J)-QEZ1(I,J)  
   10 CONTINUE     
C      
C     I=MIN  
C      
      DO 50 I=1, N2      
      HX(I,I) = QEZ1(I+1,I)+QEZ1(I,I)-2.0D0*QEZ1(I+1,I)      
      HY(I,I) = QEZ1(I+1,I)-QEZ1(I,I)      
   50 CONTINUE     
C      
C     I=N2+1, J=N2+1     
C      
      HX(N2+1,N2+1) = QEZ1(N2,N2)+QEZ1(N2+1,N2+1)      
     1   -2.0D0*QEZ1(N2,N2)    
      HY(N2+1,N2+1) = QEZ1(N2,N2)-QEZ1(N2+1,N2+1)      
C      
C      
      DO 91 J=1, N2-1, 2 
      DO 90 I=J+1, N1-(J-1)/2  
      HX(I,J) = QEZ1(I+1,J)+QEZ1(I,J)-2.0D0*QEZ1(I,J+1)      
      HY(I,J) = QEZ1(I+1,J)-QEZ1(I,J)      
   90 CONTINUE     
   91 CONTINUE     
C      
C      
      DO 101 J=2, N2-2, 2      
      DO 100 I=J+1, N1-J/2     
      HX(I,J) = QEZ1(I+1,J)+QEZ1(I,J)-2.0D0*QEZ1(I,J+1)      
      HY(I,J) = QEZ1(I+1,J)-QEZ1(I,J)      
  100 CONTINUE     
  101 CONTINUE     
C      
C      
      S12  = QEPS(1,1)*QEZ1(1,1)**2  
      SX12 = HX(1,1)**2  
      SY12 = HY(1,1)**2  
C      
      S6  = QEPS(N2+1,N2+1)*QEZ1(N2+1,N2+1)**2   
      SX6 = HX(N2+1,N2+1)**2   
      SY6 = HY(N2+1,N2+1)**2   
C      
      S4  = QEPS(N1+1,1)*QEZ1(N1+1,1)**2   
      SX4 = HX(N1+1,1)**2      
      SY4 = HY(N1+1,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     
      SX2 = SX2 + HX(I,1)**2   
      SY2 = SY2 + HY(I,1)**2   
  320 CONTINUE     
C      
      DO 340 J=3, N2-1, 2      
      I = N1-(J-3)/2     
      S2 = S2 + QEPS(I,J)*QEZ1(I,J)**2     
      SX2 = SX2 + HX(I,J)**2   
      SY2 = SY2 + HY(I,J)**2   
  340 CONTINUE     
C      
      DO 330 I=2, N2     
      S2 = S2 + QEPS(I,I)*QEZ1(I,I)**2     
      SX2 = SX2 + HX(I,I)**2   
      SY2 = SY2 + HY(I,I)**2   
  330 CONTINUE     
C      
C      
      S1 = 0.0D0   
      SX1 = 0.0D0  
      SY1 = 0.0D0  
C      
      DO 360 J=2, N2, 2  
      DO 350 I=J+1, N1-(J-2)/2 
      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      
      DO 380 J=3, N2-1, 2      
      DO 370 I=J+1, N1-(J-1)/2 
      S1 = S1 + QEPS(I,J)*QEZ1(I,J)**2     
      SX1 = SX1 + HX(I,J)**2   
      SY1 = SY1 + HY(I,J)**2   
  370 CONTINUE     
  380 CONTINUE     
C      
C      
      U = (S1*12.0 + S2*6.0 + S4*3.0 +S6*2.0 + S12)/NX/NX/8.0/PI   
      UX = (SX1*12.0 + SX2*6.0 + SX4*3.0 + SX6*2.0 + SX12)   
     1  /NX/NX/8.0/PI/OMEGA/OMEGA/DX/DX/3.0      
      UY = (SY1*12.0 + SY2*6.0 + SY4*3.0 + SY6*2.0 + SY12)   
     1  /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  
C      
      IMPLICIT REAL*8 (A-H,O-Z)      
C      
      COMMON /C1/ QEZ1(361,361), QEZ2(361,361)   
      COMMON /C2/ N1, N2 
C      
C      
      DO 100 J=1, N2+1   
      DO  90 I=J, N1+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, RD, QEPS1, QEPS2, QEPSD, NX, LX)     
C      
      IMPLICIT REAL*8 (A-H,O-Z)      
C      
      COMMON /C4/ QEPS(361,361)      
      COMMON /C5/ S3     
C      
C      
      R2 = R*R/A/A 
      RD2 = RD*RD/A/A    
      S32 = S3/2.0D0     
C      
C      
      DO 100 J=0, NX     
      Y = DFLOAT(J)/DFLOAT(NX)*S32   
C      
      DO 90  I=0, NX     
      X = DFLOAT(I)/DFLOAT(NX) + Y/S3      
C      
      IF(X*X+Y*Y.LT.R2) GO TO 50     
      IF((X-1.0D0)*(X-1.0D0)+Y*Y.LT.R2) GO TO 50 
      IF((X-0.5D0)*(X-0.5D0)+(Y-S32)*(Y-S32).LT.R2) GO TO 50 
      IF((X-1.5D0)*(X-1.5D0)+(Y-S32)*(Y-S32).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*S32     
C      
      DO 290 I=0, NX     
      X = DFLOAT(I)/NX + Y/S3  
C      
      IF(X*X+Y*Y.LT.RD2) GO TO 270   
      IF(X*X+Y*Y.LT.R2)  GO TO 280   
      GO TO 290    
C      
  270 QEPS(I+1,J+1) = QEPSD    
      GO TO 290    
C      
  280 QEPS(I+1,J+1) = QEPS2    
C      
  290 CONTINUE     
  300 CONTINUE     
C      
C      
C     DO 400 J=0, NX     
C     DO 390 I=0, NX     
C     WRITE(4,599) I, J, QEPS(I+1,J+1)     
C 390 CONTINUE     
C 400 CONTINUE     
C 599 FORMAT(1H ,I5,5X,I5,5X,E15.7)  
C      
C      
      RETURN 
      END    
