Checking programming code: fortran 95

[p$ycho]

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

 

[jtbell]

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.

Please define "didn't work."

Did the program refuse to compile?

Did it compile, but crashed when you ran it? What error messages, if any?

Did it run, but give incorrect results? What input data did you give it, what results did you expect, and what did you actually get?

 

[p$ycho]

Thanks for the reply. I've found the mistake now. Instead of using my own subroutine for matrix multiplication, matmul(a,b) is much more convenient way and it gives the answer. For God sake, my demonstrators didn't even tell me that!

Previously, it did run. The programme didn't produce the expected answer.
Thanks again for taking the trouble to read through it, although I still don't understand what's wrong with my subroutine multiply2x2.