- #1
p$ycho
- 2
- 0
Hi.
Can anyone be kind enough to check what's wrong with my source code?
Here's the project: http://www.dur.ac.uk/3h.physics/proj...nsmission.html
You can skip the introduction bits.
I'm stuck at the last question on Calculation Part I. I did the source code independently and it works for single barrier and double barrier. But when I change it to N barriers, it didn't work, not even for single barrier.
Here's my source code:
Program Project
!
! Transmission of Electronics through a Semiconductor Barrier Structure for N barriers.
!
implicit none
double precision:: V, L, r1, t1, s1
double complex,dimension(2,2):: M, F, A_i, i_A, B_i
double complex:: k_A, k_B, m_A, m_B, zi
integer:: i, N, j, a, b, p, x
double precision:: E, delta_E, L_j
!
! Define the values of m_A, m_B, k_A, k_B, V, L_j
! [L]=Angstroms=10**-10m
! 1 Bohr radius=0.53 Angstroms
!
m_A=(0.067d0,0.0d0)
m_B=(0.096d0,0.0d0)
L_j=33.9d0/0.53d0
V=0.245d0
zi=(0.0,1.0)
write(*,*)'Please enter the number of barriers, x, where x is an integer more than zero.'
read(*,*) x
N=2*x ! N is the number of interface.
open(10,file='graph_nb')
!
! Do loop. To varies the values of E from 0 to 1.0eV.
!
E=0.005d0
delta_E=0.005d0
p=199
do i=1,p
k_A=dcmplx((m_A*E/13.6057d0)**(0.5d0),0.0d0)
if (E.ge.V) then
k_B=dcmplx((m_B*(E-V)/13.6057d0)**(0.5d0),0.0d0)
else
k_B=dcmplx(0.0d0,(m_B*(V-E)/13.6057d0)**(0.5d0))
end if
!
! Computing for N interfaces.
! Define identity matrix M.
!
M(1,1)=1
M(1,2)=0
M(2,1)=0
M(2,2)=1 ! M should be identity matrix.
do b=N,1,-1
!
! Test if b is even or odd.
!
if (((dble(b/2)-0.001d0).lt.(dble(b)/2.0d0)) .and. ((dble(b)/2.0d0).lt.(dble(b/2)+0.001d0))) then
print*, 'b is even', b
call A_matrix(k_A,m_A,dble(b-1)*L_j,A_i)
call A_matrix(k_B,m_B,dble(b-1)*L_j,B_i)
call invert2x2(A_i,i_A)
call multiply2x2(i_A,B_i,F)
else
print*, 'b is odd',b
call A_matrix(k_B,m_B,dble(b-1)*L_j,A_i)
call A_matrix(k_A,m_A,dble(b-1)*L_j,B_i)
call invert2x2(A_i,i_A)
call multiply2x2(i_A,B_i,F)
end if
call multiply2x2(M,F,M)
end do
!
! Calculating Reflection coefficient,r1 and Transmission coefficient,t1.
! Calculate the resultant reflection and transmission coefficients
! Let s1 be the sum of r1 and t1, s1=r1+t1
!
r1=(abs(-M(2,1)/M(2,2)))**2
t1=(abs(M(1,1)-M(1,2)*M(2,1)/M(2,2)))**2
s1=r1+t1
write(10,10) E, r1, t1, s1
10 format(4F9.6)
E=E+delta_E
end do
end program project ! end program project
subroutine invert2x2(A,B)
!
! On return to the calling program, the matrix B contains the inverse of 2x2 matrix A
! C is the determinant of A.
!
implicit none
double complex,dimension(2,2),intent(in):: A
double complex,dimension(2,2),intent(out):: B
double complex:: C
C=A(1,1)*A(2,2)-A(1,2)*A(2,1)
if (C.ne.(0d0,0d0)) then
B(1,1)=A(2,2)/C
B(1,2)=-A(1,2)/C
B(2,1)=-A(2,1)/C
B(2,2)=A(1,1)/C
else
write(*,*) 'Subroutine invert2x2 has failed.'
STOP
end if
return
end subroutine invert2x2
subroutine multiply2x2(A,B,C)
!
! On return to the calling program, the matrix C contains the multiplication of A and B
! C=B*A
!
implicit none
double complex, dimension(2,2), intent(in):: A, B
double complex, dimension(2,2), intent(out):: C
C(1,1)=A(1,1)*B(1,1)+A(1,2)*B(2,1)
C(1,2)=A(1,1)*B(1,2)+A(1,2)*B(2,2)
C(2,1)=A(2,1)*B(1,1)+A(2,2)*B(2,1)
C(2,2)=A(2,1)*B(1,2)+A(2,2)*B(2,2)
return
end subroutine multiply2x2
Subroutine A_matrix(k_A,m_A,L_j,A_i)
implicit none
double precision:: L_j
double complex:: zi, k_A, k_B, m_A, m_B
double complex, dimension(2,2):: A_i
zi=(0.0,1.0)
print*, k_A, m_A, L_j
A_i(1,1)=exp(zi*k_A*L_j)
A_i(1,2)=exp(-zi*k_A*L_j)
A_i(2,1)=k_A*exp(zi*k_A*L_j)/m_A
A_i(2,2)=-k_A*exp(-zi*k_A*L_j)/m_A
print*,A_i(1,1),A_i(1,2),A_i(2,1),A_i(2,2)
return
end subroutine A_matrix
Can anyone be kind enough to check what's wrong with my source code?
Here's the project: http://www.dur.ac.uk/3h.physics/proj...nsmission.html
You can skip the introduction bits.
I'm stuck at the last question on Calculation Part I. I did the source code independently and it works for single barrier and double barrier. But when I change it to N barriers, it didn't work, not even for single barrier.
Here's my source code:
Program Project
!
! Transmission of Electronics through a Semiconductor Barrier Structure for N barriers.
!
implicit none
double precision:: V, L, r1, t1, s1
double complex,dimension(2,2):: M, F, A_i, i_A, B_i
double complex:: k_A, k_B, m_A, m_B, zi
integer:: i, N, j, a, b, p, x
double precision:: E, delta_E, L_j
!
! Define the values of m_A, m_B, k_A, k_B, V, L_j
! [L]=Angstroms=10**-10m
! 1 Bohr radius=0.53 Angstroms
!
m_A=(0.067d0,0.0d0)
m_B=(0.096d0,0.0d0)
L_j=33.9d0/0.53d0
V=0.245d0
zi=(0.0,1.0)
write(*,*)'Please enter the number of barriers, x, where x is an integer more than zero.'
read(*,*) x
N=2*x ! N is the number of interface.
open(10,file='graph_nb')
!
! Do loop. To varies the values of E from 0 to 1.0eV.
!
E=0.005d0
delta_E=0.005d0
p=199
do i=1,p
k_A=dcmplx((m_A*E/13.6057d0)**(0.5d0),0.0d0)
if (E.ge.V) then
k_B=dcmplx((m_B*(E-V)/13.6057d0)**(0.5d0),0.0d0)
else
k_B=dcmplx(0.0d0,(m_B*(V-E)/13.6057d0)**(0.5d0))
end if
!
! Computing for N interfaces.
! Define identity matrix M.
!
M(1,1)=1
M(1,2)=0
M(2,1)=0
M(2,2)=1 ! M should be identity matrix.
do b=N,1,-1
!
! Test if b is even or odd.
!
if (((dble(b/2)-0.001d0).lt.(dble(b)/2.0d0)) .and. ((dble(b)/2.0d0).lt.(dble(b/2)+0.001d0))) then
print*, 'b is even', b
call A_matrix(k_A,m_A,dble(b-1)*L_j,A_i)
call A_matrix(k_B,m_B,dble(b-1)*L_j,B_i)
call invert2x2(A_i,i_A)
call multiply2x2(i_A,B_i,F)
else
print*, 'b is odd',b
call A_matrix(k_B,m_B,dble(b-1)*L_j,A_i)
call A_matrix(k_A,m_A,dble(b-1)*L_j,B_i)
call invert2x2(A_i,i_A)
call multiply2x2(i_A,B_i,F)
end if
call multiply2x2(M,F,M)
end do
!
! Calculating Reflection coefficient,r1 and Transmission coefficient,t1.
! Calculate the resultant reflection and transmission coefficients
! Let s1 be the sum of r1 and t1, s1=r1+t1
!
r1=(abs(-M(2,1)/M(2,2)))**2
t1=(abs(M(1,1)-M(1,2)*M(2,1)/M(2,2)))**2
s1=r1+t1
write(10,10) E, r1, t1, s1
10 format(4F9.6)
E=E+delta_E
end do
end program project ! end program project
subroutine invert2x2(A,B)
!
! On return to the calling program, the matrix B contains the inverse of 2x2 matrix A
! C is the determinant of A.
!
implicit none
double complex,dimension(2,2),intent(in):: A
double complex,dimension(2,2),intent(out):: B
double complex:: C
C=A(1,1)*A(2,2)-A(1,2)*A(2,1)
if (C.ne.(0d0,0d0)) then
B(1,1)=A(2,2)/C
B(1,2)=-A(1,2)/C
B(2,1)=-A(2,1)/C
B(2,2)=A(1,1)/C
else
write(*,*) 'Subroutine invert2x2 has failed.'
STOP
end if
return
end subroutine invert2x2
subroutine multiply2x2(A,B,C)
!
! On return to the calling program, the matrix C contains the multiplication of A and B
! C=B*A
!
implicit none
double complex, dimension(2,2), intent(in):: A, B
double complex, dimension(2,2), intent(out):: C
C(1,1)=A(1,1)*B(1,1)+A(1,2)*B(2,1)
C(1,2)=A(1,1)*B(1,2)+A(1,2)*B(2,2)
C(2,1)=A(2,1)*B(1,1)+A(2,2)*B(2,1)
C(2,2)=A(2,1)*B(1,2)+A(2,2)*B(2,2)
return
end subroutine multiply2x2
Subroutine A_matrix(k_A,m_A,L_j,A_i)
implicit none
double precision:: L_j
double complex:: zi, k_A, k_B, m_A, m_B
double complex, dimension(2,2):: A_i
zi=(0.0,1.0)
print*, k_A, m_A, L_j
A_i(1,1)=exp(zi*k_A*L_j)
A_i(1,2)=exp(-zi*k_A*L_j)
A_i(2,1)=k_A*exp(zi*k_A*L_j)/m_A
A_i(2,2)=-k_A*exp(-zi*k_A*L_j)/m_A
print*,A_i(1,1),A_i(1,2),A_i(2,1),A_i(2,2)
return
end subroutine A_matrix
Last edited by a moderator: