Option Explicit

Private Sub cmdSquareWell_Click()

Const pi As Double = 3.1415

Const N As Long = 500 'Number of x increments

'*****Pick the appropriate number of energy states from the list below.

Const N1 As Long = 6

'Const N1 As Long = 12

'Const N1 As Long = 1

'*****

Const h As Double = 6.63E-34 'Planck

Const m As Double = 9.11E-31 'Electron mass

Const L As Double = 0.0000000001 'Well width

Const dx As Double = L / N 'x increment

Const xinit As Double = 0 'plot from left side of well

Const psi0 As Double = 1 / (2 * N1) 'Normalization psi

Dim lambda(N1) As Double 'wavelength of wave function

Dim k(N1) As Double 'wave number

Dim E(N1) As Double 'energy

Dim omega(N1) As Double 'frequency of wave function

Dim i, j 'indexes

Dim Mult As Long 'time epoch multiplier

Dim x As Double

Dim RePsiR(N1, N), ImPsiR(N1, N) As Double 'wave traveling to right

Dim RePsiL(N1, N), ImPsiL(N1, N) As Double 'wave traveling to left

Dim Repsi(N), Impsi(N) As Double 'Sums over all energies

Dim psi(N) As Double 'magnitude of resultant psi

Dim prob(N) As Double 'psi psi*

Dim t As Double 'time epoch

'Get the time multiplier from the user.

Dim Message, Title, Default, MyValue

Message = "Enter a value between 0 and 7" ' Set prompt.

Title = "Input Multiplier" ' Set title.

Default = "0" ' Set default.

Mult = InputBox(Message, Title, Default)

'Zero out various arrays.

For i = 0 To N - 1

Repsi(i) = 0

Impsi(i) = 0

psi(i) = 0

prob(i) = 0

Next i

For i = 0 To N1 - 1

For j = 0 To N - 1

RePsiR(i, j) = 0

ImPsiR(i, j) = 0

RePsiL(i, j) = 0

ImPsiL(i, j) = 0

Next j

Next i

'Set wave function parameters.

For i = 0 To N1 - 1

lambda(i) = 2 * L / (i + 1)

k(i) = 2 * pi / lambda(i)

E(i) = (h * k(i) / (2 * pi)) ^ 2 / (2 * m)

omega(i) = 2 * pi * E(i) / h

Next i

'*****Select whole or half multiples of half periods.

t = Mult * pi / omega(0)

't = (Mult + 1 / 2) * pi / omega(0)

'*****

'Now for each value of x ...

For i = 0 To N - 1

x = xinit + i * dx

'... and for each energy at that value of x ...

For j = 0 To N1 - 1

'...compute the real and imaginary parts of the right and left

'propagating waves.

RePsiR(j, i) = psi0 * Cos(k(j) * x - omega(j) * t)

ImPsiR(j, i) = psi0 * Sin(k(j) * x - omega(j) * t)

RePsiL(j, i) = -psi0 * Cos(k(j) * x + omega(j) * t)

ImPsiL(j, i) = -psi0 * Sin(-k(j) * x - omega(j) * t)

Next j

Next i

'Then sum all the energy psis.

For i = 0 To N - 1

For j = 0 To N1 - 1

Repsi(i) = Repsi(i) + RePsiR(j, i) + RePsiL(j, i)

Impsi(i) = Impsi(i) + ImPsiR(j, i) + ImPsiL(j, i)

'MsgBox ("Repsi(" & i & ") = " & Repsi(i) & "; Impsi(" & i & ") = " & Impsi(i))

Next j

Next i

'Compute the resultant psi magnitude.

For i = 0 To N - 1

psi(i) = Sqr(Repsi(i) ^ 2 + Impsi(i) ^ 2)

Next i

'Square the magnitude to get the probability density.

For i = 0 To N - 1

prob(i) = psi(i) ^ 2

Next i

'Output the probability density 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