C
C
C    One-dimensional System (uni.fort) Programmed by K. Sakoda
C                     2002.2.25
C
C     Transmittance 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)
C
      DATA N /  28 / 
      DATA AL1, AL2, AL3 / 0.5D0, 0.5D0, 0.5D0 /  
      DATA QEPS1, QEPS2 / (2.1025D0, 0.0D0), (5.29D0, 0.0D0) /
      DATA QEPS3 / (2.1025D0, 0.0D0) /
      DATA OMEGA1, OMEGA2, NOMEGA / 0.05D0, 1.0D0, 9500 /
C      
      DATA PI / 3.141592653589793D0 /  
C
      OPEN(4,FILE='uni004.dat')
C
      CALL REC1(N, AL1, AL2, AL3, QEPS1, QEPS2, QEPS3, 
     1                           OMEGA1, OMEGA2, NOMEGA)
C
C
      QETA1 = CDSQRT(QEPS1)
      QETA2 = CDSQRT(QEPS2)
      QETA3 = CDSQRT(QEPS3)
C
C
      DO 400 I = 0, NOMEGA
      OMEGA = OMEGA1 + I * (OMEGA2 - OMEGA1) / NOMEGA
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)
      CALL REC2(OMEGA, QT, QR)
C
  400 CONTINUE
C
C
      STOP
      END
C
C
C  ***********************************************************
C
C
      SUBROUTINE REC1(N, AL1, AL2, AL3, QEPS1, QEPS2, QEPS3,
     1                                  OMEGA1,OMEGA2,NOMEGA)
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) OMEGA1, OMEGA2, NOMEGA
  603 FORMAT(1H ,E15.7,5H,    ,E15.7,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
      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
      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
C
      CALL MATPRO(QCD, QBD, QAD)
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
      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

