PROGRAM Quantum
! ------------------------------------------------------------------------------
!
! Copyright 2007 Matthew Norton
!
! ------------------------------------------------------------------------------
!    This program is free software; you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation; either version 2 of the License, or
!    (at your option) any later version.
!
!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with this program; if not, write to the Free Software
!    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
!
! NOTE: This is in FORTRAN 90 and you will need a FORTRAN 90 complier to run 
! this code
! ------------------------------------------------------------------------------
!
! This Program is designed to approximate the solution to the Schr!
! ------------------------------------------------------------------------------
!
! Created by Matthew S. Norton
! Created on March 21, 2007
! Updated on March 22, 2007
! 
! ------------------------------------------------------------------------------

IMPLICIT NONE

REAL*8 Y(2), DT = 0.0, Time = 0.0, Energy = 0.0, Hbar = 0.0, mass = 0.0, Potential = 0.0, Const = 0.0
INTEGER :: I = 0

PRINT *, "What is the mass of the particle?"
READ *, Mass
PRINT *, "What is the Potential?"
READ *, Potential
PRINT *, "What is the energy (|E|<|V| for bound state)?"
READ *, Energy

! 2m(E-V)/h^2
Const = ((2 * Mass * (Energy - Potential)) / (Hbar * Hbar))   

Hbar = 1.05457E-34

     OPEN(1, File = "Finite.txt")
	      
! Initial Conditions
     Y(1) = 5.
     Y(2) = 5.
      
! Time step
     DT = 0.1
     
     DO I = 1, 100
           Time = REAL((I - 1)) * DT
           CALL RK4(Y, Time, DT)
          WRITE(*,*) Time, Y(1), Y(2)
          WRITE(1,*) Time, Y(1), Y(2)
     END DO
STOP
END PROGRAM	
! ------------------------------------------------------------------------------
SUBROUTINE FDT(Y, Time, F)
REAL*8 Y(2), F(2), Time 
F(1) = Y(1)
F(2) = -Const * Y(2)
RETURN 
END SUBROUTINE
! ----------------------------------------------------------------------------
! y(t) --> y(t+dt) using RK$ 
! ----------------------------------------------------------------------------
SUBROUTINE RK4(Y, Time, DT)
REAL*8 YH(2), Y(2), DT, Y1DT(2), Y2DT(2), Y3DT(2), Y4DT(2), Time
INTEGER L
! Calculate K1
     CALL FDT(Y, Time, Y1DT)
! ----------------------------------------------------------------------------
! Calculate K2
     DO L = 1, 2
          YH(L)=Y(L) + 0.5 * DT * Y1DT(L) 
     END DO
     CALL FDT(YH, Time + DT / 2., Y2DT)
! ----------------------------------------------------------------------------
! Calculate K3
     DO L = 1, 2
          YH(L) = Y(L) + DT * Y2DT(L)
     END DO
     CALL FDT(YH, Time + DT / 2., Y3DT)
! ----------------------------------------------------------------------------
! Calculate k4
     DO L = 1, 2
          YH(L) = Y(L) + DT + Y3DT(L)
     END DO
     CALL FDT(YH, Time + DT, Y4DT)
! ----------------------------------------------------------------------------
! RK4
     DO L = 1, 2
          Y(L) = Y(1) + (Y1DT(L) + 2 * Y2DT(L) + 2 * Y3DT(L) + Y4DT(L)) * DT / 6.
     END DO
! ----------------------------------------------------------------------------
RETURN
END SUBROUTINE