- #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: