Option Explicit

Private Sub cmdCompuateProbDens_Click()

Const pi As Double = 3.1415

Const N As Long = 500

Const c As Double = 300000000# 'speed of light

Const h As Double = 6.63E-34 'Planck constant

Const m As Double = 9.1E-31 'electron mass

Const deltax As Double = 0.000001 'light wavelength

Const x0 As Double = 0

Const dx As Double = 6 * deltax / N

Const v0 As Double = 0.01 * c 'most probable speed

Const p0 As Double = v0 * m

Const deltap As Double = h / deltax

Const dp As Double = deltap / N 'momentum increment

Const xinit As Double = -2 * deltax 'left end of x range

Const pinit As Double = p0 - N / 2 * dp 'smallest momentum in range

Const psi0 As Double = 1 / N

Dim E As Double 'particle energy at given index

Dim psi(N) As Double 'resultant psi at given x

Dim prob(N) As Double 'psi squared

Dim i, j As Long 'indexes

Dim k As Double 'wave number

Dim x As Double 'x value

Dim RePsi(N), ImPsi(N) As Double

Dim p, omega, t As Double 'variables

'Pick one of the following times for plotting the prob density.

t = 0

't = 2 * deltax / v0

'Zero out the resultant psi, the probability density and the Re/Im arrays.

For i = 0 To N - 1

psi(i) = 0

prob(i) = 0

RePsi(i) = 0

ImPsi(i) = 0

Next i

'For each momentum...

For i = 0 To N - 1

p = pinit + i * dp 'current momentum

k = 2 * pi * p / h 'current wave number

E = p ^ 2 / (2 * m) 'current particle energy

omega = 2 * pi * E / h 'associated frequency

'...and for each x in the range plotted...

For j = 0 To N - 1

x = xinit + j * dx 'x value

'... compute and update the real and imaginary parts of psi

RePsi(j) = RePsi(j) + psi0 * Cos(k * x - omega * t)

ImPsi(j) = ImPsi(j) + psi0 * Sin(k * x - omega * t)

Next j

Next i

'Set psi equal to the magnitude

For i = 0 To N - 1

psi(i) = Sqr(RePsi(i) ^ 2 + ImPsi(i) ^ 2)

Next i

'Then square each resultant psi to get the probability density.

For i = 0 To N - 1

prob(i) = psi(i) ^ 2

Next i

'Output the values of prob vs. x for plotting purposes.

Open "c:\WINMCAD\Physics\ProbDens.PRN" For Output As #1

For i = 0 To N - 1

Write #1, xinit + i * dx, prob(i)

Next i

Close

MsgBox ("Ready for plotting")

End Sub