Hi.(adsbygoogle = window.adsbygoogle || []).push({});

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 [Broken]

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

**Physics Forums | Science Articles, Homework Help, Discussion**

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Checking programming code: fortran 95

**Physics Forums | Science Articles, Homework Help, Discussion**