# Checking programming code: fortran 95

1. Jan 18, 2009

### 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 [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 Last edited by a moderator: May 3, 2017 2. Jan 19, 2009 ### jtbell ### Staff: Mentor 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? 3. Jan 19, 2009 ### 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.