code3d.f90

!!!!!!!!!!!!!!!!!!!!!!!!!!!   Program 3.D   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                       !
! Please Note:                                                          !
!                                                                       !
! (1) This computer program is written by Tao Pang in conjunction with  !
!     his book, "An Introduction to Computational Physics," published   !
!     by Cambridge University Press in 1997.                            !
!                                                                       !
! (2) No warranties, express or implied, are made for this program.     !
!                                                                       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
PROGRAM SCHR
!
! Main program for solving the eigenvalue problem of the
! one-dimensional Schroedinger equation. 
! Copyright (c) Tao Pang 1997.
!
  IMPLICIT NONE
  INTEGER, PARAMETER :: N=501,M=5,IMAX=100
  INTEGER :: I,IM,NL,NR,ISTEP
  REAL :: DL,H2M,EA,EB,E,DE,XL0,XR0,H,C
  REAL :: XL,XR,FACT,F0,F1,E1,SUM,V,VX
  REAL, DIMENSION (N) :: UL,UR,QL,QR,S
!
  DL    =  1.0E-06
  H2M   =  0.5 
  EA    =  2.4 
  EB    =  2.7
  E     =  EA
  DE    =  0.1
  XL0   = -10.0
  XR0   =  10.0
  H     =  (XR0-XL0)/(N-1)
  C     =  1.0/H2M
  UL(1) =  0.0
  UL(2) =  0.01
  UR(1) =  0.0
  UR(2) =  0.01
!
! Set up the potential q(x) and source s(x)
!
  DO I = 1, N
    XL    = XL0+(I-1)*H
    XR    = XR0-(I-1)*H
    QL(I) = C*(E-VX(XL))
    QR(I) = C*(E-VX(XR))
    S(I)  = 0.0
  END DO
!
! Find the matching point at the right turning point
!
  DO I = 1, N-1
    IF(((QL(I)*QL(I+1)).LE.0).AND.(QL(I).GT.0)) IM = I
  END DO
!    
! Numerov algorithm from left to right and vice versa
!
  NL = IM+1
  NR = N-IM+2
  CALL NMRV2 (NL,H,QL,S,UL)
  CALL NMRV2 (NR,H,QR,S,UR)
!
! Rescale the left solution
!
  FACT = UR(NR-1)/UL(IM)
  DO I = 1, NL
    UL(I) = FACT*UL(I)
  END DO
!
  F0 = UR(NR)+UL(NL)-UR(NR-2)-UL(NL-2)
  F0 = F0/(2.0*H*UR(NR-1))
!
! Bisection method for the root
!
  ISTEP = 0
  DO WHILE ((ABS(DE).GT.DL).AND.(ISTEP.LT.IMAX))
    E1 = E
    E  = (EA+EB)/2.0
    DO I = 1, N
      QL(I) = QL(I)+C*(E-E1)
      QR(I) = QR(I)+C*(E-E1)
    END DO
!
! Find the matching point at the right turning point
!
    DO I = 1, N-1
      IF(((QL(I)*QL(I+1)).LE.0).AND.(QL(I).GT.0)) IM = I
    END DO
!    
! Numerov algorithm from left to right and vice versa
!
    NL = IM+1
    NR = N-IM+2
    CALL NMRV2 (NL,H,QL,S,UL)
    CALL NMRV2 (NR,H,QR,S,UR)
!
! Rescale the left solution
!
    FACT = UR(NR-1)/UL(IM)
    DO I = 1, NL
          UL(I) = FACT*UL(I)
    END DO
!
    F1 = UR(NR)+UL(NL)-UR(NR-2)-UL(NL-2)
    F1 = F1/(2.0*H*UR(NR-1))
!
    IF ((F0*F1).LT.0) THEN
      EB = E 
      DE = EB-EA
    ELSE
      EA = E
      DE = EB-EA
      F0 = F1
    END IF
    ISTEP = ISTEP+1
  END DO
!
  SUM = 0.0
  DO I = 1, N
    IF(I.GT.IM) UL(I) = UR(N-I)
    SUM = SUM+UL(I)*UL(I)
  END DO
!
  WRITE(6,"(2I4)") ISTEP,IMAX
  WRITE(6,"(4F20.8)") E,DE,F0,F1
!
  SUM=SQRT(H*SUM)
  DO I = 1, N, M
    XL = XL0+(I-1)*H
    UL(I) = UL(I)/SUM
    WRITE(15,"(4F20.8)") XL,UL(I)
    WRITE(16,"(4F20.8)") XL,VX(XL)
  END DO
END PROGRAM SCHR
!
FUNCTION VX (X) RESULT (V)
  REAL :: A,B,X,V
!
  A = 1.0
  B = 4.0
  V = 3.0-A*A*B*(B-1.0)/(COSH(A*X)**2)/2.0
END FUNCTION VX
!
SUBROUTINE NMRV2 (N,H,Q,S,U)
!
! The Numerov algorithm for the equation u"(x)+q(x)u(x)=s(x)
! as given in Eqs. (3.82)-(3.85) in the book.
! Copyright (c) Tao Pang 1997.
!
  IMPLICIT NONE
  INTEGER, INTENT (IN) :: N
  INTEGER :: I
  REAL,INTENT (IN) :: H
  REAL :: G,C0,C1,C2,D,UTMP
  REAL, INTENT (IN), DIMENSION (N) :: Q,S
  REAL, INTENT (INOUT), DIMENSION (N) :: U
!
  G = H*H/12.0
!
  DO I = 2, N-1
    C0 = 1.0+G*Q(I-1)
    C1 = 2.0-10.0*G*Q(I)
    C2 = 1.0+G*Q(I+1)
    D  = G*(S(I+1)+S(I-1)+10.0*S(I))
    UTMP   = C1*U(I)-C0*U(I-1)+D
    U(I+1) = UTMP/C2
  END DO
END SUBROUTINE NMRV2