PROGRAM TDSE ! 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 models particles falling in an electric field using quantum ! mechanics. This program outputs the spread in the momentum at t = 0, and ! also provides the mean and standard deviation for the position and ! momentum. ! ! Programmed by: ! Matthew Norton ! !---------------------------------------------------------------------------- IMPLICIT NONE REAL*8, PARAMETER :: E0 = 1, Dx0 = 1.0, dX = 0.1, Pie = -1.0, TMax = 51.20 INTEGER, PARAMETER :: NDIMX = 512 INTEGER :: I, SIGN REAL*8 :: Pi, dw, Xn(NDIMX), PhiK(NDIMX), Psum = 0.0, P2 = 0.0, PStd = 0.0, K(NDIMX) REAL*8 :: XSum = 0.0, X2 = 0.0, Ps(NDIMX), Sum = (0.0,0.0) COMPLEX*8 :: PPsi(NDIMX)=(0,0), Psi(NDIMX)=(0,0) Pi = ACOS(Pie) DW = (2. * PI) / TMax OPEN (UNIT = 30, FILE = "Psi(k0).txt") DO I = 1, NDIMX Xn(I) = (- NDIMX / 2.0 - 0.5 + I) * dX Psi(I) = DEXP(-(Xn(I)/Dx0) ** 2 / 4.0) END DO CALL norm(Psi) DO I = 1, NDIMX XSum = XSum + Xn(I) * REAL(Psi(I) * CONJG(Psi(I))) X2 = X2 + Xn(I) ** 2 * REAL(Psi(I) * CONJG(Psi(I))) END DO PRINT *, "XSum", XSum, "X^2", X2, "STDEV", SQRT(ABS(X2 -(XSum * XSum))) SIGN = -1 CALL FFTBYHAND(SIGN, Psi, PPsi) Sum = (0.0,0.0) DO I = NDIMX / 2 + 1, NDIMX K(I) = (-NDIMX - 1 + I) * Dw END DO DO I = 1, NDIMX/2 k(I) = (I-1) * Dw END DO CALL norm(PPsi) DO I = 1, NDIMX ! here the negative momenta PhiK(I) = ABS(PPsi(I) * CONJG(PPsi(I))) END DO DO I = 1, NDIMX PSum = PSum + k(I) * PhiK(I) P2 = P2 + K(I) ** 2 * PhiK(I) END DO ! PSum = Psum / REAL(NDIMX) ! P2 = P2 / REAL(NDIMX) PStd = SQRT(ABS(P2 - (PSum * PSum))) PRINT *, "Psum", Psum, "P^2", P2, "STDEV", PStd ! ------------------------------------------------------------------------------------------------------- CONTAINS ! ------------------------------------------------------------------------------------------------------- ! Subroutine for slow fourier transform SUBROUTINE FFTBYHAND(SIGN, HIN, HOUT) COMPLEX*8 E1, SUM, HIN(NDIMX), HOUT(NDIMX), TPII INTEGER I, I1, SIGN TPII = 2. * Pi * (0.0,1.0)! 2PI DO I = 1, NDIMX SUM = (0.0,0.0) DO I1 = 1, NDIMX E1 = EXP(SIGN * TPII * (I1 - 1) * (I - 1.) / (NDIMX + 0.0)) SUM = SUM + E1 * HIN(I1) END DO HOUT(I) = SUM IF (SIGN == 1) HOUT(I) = SUM / (NDIMX + 0.) ! Normalization END DO END SUBROUTINE FFTBYHAND ! ------------------------------------------------------------------------------------------------------ SUBROUTINE norm(hin) complex*8 hin(512),sum INTEGER L sum=0. DO l=1,512 sum=sum+abs( hin(l))**2 END DO DO l=1,512 hin(l)=hin(l)/sqrt(sum) END DO END SUBROUTINE norm END PROGRAM