C
C
C    One-dimensional System (uni.fort) Programmed by K. Sakoda
C                     2002.2.25
C
C     Field Distribution of One-Dimensional Structure with a Defect
C        Transfer Matrix Method
C
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)
C
      DIMENSION Q(2, 2)
      COMMON /C1/ E(10000, 3)
C
      DATA N / 28 / 
      DATA AL1, AL2, AL3 / 0.5D0, 0.5D0, 1.2D0 /  
      DATA QEPS1, QEPS2 / (2.1025D0, 0.0D0), (5.29D0, 0.0D0) /
      DATA QEPS3 / (2.1025D0, 0.0D0) /
      DATA OMEGA / 0.2651100D0 /
      DATA NX, NA / 50, 10 /
C
C     N: Number of regular layers
C     AL1: Thickness of regular layer 1 in unit of lattice constant
C     AL2: Thickness of regular layer 2 in unit of lattice constant
C     AL3: Thickness of defect layer in unit of lattice constant
C     QEPS1: Dielectric constant of layer 1
C     QEPS2: Dielectric constant of layer 2
C     QEPS3: Dielectric constant of defect layer
C     OMEGA: Angular frequency of incident plane wave
C     NX: Number of data points in a unit cell
C     NA: Number of lattice at entrance and exit
C      
      DATA PI / 3.141592653589793D0 /  
C
      OPEN(4,FILE='uni053.dat')
C
      CALL REC1(N, AL1, AL2, AL3, QEPS1, QEPS2, QEPS3, OMEGA, NX, NA)
C
C
      QETA1 = CDSQRT(QEPS1)
      QETA2 = CDSQRT(QEPS2)
      QETA3 = CDSQRT(QEPS3)
C
C     QETA1: Refractive index of layer1
C     QETA2: Refractive index of layer2
C     QETA3: Refractive index of defect layer
C
C
      QK0 = 2.0 * PI * OMEGA
      QK1 = QK0 * QETA1
      QK2 = QK0 * QETA2
      QK3 = QK0 * QETA3
C
      CALL CALMAT(N, AL1, AL2, AL3, QK1, QK2, QK3, QK0, Q) 
      CALL TANDR(QT, QR, Q)
C     CALL REC2(OMEGA, QT, QR)
C
      CALL CALE(N, AL1, AL2, AL3, QK1, QK2, QK3, QK0, NX, NA, QR) 
C
C
      STOP
      END
C
C
C  ***********************************************************
C
C
      SUBROUTINE REC1
     1       (N, AL1, AL2, AL3, QEPS1, QEPS2, QEPS3, OMEGA, NX, NA)
C
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)
C
      WRITE(4,601) N, AL1, AL2, AL3
  601 FORMAT(1H ,I5,5H,    ,E15.7,5H,    ,E15.7,5H,    ,E15.7)
C
      WRITE(4,602) QEPS1, QEPS2, QEPS3
  602 FORMAT(1H ,E15.7,5(5H,    ,E15.7))
C
      WRITE(4,603) OMEGA, NX, NA
  603 FORMAT(1H ,E15.7,5H,    ,I5,5H,    ,I5) 
C
C
      RETURN
      END
C
C
C  ********************************************************
C
C
      SUBROUTINE REC2(OMEGA, QT, QR)
C
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q) 
C 
C     T = DABS(QT)**2
C     R = DABS(QR)**2
C
      WRITE(4,600) OMEGA, QT, QR
  600 FORMAT(1H ,E15.7,4(3H,  ,E15.7))
C
  400 RETURN 
      END
C
C
C  *************************************************************
C
C
      SUBROUTINE CALMAT(N, AL1, AL2, AL3, QK1, QK2, QK3, QK0, QQ) 
C
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)
C
      DIMENSION Q(2, 2), QQ(2, 2), QM1(2, 2), QT1INV(2, 2)
      DIMENSION QAD(2, 2), QBD(2, 2), QCD(2, 2), QTN(2, 2), QMNP(2, 2)
      DIMENSION QDD(2, 2), QED(2, 2)
      DATA QI / (0.0D0, 1.0D0) /
      DATA Q0 / (0.0D0, 0.0D0) /
C
C
      QM1(1, 1) = (QK1+QK0)/2.0/QK1
      QM1(1, 2) = (QK1-QK0)/2.0/QK1
      QM1(2, 1) = (QK1-QK0)/2.0/QK1
      QM1(2, 2) = (QK1+QK0)/2.0/QK1
C
C     QM1: Matrix for region 0 to 1
C
C
      QT1INV(1, 1) = CDEXP( QI * QK1 * AL1 )
      QT1INV(1, 2) = Q0
      QT1INV(2, 1) = Q0
      QT1INV(2, 2) = CDEXP( - QI * QK1 * AL1 )
C
      Q1 = (QK1+QK2)/2.0/QK1 
      Q2 = (QK1-QK2)/2.0/QK1 
      Q3 = CDEXP( QI * QK1 * AL1 )
      Q4 = CDEXP( - QI * QK1 * AL1 )
      QAD(1, 1) = Q1 * Q3
      QAD(1, 2) = Q2 * Q3
      QAD(2, 1) = Q2 * Q4
      QAD(2, 2) = Q1 * Q4
C     QAD: Matrix from layer 2 to layer1
C
C
      Q5 = (QK2+QK1)/2.0/QK2 
      Q6 = (QK2-QK1)/2.0/QK2 
      Q7 = CDEXP( QI * QK2 * AL2 )
      Q8 = CDEXP( - QI * QK2 * AL2 )
      QBD(1, 1) = Q5 * Q7
      QBD(1, 2) = Q6 * Q7
      QBD(2, 1) = Q6 * Q8
      QBD(2, 2) = Q5 * Q8
C     QBD: Matrix for layer1 to layer2
C
      CALL MATPRO(QCD, QBD, QAD)
C     QCD: Matrix for layer2 via layer1 to layer2
C
C
      ALN = N/2.0 + AL3
      QTN(1, 1) = CDEXP( - QI * QK1 * ALN )
      QTN(1, 2) = Q0
      QTN(2, 1) = Q0
      QTN(2, 2) = CDEXP( QI * QK1 * ALN )
C
      Q1 = (QK0+QK1)/2.0/QK0
      Q2 = (QK0-QK1)/2.0/QK0
      Q3 = QI * (QK0-QK1) * ALN
      Q4 = QI * (QK0+QK1) * ALN
      QMNP(1, 1) = Q1 * CDEXP(-Q3)
      QMNP(1, 2) = Q2 * CDEXP(-Q4)
      QMNP(2, 1) = Q2 * CDEXP(Q4)
      QMNP(2, 2) = Q1 * CDEXP(Q3)
C     QMNP: Matrix for region n to outside
C
C
      CALL MATPRO(QQ, QT1INV, QM1)
      CALL MATPRO(Q, QBD, QQ)
C
      NN = (N-4)/8
      DO 200 I = 1, NN
      CALL MATPRO(QQ, QCD, Q)
      CALL MATPRO(Q, QCD, QQ)
  200 CONTINUE
C
      Q1 = (QK3+QK2)/2.0/QK3 
      Q2 = (QK3-QK2)/2.0/QK3 
      Q3 = CDEXP( QI * QK3 * AL3 )
      Q4 = CDEXP( - QI * QK3 * AL3 )
      QDD(1, 1) = Q1 * Q3
      QDD(1, 2) = Q2 * Q3
      QDD(2, 1) = Q2 * Q4
      QDD(2, 2) = Q1 * Q4
C
      CALL MATPRO(QQ, QDD, Q)
C
      Q5 = (QK2+QK3)/2.0/QK2 
      Q6 = (QK2-QK3)/2.0/QK2 
      Q7 = CDEXP( QI * QK2 * AL2 )
      Q8 = CDEXP( - QI * QK2 * AL2 )
      QED(1, 1) = Q5 * Q7
      QED(1, 2) = Q6 * Q7
      QED(2, 1) = Q6 * Q8
      QED(2, 2) = Q5 * Q8
C
      CALL MATPRO(Q, QED, QQ)
C
      NN = (N-4)/8
      DO 300 I = 1, NN
      CALL MATPRO(QQ, QCD, Q)
      CALL MATPRO(Q, QCD, QQ)
  300 CONTINUE
C
C
      CALL MATPRO(QQ, QAD, Q)
      CALL MATPRO(Q, QTN, QQ)
      CALL MATPRO(QQ, QMNP, Q)
C
C
      RETURN
      END
C
C
C  *************************************************************
C
C
      SUBROUTINE MATPRO(QA, QB, QC)
C
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)
C
      DIMENSION QA(2,2), QB(2,2), QC(2,2)
C
      DO 20 I=1, 2
      DO 10 J=1, 2
      QA(I, J) = QB(I, 1)*QC(1, J) + QB(I, 2)* QC(2, J)
   10 CONTINUE
   20 CONTINUE
C
      RETURN
      END
C
C
C  ***********************************************************
C
C
      SUBROUTINE TANDR(QT, QR, Q)
C
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)
C
      DIMENSION Q(2,2)
C
      QR = - Q(2, 1)/Q(2,2)
      QT = Q(1, 1) + Q(1, 2) * QR
C
      RETURN
      END
C
C
C  ***********************************************************
C
C
      SUBROUTINE CALE(N, AL1, AL2, AL3, QK1, QK2, QK3, QK0, NX, NA, QR)
C
      IMPLICIT REAL*8 (A-H,O-P,R-Z), COMPLEX*16 (Q)
C
      DIMENSION AL(100), QK(100)
      COMMON /C1/ E(10000,3)
C
      DATA QI / (0.0D0, 1.0D0) / 
C
C
      J = 0
      M1 = N/4
      M2 = NX*NA
C
C
      QA = 1.0D0
      QB = QR
C
      DO 10 I = M2, 0, -1
      J = J + 1
      X = - DFLOAT(I) / NX
      QX = QI * QK0 * X
      QE = QA*CDEXP(QX) + QB * CDEXP(-QX)
      E(J, 1) = X
      E(J, 2) = DREAL(QE)
      E(J, 3) = DIMAG(QE)
   10 CONTINUE
C
      QA0 = QA
      QB0 = QB
C
C
      ALX = 0.0D0
      AL(1) = ALX
      QK(1) = QK0
C
      DO 30 I = 1, M1
      I2 = 2*I
      ALX = ALX + AL1
      AL(I2) = ALX
      QK(I2) = QK1
      ALX = ALX + AL2
      AL(I2+1) = ALX
      QK(I2+1) = QK2
   30 CONTINUE
C
      I2 = N/2 + 2
      ALX = ALX + AL3
      AL(I2) = ALX
      QK(I2) = QK3
C
      DO 50 I = 1, M1
      I2 = 2*I + N/2 +1
      ALX = ALX + AL2
      AL(I2) = ALX
      QK(I2) = QK2
      ALX = ALX + AL1
      AL(I2+1) = ALX
      QK(I2+1) = QK1
   50 CONTINUE
C
      QK(N+3) = QK0
C
C
      DO 200 I = 1, N+1
C
      Q1 = ( QK(I+1) + QK(I) )/2/QK(I+1)
      Q2 = ( QK(I+1) - QK(I) )/2/QK(I+1)
      Q3 = CDEXP( QI*( QK(I+1) + QK(I) ) * AL(I) )
      Q4 = CDEXP( QI*( QK(I+1) - QK(I) ) * AL(I) )
C
      Q11 = Q1 /Q4
      Q12 = Q2/ Q3
      Q21 = Q2*Q3
      Q22 = Q1*Q4
C
      QA = Q11*QA0 + Q12*QB0
      QB = Q21*QA0 + Q22* QB0
C
      NX1 = NX * (AL(I+1)-AL(I)) / (AL1 + AL2)
C
      DO 180 I1 = 0, NX1
      J = J + 1
      X = AL(I) + (AL(I+1)-AL(I)) * I1 / NX1
      QX = QI * QK(I+1) * X
      QE = QA*CDEXP(QX) + QB * CDEXP(-QX)
      E(J, 1) = X
      E(J, 2) = DREAL(QE)
      E(J, 3) = DIMAG(QE)
  180 CONTINUE
C
      QA0 = QA
      QB0 = QB
C
  200 CONTINUE
C
C
      Q1 = ( QK0 + QK1 )/2/QK0
      Q2 = ( QK0 - QK1 )/2/QK0
      Q3 = CDEXP( QI*( QK0 + QK1 ) * AL(N+2) )
      Q4 = CDEXP( QI*( QK0 - QK1 ) * AL(N+2) )
C
      Q11 = Q1 /Q4
      Q12 = Q2/ Q3
      Q21 = Q2*Q3
      Q22 = Q1*Q4
C
      QA = Q11*QA0 + Q12*QB0
C     QB = Q21*QA0 + Q22* QB0
C
      DO 210 I = 0, M2
      J = J + 1
      X = AL(N+2) + DFLOAT(I) / NX
      QX = QI * QK0 * X
C     QE = QA*CDEXP(QX) + QB * CDEXP(-QX)
      QE = QA*CDEXP(QX)
      E(J, 1) = X
      E(J, 2) = DREAL(QE)
      E(J, 3) = DIMAG(QE)
  210 CONTINUE
C
C
      AM = AL(N/2+1) + AL3/2.0
C
      DO 300 I = 1, J
      E(I, 1) = E(I, 1) - AM
  300 CONTINUE
C
C
      WRITE(4, 600) J
  600 FORMAT(1H ,I5)
C
      WRITE(4, 601) (E(I,1), E(I, 2), E(I,3), I=1, J)
  601 FORMAT(1H ,E15.7,5H,    ,E15.7,5H,    ,E15.7)
C
C
      RETURN
      END

