Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Checking programming code: fortran 95

  1. Jan 18, 2009 #1
    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
     
  2. jcsd
  3. Jan 19, 2009 #2

    jtbell

    User Avatar

    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?
     
  4. Jan 19, 2009 #3
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Checking programming code: fortran 95
  1. Fortran 95 (Replies: 12)

Loading...